Self-purification and entanglement revival in lambda matter

Quantum entanglement plays a pivotal role in quantum information processing, specifically in quantum communication [1, 2], quantum simulations [3], and quantum computing [4]. However, environmental noise poses a significant challenge to preserving entanglement and coherence. For the majority of systems, decoherence results in an irreversible loss of entanglement, becoming a major hindrance in advancing quantum information technologies [5]. Despite this common behavior, certain exceptional systems exhibit a fascinating phenomenon known as entanglement revival, where entanglement increases during their dissipative evolution. Moreover, a related scenario called 'entanglement sudden birth' has been explored, wherein entanglement arises from a separable state after a specific time interval. Entanglement revival and sudden birth have been extensively studied in a variety of discrete and continuous variable systems [6–20]. This non-monotonic behavior relies on a delicate interplay between entanglement generation, induced by internal system interactions, and dissipation from the external environment [6, 7, 9]. Generally speaking, the entanglement revival shows a sensitive dependence on the initial state [12–14, 17, 18]. In some systems, it was suggested that revivals might be induced by the transfer of entanglement between different subsystems [21–23].

In this work, we present a comprehensive investigation of a robust mechanism for entanglement revival, centered around the utilization of special metastable states. These unique states act as entanglement reservoirs during the later stages of time evolution and can be deliberately engineered by leveraging specific (exact or approximate) symmetries inherent in the underlying quantum system. Our focus lies on an ensemble of qutrits exhibiting a lambda-type level structure interacting with a single bosonic mode, a model pertinent to various platforms such as cavity quantum electrodynamics (QED) [24, 25], trapped ions [26], circuit QED [27], and atomic ensembles [28]. The dissipative evolution in our study follows a universal two-stage dynamics. In the initial stage, we observe the expected rapid decay of entanglement at short times. Subsequently, in the second stage leading to entanglement revival, the system undergoes a slower evolution caused by relaxation, ultimately converging to a special dark state. This self-purification process plays a crucial role in restoring entanglement. Importantly, the entire process is closely connected to an antisymmetry (rather than symmetry) associated with the parity operator of the system. Furthermore, we provide a comprehensive description of the quantum states involved, enabling a clear identification of the underlying physical processes at play. Armed with this in-depth understanding, we anticipate that the essential features of the revival dynamics can be realized in a broader class of systems.

This paper is structured as follows: In section 2, we introduce the system, detailing its fundamental physical processes and associated properties. Section 3 explores the observation of entanglement revival behavior, illustrating the dissipative evolution characterized by universal two-stage dynamics. Section 4 provides a semiclassical perspective on the system, while sections 5 and 6 elaborate on the detailed decay processes based on this semiclassical description. Finally, section 7 summarizes the findings and presents an outlook for future research.

We consider n identical three-level systems (qutrits, for short) coupled with single-mode bosons. As shown in figure 1, each qutrit has a Λ-type structure with two ground-state levels, denoted by $|0\rangle$ and $|2\rangle$, and one excited-state level, $|1\rangle$. The bosonic mode induces the transition $|{0}\rangle\leftrightarrow|1\rangle$ whereas the transition $|{1}\rangle\leftrightarrow|{2}\rangle$ is driven resonantly by an external classical field. In the interaction picture, the system is governed by the Hamiltonian

$$\begin{align} \hat{H} = g\sum_{k = 1}^n\hat{c}\,|1\rangle_k\langle0| + \Omega\sum_{k = 1}^n|1\rangle_k\langle2| + \mathrm{h.c.} , \end{align} \tag{ 1 }$$

where $\hat{c}$ is the bosonic annihilation operator, $|\mu\rangle_k$ ($\mu = 0,1,2$) are the quantum states of the kth qutrit, g is the qutrit-boson coupling, and Ω is the Rabi transition amplitude. Note that, in this work, we assume uniform g and Ω for all qutrits. The operator $\hat{N} = \hat{c}^\dagger \hat{c} +\sum_{k = 1}^n(|1\rangle_k\langle 1| + |2\rangle_k\langle 2|)$ is conserved by $\hat H$, and we will usually refer to its eigenvalues (i.e. the total number of excitations) as p.

Zoom In Zoom Out Reset image size

Any realistic system is subject to the influence of quantum noise, thus its dynamics is not unitary any longer. Dissipative dynamics is typically modeled by a quantum master equation, which we take of the following form

