Generating Entangled Steady States in Multistable Open Quantum Systems via Initial State Control

TL;DR

This work addresses how initial-state preparation controls entangled steady states in open quantum systems described by the Lindblad equation. By decomposing the steady-state manifold into kernel vectors of the Liouvillian and introducing a weight that depends on the initial state and a non-Euclidean metric, the authors provide analytic formulas to predict ρ_{SS} without simulating dynamics. They highlight a special Hermitian case where steady states are determined solely by kernel overlaps, and apply the framework to spin ensembles to design metrologically useful entangled steady states using balanced collective decay. The results offer a computationally efficient route to optimize initial states for robust long-time quantum resources, with practical implications for differential sensing and quantum metrology in cavity QED settings.

Abstract

Entanglement underpins the power of quantum technologies, yet it is fragile and typically destroyed by dissipation. Paradoxically, the same dissipation, when carefully engineered, can drive a system toward robust entangled steady states. However, this engineering task is nontrivial, as dissipative many-body systems are complex, particularly when they support multiple steady states. Here, we derive analytic expressions that predict how the steady state of a system evolving under a Lindblad equation depends on the initial state, without requiring integration of the dynamics. These results extend the frameworks developed in Refs. [Phys. Rev. A 89, 022118 (2014) and Phys. Rev. X 6, 041031 (2016)], showing that while the steady-state manifold is determined by the Liouvillian kernel, the weights within it depend on both the Liouvillian and the initial state. We identify a special class of Liouvillians for which the steady state depends only on the initial overlap with the kernel. Our framework provides analytical insight and a computationally efficient tool for predicting steady states in open quantum systems. As an application, we propose schemes to generate metrologically useful entangled steady states in spin ensembles via balanced collective decay.

Figures (5)