$$\begin{align} \dot{\hat\rho}\left(t\right) = -i\left[\hat{H},\hat\rho\right] + \kappa\mathcal{L}\left[\hat{c}\right]\hat\rho + \Gamma_0\mathcal{L}\left[\hat{L}_0\right]\hat\rho + \Gamma_2\mathcal{L}\left[\hat{L}_2\right]\hat\rho , \end{align} \tag{ 2 }$$

where $\hat\rho$ is the density operator of the system and $\mathcal{L}[\hat L]\hat\rho$ is the Lindblad superoperator, defined by $\mathcal{L}[\hat L]\hat\rho: = \hat L\hat\rho\hat L^\dagger-(\hat L^\dagger\hat L\hat\rho+\hat \rho\hat L^\dagger\hat L)/2$ for the linear operator $\hat{L}$ associated with a quantum decoherence process. Specifically, $\mathcal{L}[\hat{c}]$ is responsible for the loss of bosons while $\mathcal{L}[\hat{L}_\mu]$, with $\hat{L}_\mu : = \sum_{k = 1}^n{|{\mu}\rangle_k}{\langle{1}|},$ describes the collective spontaneous decay of qutrits from the excited-state level $|{1}\rangle$ to the ground-state level $|\mu\rangle$ ($\mu = 0,2$); κ and $\Gamma_\mu$ are the corresponding rates. As implied by the specific form of the quantum master equation (2), we mainly concentrate on collective decay of qutrits. However, as we discuss briefly later and demonstrate in detail in the appendix B, our main results are not affected qualitatively by individual qutrit decay.

In this work, we investigate an unusual dynamics of entanglement in the presence of quantum decoherence. The effect of our interest is most pronounced in the zero-energy subspace of $\hat H$ with p excitations and totally symmetric wave functions (with respect to the exchange of qutrits). Therefore we mainly focus on this subspace, denoted as $\mathcal{Z}_p$. Interestingly, it was pointed out [29] that the $\mathcal{Z}_p$ subspace features a decoherence-free nature [30–32], and is always degenerate regardless of the parameters, allowing one to geometrically manipulate quantum states within the subspace. This property is due to three symmetry properties of the system: the conservation of total number of excitations $\hat N$, the exchange symmetry of qutrits, and the anti-symmetry, $\{\hat\Pi_1,\hat{H}\} = 0,$ of $\hat\Pi_1 : = \exp\left(i\pi\sum_{k = 1}^n|1\rangle_k\langle1|\right).$

We denote the basis states of $\mathcal{Z}_p$ by $\mathinner{\lvert{\textstyle Z_p^i}\rangle}$ for $i = 0,1,2,\ldots,[p/2]$ ⁴ , where $[x]$ indicates the integer part of the real number x. The first zero-energy eigenstate can always be chosen of the following form:

$$\begin{align} |Z\,_{p}^0\rangle = \sum_{k = 0}^{p} \frac{\left(-1\right)^{p-k}}{\sqrt{\left(p-k\right)!}g^{p-k} \Omega^{k}}|\Phi_n^k\rangle_Q|p-k\rangle_c, \end{align} \tag{ 3 }$$

where $\mathinner{\lvert{\textstyle\Phi_n^k}\rangle}_Q = \sum_\mathcal{P}\mathcal{P}\mathinner{\lvert{\textstyle2}\rangle}^{\otimes k}\mathinner{\lvert{\textstyle0}\rangle}^{\otimes(n-k)}$ ($\mathcal{P}$ are permutations of qutrits) is a symmetric Dicke state and $\mathinner{\lvert{\textstyle k}\rangle}_c$ is a Fock state with k bosons. Equation (3) describes a special zero-energy state which we call master dark state. It has several interesting properties useful for quantum-state engineering applications, including the generation of arbitrary symmetric Dicke states [29, 33]. We will see below that this master dark state plays a key role in the entanglement-revival behavior as well.

Figure 2 illustrates the main phenomenon in the focus of this work. For four qutrits (n = 4) with three excitations (p = 3), which we take as a prototypical example, there are two zero-energy states, $\mathinner{\lvert{\textstyle Z\,_p^0}\rangle}$ and $\mathinner{\lvert{\textstyle Z_p^1}\rangle}$, where $\mathinner{\lvert{\textstyle Z_p^1}\rangle}$ has a relatively large qutrit-boson entanglement. In this work, we quantify the entanglement by the logarithmic negativity [34–36], defined as

$$\begin{align} E_N\left(\rho_{AB}\right): = \log_2 \left\|\rho_{AB}^{T_B} \right\|, \end{align} \tag{ 4 }$$

where $\|\cdots\|$ denotes the trace norm and $\rho^{T_B}$ indicates the partial transposition of ρ with respect to subsystem B. From the quantum dissipative dynamics of equation (2), we find in figure 2(a) that the entanglement content of $\mathinner{\lvert{\textstyle Z_p^1}\rangle}$ decreases fast initially, as usually expected. Surprisingly, however, the entanglement increases again to a certain level as time evolves. This type of two-stage dynamics and entanglement-revival behavior are universally found in a wide range of parameters.

Zoom In Zoom Out Reset image size

In search for clues to explain the above behavior, we consider the probabilities $P_p^k(t) : = \mathinner{\langle{\textstyle Z_p^k}\rvert}\hat\rho(t)\mathinner{\lvert{\textstyle Z_p^k}\rangle}$ of being in a zero-energy state ($p\unicode{x2A7D} 3$). We find that, among all zero-energy states, the key role is played by $\mathinner{\lvert{\textstyle Z_3^1}\rangle}$, $\mathinner{\lvert{\textstyle Z\,_1^0}\rangle}$, and $\mathinner{\lvert{\textstyle Z\,_0^0}\rangle}$, whose populations are shown in in figure 2(b). In the first stage of the dynamics, the initial state $|Z_3^1\rangle$ decays rapidly (over a time scale of order $1/\kappa$ or $1/\Gamma_\mu$) to the p = 1 master dark state $|Z\,_1^0\rangle$. Once the master dark state $|Z\,_1^0\rangle$ becomes dominant, the second stage of dynamics kicks in, and $\mathinner{\lvert{\textstyle Z\,_1^0}\rangle}$ decays to the trivial state $\mathinner{\lvert{\textstyle Z\,_0^0}\rangle}\equiv\mathinner{\lvert{\textstyle0}\rangle}_1\otimes\cdots\otimes\mathinner{\lvert{\textstyle0}\rangle}_n\otimes\mathinner{\lvert{0}\rangle}_c.$ The latter process is much slower, since the master dark states are extremely robust to decoherence [29]. As inferred from equation (3), they are not affected by spontaneous decay of qutrits (neither individual or collective) from the excited-state to ground-state levels, and become immune to boson losses in the limit of large $g/\Omega$. Between the two stages there is a certain transition period (shaded region in figure 2) characterized by a low purity, see the dashed curve in figure 2(b). Importantly, the probabilities $P_3^1, P_1^0, P_0^0$ do not sum up to unity during this transition period. This means that some other states give a significant contribution to $\hat\rho(t)$. Indeed, as we discuss in detail below, the decay of the initial state $\mathinner{\lvert{\textstyle Z_3^1}\rangle}$ to $\mathinner{\lvert{\textstyle Z\,_1^0}\rangle}$ occurs via intermediate states with finite energies, accounting for the missing population $1-P_3^1-P_1^0-P_0^0$. The important point here is that the mixture of those states dramatically suppresses the overall entanglement to a level that is much smaller than the entanglement of each component (pure) state.

We conclude that the entanglement revival can be attributed to an initial suppression of entanglement due to strong mixture, followed by a recovery of entanglement when the system evolves towards the (relatively) pure master dark state which, like the initial state, is significantly entangled. This self-purification is possible because the master dark state is stable against the decay processes. As noted already, the effect of boson loss gets suppressed by increasing $g/\Omega$. However, this also leads to a master dark state with smaller entanglement, implying that the visibility of the entanglement revival can be optimized over the coupling ratio in appendix A.

We now quantitatively analyze the two-stage dissipative dynamics and the corresponding entanglement-revival behavior, observed numerically above. To do this, we first simplify our model through a semiclassical approximation, recalling that we are working in the subspace of totally symmetric wave functions. As long as the quantum master equation (2) contains only collective decoherence of qutrits and boson loss, the system remains in this subspace. Then, it is convenient to describe the ensemble of qutrits using bosonic operators $\hat{a}_\mu$ associated with the levels $\mathinner{\lvert{\textstyle\mu}\rangle}$. Expressed in terms of these bosonic operators, $\hat{H}$ reads

$$\begin{align} \hat{H} = g\hat{c}^\dagger \hat a_0^\dagger\hat a_1 + \Omega\hat{a}_2^\dagger\hat{a}_1 + \mathrm{h.c.}, \end{align} \tag{ 5 }$$