• Figure 1: Role of initial states in multistable systems:(a) In the case of a Hermitian Liouvillian $\mathcal{L}^{\dagger} = \mathcal{L}$, the steady state of the system $\bm{\rho}_{\rm SS}$ is completely determined by the overlap between the initial state $\bm{\rho}_0$ and the elements of the kernel of the Liouvillian [see Eq. \ref{['result:specialcase']}]. For this and panel (c), the yellow region corresponds to all the information contained in the Liouvillian. A smaller set of this area corresponds to the kernel, spanned by all vectors $\bm{\rho}$ that follow $\mathcal{L} \bm{\rho}=\bm{0}$, denoted by the purple region. The initial state is represented by the dark blue region, while the white region denotes the steady state. (b) For the case described in (a), the steady state is simply the sum of the projections of the initial state $\bm{\rho}_0$ into each of the vectors $\bm{\rho}_{{\rm SS},j}$ spanning the kernel. We note that in this case, since $T$ is unitary, the vectors $\bm{\rho}_{{\rm SS},j}$ can always be constructed to be orthogonal to each other. Moreover, in all examples studied here, these vectors correspond to valid density matrices, up to a constant factor, as seen in Eq. \ref{['vectorsBalanced']}, for example. (c) For a general Liouvillian, the steady state is not fully determined by the overlap of the initial state and the elements of the kernel (enclosed with a dashed line), and additional information from the Liouvillian is required to predict the steady state. This information is encoded in the transformation $T$ for Eq. \ref{['eq:result1']} and in the elements of the kernel of $\mathcal{L}^{\dagger}$ ($\bm{\mathcal{U}}_j$) for Eq. \ref{['eq:Albert']}. (d) For the more general case described in (c), the steady state is still a weighted sum of the vectors $\bm{\rho}_{{\rm SS},j}$ [Eq. \ref{['eq:result1']}], however, the weights $c_j$ are now defined by a different inner product $\bm{\rho}_{{\rm SS},j}^{\dagger} \mathcal{M}\, \bm{\rho}_0$, with $\mathcal{M}=(T^{-1})^\dagger T^{-1}$. We can think of $\mathcal{M}$ being a non-Euclidean metric that defines the inner product of vectors $\bm{\rho}_0$ and $\bm{\rho}_{{\rm SS},j}$. The non-Euclidean behavior is illustrated by the projections not being orthogonal (see dashed lines). For this case, the vectors $\bm{\rho}_{{\rm SS},j}$ are not necessarily orthogonal nor valid density matrices.
• Figure 2: Two qubits under collective dissipation:(a) Two qubits experience collective dissipation through jump operators $L_1=S_-=\sigma_-^1+\sigma_-^2$ and $L_2 =S_+ =\sigma_+^1+\sigma_+^2$. The two processes have the same rates $\gamma_1=\gamma_2=\gamma$. (b) Concurrence in the steady state as a function of the overlap $c_2 = \bm{\rho}_{{\rm SS},2}^\dagger \bm{\rho_0}$ for the setup in (a). Data represents 300 initial states $\bm{\rho}_0$ sampled at random. The steady state is computed using Eq. \ref{['eq:result2']} with $\epsilon = 0.0001$ (solid line), and Eq. \ref{['eq:result1']} (markers). (c) Same as in (a), but the qubits only undergo collective decay with jump operator $L_1$. (d) Concurrence in the steady state as a function of the overlap $\tilde{c}_2 = \bm{\mathcal{U}}_2^{\dagger}\bm{\rho}_0$ for the setup in (c). Data represents 300 initial states $\bm{\rho}_0$ sampled at random. The steady state is computed using Eq. \ref{['eq:result2']} with $\epsilon = 0.0001$ (solid line), and Eq. \ref{['eq:Albert']} (markers).
• Figure 3: Two ensembles of qubits under collective decay:(a) The system is initially prepared with all qubits in $A$ pointing up and all qubits in $B$ pointing down. This state can be written in the total angular momentum basis as $\ket{\psi_{\rm dif}}=\sum_{S=0}^{N/2}\sqrt{p(S)}\ket{S,0}$, where the distribution $p(S)$ is defined by the Clebsch-Gordan coefficients $p(S) = \vert \langle S,0 \vert N/4, +N/4, N/4, -N/4\rangle \vert^2$. We consider the two ensembles to be inside an optical cavity such that they can undergo collective decay characterized by the jump operator $L_1 = S_- = S_-^A + S_-^B$ and rate $\gamma$. (b) The collective decay jump operator $L_1 = S_-$ causes the population on each sector of $S$ to decay into the state $\vert S,-S\rangle$ with the final population in each of these states given by $p(S)$ as indicated by the size of the circles. The steady state corresponds to a mixture of all such states $\rho_{{\rm SS},L_1} = \sum_{S=0}^{N/2}p(S) \vert S,-S\rangle \langle S \vert$, since the initial state does not have any coherences of the form $\ket{S,-S} \bra{S',-S'}$ for $S\neq S'$ (see main text). (c) The quantum Fisher information of each state $\ket{S,-S}$ is plotted as a function of $S$ showing a rapid decay as $S$ increases. The inset shows the distribution $p(S)$ as a function of $S$, which peaks close to $\sqrt{N}$ as discussed in the main text. Here $N=60$.
• Figure 4: Enhancing the QFI using balanced decay processes:(a) The action of balanced collective jump operators $L_1 = S_-$ and $L_2 = S_+$ causes an initial Dicke state $\ket{S,M}$ to evolve into state $\rho_{{\rm B},S}$ which is an equal mixture of all Dicke states $\ket{S,M}$ with fixed $S$ but different $M$. (b) The quantum Fisher information of each state $\rho_{{\rm B},S}$ is plotted as a function of $S$. The numerical values (markers) computed with Eq. \ref{['eq:QFI2']} agree with the analytical expression found in \ref{['app_D']}. Here $N=60$. (c) Schematics of a single run of the protocol described in the main text. First, the system is initialized in the state $\ket{\psi_{\rm dif}}$. Subsequently, it decays collectively with the jump operator $L_1$ until it reaches the steady state $\ket{S,-S}$ with probability $p(S)$. Finally, we turn on the $L_2$ process and let the system evolve under the action of the balanced jump operators $L_1$ and $L_2$, reaching the steady state $\rho_{{\rm B},S}$. Additionally, if we wanted to do some post-processing to only keep evolutions that lead to low values of $S$ (high QFI), we could include measurements during the first evolution (dotted line). In a cavity setup, this could be done by continuously collecting the photons leaking out of the cavity while the system evolves under $L_1$. By doing so, we could know the specific value of $S$ and then discard the experimental runs corresponding to a large $S$ value. (d) Quantum Fisher information of different steady states as a function of the system size $N$. The markers represent steady states obtained numerically with Eq. \ref{['eq:result2']} using $\epsilon = 0.00001$. The solid line is the analytic expression for the QFI obtained for the protocol after many runs, while the dashed curve represents the QFI for the state $\rho_{{\rm SS},L_1}$ obtained for evolution under $L_1$ only [Eq. \ref{['ssUnb']}]. Details on how these two curves are obtained are provided in \ref{['app_D']}. Purple triangles correspond to the steady state obtained when the initial state $\ket{\psi_{\rm dif}}$ is evolved under both $L_1$ and $L_2$.
• Figure 5: Effects of population imbalance:(a) The distribution $p(S)$ corresponding to the initial state $\ket{\psi_{\rm dif}}$ is shown for different imbalance values $\eta$. (b) The quantum Fisher information of each state $\ket{S,-S}$ is plotted as a function of $S$. Different markers correspond to different values of the imbalance. (c) The quantum Fisher information of each state $\rho_{{\rm B},S}$ is plotted as a function of $S$. Different markers correspond to different values of the imbalance. (d) Quantum Fisher information of different steady states as a function of the imbalance $\eta$. All the data was obtained with Eq. \ref{['eq:result2']} using $\epsilon = 0.00001$. The circle markers signal the QFI of the state $\rho_{{\rm SS},L_1}$ obtained after the evolution under $L_1=S_-$, while the diamonds represent the QFI of the state at the end of the protocol. The solid purple line corresponds to the ratio of these two quantities. (e) Dynamics of the observables $\langle S_z^A \rangle$ and $\langle S_z^B \rangle$ when the system evolves under the jump operator $L_1$ and rate $\gamma$ into the steady state $\rho_{{\rm SS},L_1}$. Orange curves correspond to $\langle S_z^A \rangle$, while the pink ones correspond to $\langle S_z^B \rangle$. The solid lines represent the dynamics when the system is initialized in $\ket{\psi_{\rm dif}}$, while the dashed lines correspond to the initial state $\ket{\psi_{{\rm dif}_2}}$. For all panels, we consider $N=50$.