which preserves the total number of qutrits $\hat{a}_0^\dagger\hat{a}_0 + \hat{a}_1^\dagger\hat{a}_1 + \hat{a}_2^\dagger\hat{a}_2 = n$ as well as the total number of excitations $\hat{a}_1^\dagger\hat{a}_1 + \hat{a}_2^\dagger\hat{a}_2 + \hat{c}^\dagger\hat{c} = p.$ Likewise, equation (2) becomes:

$$\begin{align} \dot{\hat\rho}\left(t\right) = -i\left[\hat{H},\hat\rho\right] + \kappa\mathcal{L}\left[\hat{c}\right]\hat\rho + \Gamma_0\mathcal{L}\left[\hat{a}_0^\dagger\hat{a}_1\right]\hat\rho + \Gamma_2\mathcal{L}\left[\hat{a}_2^\dagger\hat{a}_1\right]\hat\rho.\end{align} \tag{ 6 }$$

So far, everything is exact. We now make a semiclassical approximation $g\hat{c}^\dagger\hat a_0^\dagger\hat a_1 \approx g\hat{c}^\dagger \hat a_1\sqrt{n},$ which is valid in the limit $n\gg p$. Then, equation (5) gives:

$$\begin{align} \hat{H} \approx \hat{H}_\mathrm{sc} : = g\sqrt{n}\,\hat{c}^\dagger\hat{a}_1 + \Omega\,\hat{a}_2^\dagger\hat{a}_1 + \mathrm{h.c.} \end{align} \tag{ 7 }$$

This semiclassical Hamiltonian is quadratic in the bosonic operators, hence can be solved exactly. By introducing three new bosonic operators

$$\begin{align} \hat{C}_{0} &: = \cos\theta\,\hat{c} - \sin\theta\,\hat{a}_{2}, \end{align} \tag{ 8a }$$

$$\begin{align} \hat{C}_{\pm} &: = \left(\sin\theta\,\hat{c}+\cos\theta\,\hat{a}_{2} \mp \hat{a}_{1}\right)/\sqrt2, \end{align} \tag{ 8b }$$

where $\tan\theta: = g\sqrt{n}/\Omega$, the semiclassical Hamiltonian takes the following form:

$$\begin{align} \hat{H}_\mathrm{sc} = \epsilon\left(\hat{C}_+^\dagger \hat{C}_+-\hat{C}_-^\dagger\hat{C}_-\right), \end{align} \tag{ 9 }$$

where $\epsilon: = \sqrt{g^2n+\Omega^2}$. We note from equation (9) that $\hat{C}_0$ is a zero-frequency bosonic eigenmode, while $\hat{C}_\pm$ have opposite mode-frequencies ±ε. This interesting property is due to the anti-symmetry $\{\hat\Pi_1,\hat{H}\} = 0,$ with $\hat\Pi_1$ expressed now as $\hat\Pi_1 = \exp\left(i\pi\hat{a}_1^\dagger\hat{a}_1\right).$ The semiclassical Hamiltonian immediately yields the eigenenergies $(k_+-k_-)\epsilon$, with the corresponding eigenstates:

$$\begin{align} |E_{k_{0};k_{+}k_{-}}\rangle = \frac{1}{\sqrt{k_0! k_+! k_-!}} \left(\hat{C}_{0}^{\dagger}\right)^{k_{0}} \left(\hat{C}_{+}^{\dagger}\right)^{k_{+}} \left(\hat{C}_{-}^{\dagger}\right)^{k_{-}}\mathinner{\lvert{\textstyle\;}\rangle}, \end{align} \tag{ 10 }$$

where $\mathinner{\lvert{\textstyle\;}\rangle}$ is the vacuum state. Of particular importance are the semiclassical zero-energy states with p excitations (corresponding to the subspace $\mathcal{Z}_p$), which take the general form $\mathinner{\lvert{\textstyle E_{(p-2k);kk}}\rangle}.$ The case k = 0 gives the master dark states of equation (3), i.e. $\mathinner{\lvert{\textstyle E_{p;00}}\rangle} \approx \mathinner{\lvert{\textstyle Z\,_p^0}\rangle}$. For p = 3, applicable to figure 2, the two semiclassical zero-energy states are $|E_{3;00}\rangle\approx|Z_3^0\rangle$ and $|E_{1;11}\rangle\approx|Z_3^1\rangle$.

In the physically relevant regime of $\kappa,\Gamma_\mu\lesssim g,\Omega$, the relative phase between any pair of semiclassical eigenstates with different energies oscillates quickly (the energy difference is of order $\epsilon\gg\kappa,\Gamma_\mu$). Therefore, by ignoring the phase coherence between such eigenstates, we consider a semiclassical solution of the following form

$$\begin{align} \hat \rho_\mathrm{sc}\left(t\right) \approx \sum_\alpha P_\alpha \left(t\right) |E_\alpha\rangle \langle E_\alpha|, \end{align} \tag{ 11 }$$

where $\alpha\equiv(k_0;k_+k_-)$ collectively denotes the semiclassical quantum numbers k₀ and $k_\pm$. Substituting this ansatz into equation (6), we obtain the (classical) equations

$$\begin{align} \frac{{\mathrm {d}}P_\alpha}{{\mathrm {d}}t} = \sum_{\beta}\gamma_{\alpha\beta}P_\beta\left(t\right) \end{align} \tag{ 12 }$$

for the probabilities $P_\alpha(t)$, where $\gamma_{\alpha\beta}$ is the transition rate between semiclassical states $|\alpha\rangle$ and $|\beta\rangle$.

We now apply the semiclassical description of equations (11) and (12) to analyze in detail the dynamical process. We first consider the effects of boson loss and, for simplicity, we focus on the particular example with p = 3. Similar results apply to arbitrary p.

In the semiclassical limit ($n\gg p$, hence $\tan\theta\gg 1$), the rates $\gamma_{\alpha\beta} : = \kappa |\langle E_\alpha|\hat{c}|E_{\beta}\rangle|^2$ of the four processes $\mathinner{\lvert{\textstyle E_{1;11}}\rangle} \to \mathinner{\lvert{\textstyle E_{1;10}}\rangle}$, $\mathinner{\lvert{\textstyle E_{1;11}}\rangle} \to \mathinner{\lvert{\textstyle E_{1;01}}\rangle}$, $\mathinner{\lvert{\textstyle E_{1;10}}\rangle} \to \mathinner{\lvert{\textstyle E_{1;00}}\rangle}$, and $\mathinner{\lvert{\textstyle E_{1;01}}\rangle} \to \mathinner{\lvert{\textstyle E_{1;00}}\rangle},$ are much larger (by a factor of order $g^2n/\Omega^2\gg 1$) then other transition rates. Therefore, as illustrated in figure 3, one can identify two dominant (incoherent) decay paths of $\mathinner{\lvert{\textstyle Z_3^1}\rangle}\approx\mathinner{\lvert{\textstyle E_{1;11}}\rangle}$:

$$\begin{align} & \mathinner{\lvert{\textstyle E_{1;11}}\rangle} \to \mathinner{\lvert{\textstyle E_{1;10}}\rangle} \to \mathinner{\lvert{\textstyle E_{1;00}}\rangle}, \end{align} \tag{ 13 }$$

$$\begin{align} \notag\\[-6pt] & \mathinner{\lvert{\textstyle E_{1;11}}\rangle} \to \mathinner{\lvert{\textstyle E_{1;01}}\rangle} \to \mathinner{\lvert{\textstyle E_{1;00}}\rangle}. \end{align} \tag{ 14 }$$

We also include $\mathinner{\lvert{\textstyle E_{1;00}}\rangle} \to \mathinner{\lvert{\textstyle E_{0;00}}\rangle}$ among the major transitions of figure 3(solid arrows). The latter process has a much smaller rate, but is the only allowed transition once the system has reached $\mathinner{\lvert{\textstyle E_{1;00}}\rangle}$.

Zoom In Zoom Out Reset image size

The above remarks allow us to further simplify our semiclassical description, by only including in the ansatz $\hat\rho_\mathrm{sc}(t)$ of equation (11) the five most relevant states $\mathinner{\lvert{\textstyle E_{1;11}}\rangle}$, $\mathinner{\lvert{\textstyle E_{1;10}}\rangle}$, $\mathinner{\lvert{\textstyle E_{1;01}}\rangle}$, $\mathinner{\lvert{\textstyle E_{1;00}}\rangle}$, and $\mathinner{\lvert{\textstyle E_{0;00}}\rangle}$, which appear in the major transition paths. As figure 4(a) demonstrates, the approximate evolution of the five populations $P_\alpha(t)$ agrees very well with the full numerical solution, obtained from equation (6). In figure 4(b) we further compare the entanglement content of the dominant semiclassical state $\hat\rho_\mathrm{sc}(t)$ (red dashed line) to the full solution $\hat\rho(t)$ (blue empty circles). Here the agreement is worse than for the probabilities of panel (a); apparently, the entanglement content is more sensitive to the detailed form of $\hat\rho(t)$. However, the discrepancy can be easily corrected by including all semiclassical eigenstates $\mathinner{\lvert{\textstyle E_{k_0;k_+k_-}}\rangle}$ in $\hat\rho_\mathrm{sc}(t)$. In this case, the possible transition paths include those marked by dashed arrows in figure 3. The black dashed curve of figure 4(b) plots the evolution of the entanglement content within the full semiclassical description, and shows an excellent agreement with the full solution.

Zoom In Zoom Out Reset image size

We have shown above that, as expected, the semiclassical approximation works very well when $n\gg p$. However, it also gives qualitatively correct predictions of the dynamical behavior for values of n as small as n = 4 (with p = 3). Many features can be even described quantitatively, as long as we replace $\hat{a}^\dagger_0\hat{a}_0$ by the average occupation n₀ of state $|0\rangle$.

The semiclassical description is still applicable when collective decay of qutrits is the dominant decoherence mechanism. Assuming $\kappa = \Gamma_2 = 0$, we illustrate in figure 5(a) the transitions induced by the $\hat{a}_0^\dagger \hat{a}_1 \simeq \sqrt{n} \hat{a}_1$ process of equation (2). This semiclassical diagram is similar to the one in figure 3, except that now mode $\hat{C}_0$ is not affected by dissipation. Therefore, all the slow transitions of figure 3 do not take place (including $|E_{1;00}\rangle \to |E_{0;00}\rangle$). The system will ultimately reach the master dark state $|E_{1;00}\rangle$, instead of the vacuum state. In figure 5(a), the transition rates are all of order $\Gamma = \sqrt{n} \Gamma_0$, i.e. collectively enhanced by the occupancy of state $|0\rangle$. If, on the other hand, the $\hat{a}_2^\dagger \hat{a}_1$ process dominates ($\kappa = \Gamma_0 = 0$) we find the complete semiclassical diagram shown in figure 5(b). Since the decay from $|1\rangle$ to $|2\rangle$ does not alter the excitation number, the dynamics is confined to the p = 3 subspace, and terminates in the master dark state $|E_{3;00}\rangle$.

Zoom In Zoom Out Reset image size

As before, both diagrams of figure 5 lead naturally to a two-stage dynamics where, after an initial loss of purity, the system evolves towards a master dark state. As a consequence, revivals are observed in the time dependence of the logarithmic negativity, shown in figure 6. Since the master dark states are completely immune to the decay of qutrits, at long times the entanglement content saturates to a finite value.

Zoom In Zoom Out Reset image size

When both types of qutrit decay are present, the semiclassical diagram is more complex but can be obtained in a similar manner. An entanglement revival is still found with the parameters of figure 6 (see the lowest curve). However, now the system relaxes to a mixture of three master dark states, $|E_{3;00}\rangle$, $|E_{1;00}\rangle$, and $|E_{2;00}\rangle$, leading to a general reduction of entanglement. As seen, in all three cases of figure 6 the semiclassical approximation (dashed curves) is in close agreement with the exact evolution (empty circles).

Finally, we comment on the effect of individual decay of qutrits. This process is more complex to describe, as it causes transitions between different symmetry sectors, expanding the dynamical process beyond the totally symmetric subspace. Nevertheless, we find that the dynamics is qualitatively similar in appendix B. The system can still reach a distinct master dark state, causing the entanglement to revive. Only when the decay rates from $|1\rangle$ to $|0\rangle$ and from $|1\rangle$ to $|2\rangle$ are comparable, a larger number of states is involved and this might prevent the self-purification process in appendix B. The analysis of individual qubit decay suggests that the revival is robust to other types of local perturbations, e.g. small deviations of $g,\Omega$ from the homogenous limit. The entanglement revival also survives in the presence of detuning in appendix C.

By analyzing the evolution of an ensemble of qutrits interacting with a bosonic mode, we have identified a robust physical mechanism for entanglement revivals. We show that starting from a highly-entangled state, decoherence leads to a universal two-stage dynamics which can be accurately captured by a semiclassical approximation. The entanglement revival is attributed to the self-purification of the quantum state as it relaxes towards a special dark state. While the specific model of our study has direct applicability to various experimental setups [24–28], the revival mechanism boasts broad applicability and bears some analogy to other extensively studied protocols to generate entanglement through dissipation such as reservoir engineering [37–41]. However, our primary focus here is on inducing an unconventional form of quantum evolution.

S C acknowledges support from the Innovation Program for Quantum Science and Technology (Grant No. 2021ZD0301602), the National Science Association Funds (Grant No. U2230402), and the National Natural Science Foundation of China (Grant Nos. 11974040 and 12150610464). C D and M S C has been supported by the National Research Function (NRF) of Korea (Grant Nos. 2022M3H3A106307411 and 2023R1A2C1005588) and by the Ministry of Education through the BK21 Four program.

All data that support the findings of this study are included within the article (and any supplementary files).

Here, we analyze the optimal value of the coupling ratio $\Omega/(g\sqrt{n})$, which is determined by the competing requirements of having a dark state with large entanglement and being weakly affected by boson decay.

For simplicity and to focus on the essential point, we assume no qutrit decay (the effect of which does not affect the qualitative feature). In this case, the semiclassical dynamics is described by figure 3 of the main text. We also recall that the revival of entanglement is due to a self-purification of the system, as it relaxes to the special master dark state $|E_{1;00}\rangle$. This process takes place on a timescale $\sim\kappa^{-1}$. Later, boson decay causes a slow leakage towards the vacuum state, $|E_{1;00}\rangle \to |E_{0;00}\rangle$. The rate of this process is of order $\kappa/\tan\theta$ and can be suppressed by making $\tan\theta$ as large as possible. In this limit the probability of the master dark state can approach one (ideal purification), but the entanglement between the qutrits and the cavity becomes vanishingly small. We see from equation (3) of the main text that $|Z\,_p^0\rangle \simeq |\Phi_n^p\rangle_Q|0\rangle_c$ when $\tan\theta = g\sqrt{n}/\Omega \gg 1$.

Because of this competition, the entanglement revival will at first become more visible by increasing the ratio of the two coupling strengths $1/\tan\theta = \Omega/(g\sqrt{n})$ but gradually disappears at larger values. This behavior is illustrated in panel (a) of figure A1, where we plot the time dependence of the entanglement with various coupling strength ratios. The height of the revival as function of $1/\tan\theta$ is plotted in figure A1(b), showing that the optimal ratio is in the range of ${\textrm{0}}{\textrm{.05 - - 0}}{\textrm{.1}}$.

Zoom In Zoom Out Reset image size

This optimal ratio can be estimated through simple analytical calculations. At long times $t \gg \kappa^{-1}$, the initial state $|E_{1;11}\rangle$ has almost fully decayed to a mixture of $|E_{1;00}\rangle$ and $|E_{0;00}\rangle$. The approximate density matrix, denoted by $\tilde{\rho} (t)$, can be expressed in the following form:

$$\begin{align} \tilde{\rho}\left(t\right) = P_{1;00}\left(t\right)|E_{1;00}\rangle\langle E_{1;00}|+\left(1-P_{1;00}\left(t\right)\right)|E_{0;00}\rangle\langle E_{0;00}|, \end{align} \tag{ A.1 }$$

where the probability can be derived through the semiclassical rate equation, equation (11) of the main text, giving

$$\begin{align} P_{1;00}\left(t\right) = e^{-\kappa t}\left(1-e^{\frac{\sin\theta}{2\left(\sin\theta+\cos\theta\right)}\kappa t}\right)^2. \end{align} \tag{ A.2 }$$

As a reference, we consider the typical timescale $\kappa t = 10$ and compute the entanglement of the approximate density operator $\tilde{\rho}(10)$. The entanglement obtained in this manner is plotted in figure A2, from where we estimate an optimal coupling ratio $1/\tan\theta\approx 0.13$, which is close to the value from the full simulations of figure A1.

Zoom In Zoom Out Reset image size

To incorporate the effect of individual decay of qutrits, we extend the master equation of the main text as follows:

$$\begin{align} \dot{\hat{\rho}} = -i \left[\hat{H},\hat{\rho}\right] +\kappa\mathcal{L}\left[\hat{c}\right]\hat{\rho} +\sum_{\mu = 0,2} \left(\Gamma_{\mu}\mathcal{L}\left[\hat L_\mu\right]\hat{\rho}+ \sum_{j = 1}^n \Gamma_{1\mu}\mathcal{L}\left[\hat{A}^j_{1\mu}\right]\hat{\rho} \right), \end{align} \tag{ B.1 }$$

where $\hat L_\mu: = \sum_{k = 1}^n|\mu\rangle_k\langle1|$ describes the collective decay of qutrits from state $|1\rangle$ to state $|\mu\rangle$ and $\hat{A}_{1\mu}^j = |\mu\rangle_j\langle 1|$ describes the individual decay of qutrits. Due to this additional dissipation channel, the dynamics become more complex. It extends beyond the totally symmetric subspace, thus the semiclassical description is not applicable anymore. Despite these technical difficulties, the qualitative behavior is remarkably similar to the evolution with only bosonic and/or collective qutrit decay.

We first only include the individual decay of qutrits ($\kappa = \Gamma_0 = \Gamma_2 = 0$). Starting from $|E_{1;11}\rangle$, the allowed transitions between different symmetry sectors (Young diagrams) are represented in figure B1(a). The individual decay $|1\rangle\to|0\rangle$ gives the 'horizontal' transitions, which reduce the number of excitations p. On the other hand, the $|1\rangle\to|2\rangle$ decay induces 'vertical' transitions, which do not change p. If we restrict ourselves to the former process (setting $\Gamma_{12} = 0$), after two decays the systems is one of the $|E_1^m\rangle$ eigenstates of the p = 1 subspace. We can estimate the probabilities of these states as follows:

$$\begin{align} P_m \simeq \sum_k \frac{\Gamma\left(E_2^k \to E_1^m\right)}{\sum_{m^{^{\prime}}} \Gamma\left(E_2^k \to E_1^{m^{^{\prime}}}\right)}\frac{\Gamma\left(E_{1;11} \to E_2^{k}\right)}{\sum_{k^{^{\prime}}} \Gamma\left(E_{1;11} \to E_2^{k^{^{\prime}}}\right)}, \end{align} \tag{ B.2 }$$

where $\Gamma(\alpha \to \beta) = \sum_i \Gamma_{10} |\langle \beta|\hat A_{10}^i|\alpha \rangle|^2$ and the sum over k runs over all the intermediate states with p = 2. In figure B1(b) we show that the largest P_m is for the state $|E_{1;00}\rangle$, indicating that the decay to the master dark state is the dominant process. Thus, the same type of two-stage dynamics discussed already should take place.

Zoom In Zoom Out Reset image size

In agreement with this argument, we see a clear revival in figure B2(a), where the decay $|1\rangle \to |0\rangle$ is dominant. If the $|1\rangle \to |2\rangle$ decay is dominant, the dynamics is also restricted, as it remains approximately confined to the p = 3 subspace. As seen figure B2(b), a strong revival exists in this case. Instead, when the effects of $\hat{A}_{10}^j$ and $\hat{A}_{12}^j$ are comparable, all the processes shown in figure B1(a) occur on a similar timescale and it is difficult to achieve the self-purification process. As a consequence, the revival is not observed in panel (c) of figure B2.

Zoom In Zoom Out Reset image size

Lastly, we investigate the entanglement revival when both bosonic decay and individual decay of the qutrits are present. As shown in figure B3(a), where all the decay rates are comparable, the entanglement revival can persists in this scenario as well.

Zoom In Zoom Out Reset image size

Until now, we have treated the resonant case. Considering a finite detuning Δ, the effective Hamiltonian in the bosonic representation reads:

$$\begin{align} \hat{H}_{\Delta} = \Delta \hat{a}_1^\dagger\hat{a}+g\sqrt{n}\left(\hat{a}_1^\dagger\hat{a}_0+\hat{a}_0^\dagger\hat{a}_1\right)+\Omega\left(\hat{a}_1^\dagger\hat{a}_2+\hat{a}_2\hat{a}_1^\dagger\right). \end{align} \tag{ C.1 }$$

Simulations for different detunings are presented in figure B3(b). If the detuning is small compared to the effective coupling strength $g\sqrt n$, the entanglement revival behavior survives because the master dark state $|Z\,_p^0\rangle$ still satisfies $\hat{H}_{\Delta}|Z\,_p^0\rangle = 0$. If the detuning is too large, the oscillations caused by detuning and the possibility of transitions from a master dark state $|Z\,_p^0\rangle$ to other excited states can destroy the entanglement revival behavior.