Information preservation of two qubits in a structured environment

The environment-induced decoherence of a quantum open system makes it fundamentally import to preserve the initial quantum information of the system in its steady state. Here we study information preservation of two maximally entangled qubits lying inside a photonic-crystal waveguide with semi-infinite cavity-array structure. We generalize our study to arbitrary position and arbitrary frequency detuning of the qubits. We find that for weak qubits-waveguide couplings, the information preservation greatly depends on the position and the frequency detuning of the qubits, while for strong couplings, both of these dependence is significantly weakened. Interestingly, by suitably choosing the position and the frequency of the qubits, high information preservation could be achieved for both weak and strong couplings, irrespective to Markovian or non-Markovian dynamics. Physically, we analytically verify that the ability of information preservation is indeed determined by the existence of the bound states of the entire system, but the probability of information preservation is closely related to the probability of the initial state of the qubits in the bound states. Our results provide an alternative route getting high information preservation without any external controls of the system.


Introduction
The environment-induced decoherence of a quantum open system is not only a fundamental topic in the quantum physics, but also an essential obstacle in quantum information and computation [1]. Thus the study of the decoherence dynamics attracts large amounts of attentions, including the Markovian dynamics [2][3][4][5] and the non-Markovian dynamics [6][7][8][9][10][11]. Traditionally, the Markovian memory less environment may result in the completely loss of the initial quantum information of the system in the long-time limit [12][13][14][15][16][17]. This happens for weak system-environment couplings. For strong couplings, the non-Markovian environment memory will feedback some information to the system, which can partially preserve the initial quantum information of the system [11,[15][16][17][18][19]. In return, this information preservation greatly inspires people to investigate the non-Markovian dynamics in many quantum systems in physics [20][21][22][23][24][25][26][27][28][29][30], biology and chemistry [31][32][33][34][35][36][37]. Indeed, the non-Markovian effect is only one essential ingredient for information preservation, another key ingredient is the existence of the bound states (it is the localized modes in the continuous-variable systems) in the entire system consisting of both the system and the environment [19,38,39]. As a result, the formation of bound states (localized modes) plays a crucial role in quantum information preservation [40][41][42][43][44][45][46][47][48][49][50]. However, all these above studies were mostly concentrated on the Lorentz-type or Ohmic-type spectral density, where it was shown that strong non-Markovian effect could induce high information preservation.
In quantum optics, the qubits (atoms or impurities) coupled to a field within a photonic crystal exhibit strong non-Markovian dynamics [11,16,19,29,39,[51][52][53][54][55][56][57][58][59][60][61][62][63][64]. This is because the band-gap or finite-band structure of the spectral density brings bound states in the band gaps and the dissipationless non-Markovian dynamics of the system in the long-time scale [11,16,19,39,[56][57][58][59]. One typical material is the photonic-crystal waveguide with semi-infinite cavity-array structure [11,16,19,54,[56][57][58][59]61], which is of experimental advances since both the qubit-waveguide coupling strength and the cavity-cavity interaction in the waveguide are easily controllable by adjusting the geometry structure [65]. It provides a spectral density of finite-band type, which may lead to the oscillating dissipationless dynamics of the qubits due to the competitions of the bound states. This oscillating dynamics is not observed for other spectral densities, e.g. the lorentz type [11,19,22], the Ohmic type [11,16,18,19], since there is no more than one bound state. Anyhow, photonic-crystal waveguide with semi-infinite cavity-array structure provide a rich platform for the investigation of decoherence dynamics and information preservation. However, the previous studies [11,16,19,54,[56][57][58][59]61] were mostly focused on the cases where the system frequency is fixed (mostly it is resonant with the optical field of the cavity) and the system position is fixed (generally it is only placed inside the edge cavity of the waveguide). In these cases, it was also showed that strong non-Markovian effect brings high information preservation.
In this work, we study two initially maximally entangled qubits coupled to a photonic-crystal waveguide with semi-infinite cavity-array structure. Here, we will extend the previous studies to arbitrary frequency detuning and arbitrary position of the qubits to further explore information preservation. We concentrate on the maximum probability of the initial state of the system preserved in the steady state. Interestingly, due to the translational symmetry breaking of the semi-infinite cavity-array structure, different positions of the qubits would bring totally different shapes of the spectral density, which will induce more zero points of the spectral density inside the band similar to those in the band gaps. In terms of numerical analyses and analytical explanations, we find that the information preservation of the qubits is both of strong position-dependence and detuning-dependence for weak couplings, while for strong couplings, these dependence is remarkably suppressed. Amusingly, except for high information preservation induced by the strong couplings, when the frequency detuning of the qubits locates at the zero points of the spectral density, the initial information of the qubits could also be highly preserved for weak couplings. The former corresponds to strong non-Markovian effect, the latter may be Markovian or weak non-Markovian effect. Physically, we verify that the ability of the information preservation is indeed decided by the existence of bound states, but the probability of information preservation is closely related to the probability of the initial state in the bound states.
The paper is organized as follows. In section 2, we introduce the Hamiltonian of the system. In section 3, the exact dynamical equation of the system is derived. Then in section 4, we give the numerical analyses and analytical explanations of the exact results. Finally, conclusions and discussions about the applications of our model and the relevant experimental implementations are given in section 5.

Hamiltonian
We consider an open quantum system of two qubits (e.g. two-level atoms) embedded in the qth cavity of the photonic-crystal waveguide. The photonic-crystal waveguide consists of semi-infinite cavities, which is regraded as the environment of the qubits. A schematic diagram is shown in figure 1. The Hamiltonian of the entire system reads (set ℏ = 1) [54,56] Here the first term describes the free Hamiltonian of the qubits. For the ith qubit, σ + q i = 1⟩ q i ⟨0 and σ − q i = 0⟩ q i ⟨1 are the transition operators between the ground state |0⟩ q i and the excited state |1⟩ q i with the transition frequency ω S . The index q denotes that the qubits are located in the q-th cavity of the cavity array, as shown in the last term. The second term and the third term represent the Hamiltonian of the photonic-crystal waveguide with semi-infinite cavities in one dimension where j labels the position of the cavities in the coordinate space. We assume each of the cavities has a single-mode field with the same frequency ω 0 . a † j and a j are the creation and annihilation operators of the jth cavity mode. The interaction between the cavities is of tight-binding type, with the neighboring interaction strength ξ 0 . Finally, the last term gives the interaction of the qubits with the cavity where ξ q is the coupling strength between them. For convenience, we would like to transform the Hamiltonian (1) from the coordinate space to the momentum space. By making the sine-type Fourier transformation that a j = 2 π +∞ k=1 sin(kj)b k , the Hamiltonian becomes a standard spin-boson model [3,5] where k represents the momentum of the cavities in the momentum space. In this representation, the dispersion relation of the cavity array is ω k = ω 0 − 2ξ 0 cos k. Moreover, the qubits interact with the cavity of the momentum k taking the coupling strength Here we see that due to the translational symmetry breaking of the semi-infinite cavities, the effective interaction strength of this structured environment depends strongly on the position q of the qubits. Thus the dynamics of the qubits will closely depend on the position, which did not appear in the previous works [16,19,[57][58][59]61].

Exact dynamics of the system
We assume initially the qubits are maximally entangled and the waveguide is vacuum, i.e. |ψ SE (0)⟩ = |ψ S (0)⟩ ⊗ |ψ E (0)⟩ with |ψ S (0)⟩ = |ϕ⟩ ≡ 1 √ 2 (|1, 0⟩ + |0, 1⟩) and |ψ E (0)⟩ = |{0}⟩ E . Due to the total particle number of the qubits and its environmentN = is conserved and the two qubits are identical, one can write the time evolution state of the entire system as Here the time-dependent coefficient c q (t) represents the probability amplitude of the entire system in its initial entangled state, and the coefficient c k (t) is the probability amplitude when the single excitation emits to the cavity with the momentum k.
In addition, by tracing over the degrees of the environment, the evolution state of the qubits could be obtained as where P(t) ≡ c q (t) 2 is the probability of the qubits maintained in their initial entangled state. This is closely related to the fidelity [66,67] F (t) = ⟨ψ S (0) |ρ S (t) |ψ S (0)⟩ = c q (t) = P(t). So the probability P(t) could be safely used to quantitatively reflect the information preservation of the system. By the way, take a derivative of the density matrix ρ S (t), one can obtain the following differential equation (see appendix A) q Jq(∆ω) zero points The Lindblad form of equation (6) ensures that the coefficient γ(t) is just the dynamical decay rate of the qubits. Thus if γ(t) ⩾ 0 ( dP(t) dt ⩽ 0 ), the qubits undergo a Markovian dynamics, else γ(t) < 0 ( dP(t) dt > 0) indicates a non-Markovian dynamics [6][7][8][9][10]. Qualitatively speaking, the larger positive dP(t) dt is, the stronger non-Markovian effect is. It should be pointed out that γ(t) = 0 ( dP(t) dt = 0) means that the system undergoes a dissipationless Markovian dynamics. In this work, we do not pay much attention on the quantitative property of the non-Markovianity of the system dynamics, we mainly focus on how much initial information could be preserved after the system reaches its steady state.
Since the entire system of the qubits and the waveguide is a close one, its dynamics is determined by the Schrödinger equation. Thus the equations of motion of c q (t) is derived as Here we definec q (t) ≡ c q (t)e iωSt to clearly see the detuning dependence of the dynamics of the qubits. The memory kernel of the waveguide takes the form where we set ∆ω ≡ ω−ω0 2ξ0 . And ∆ω S ≡ ωS−ω0 2ξ0 represents the frequency detuning between the qubits and the cavity field in the waveguide. Here J q (ω) ≡ 2π +∞ k=1 |V qk | 2 δ(ω − ω k ) is the spectral density of the system-environment interaction. From the dispersion relation and the coupling strength in the Hamiltonian (2), we get the spectral density as with η = ξq ξ0 a dimensionless coupling strength between the qubits and the waveguide. Obviously, if q = 1, it reduces to the previous result J 1 (∆ω) = η 2 2ξ 0 1 − (∆ω) 2 [54,56]. For q = 2, 3, 4, the expressions of J q (∆ω) are shown in table 1. For arbitrary q, the expression of J q (∆ω) could be obtained from the general expansion of sin(qθ) ( set θ = arccos(∆ω)) Obviously, the key equation (8) is an integrodifferential equation. The integration in the kernelf (t − τ ) reflects the non-Markovian memory back action effect of the waveguide on the qubits. Clearly, from equations (9) and (10), one find that except for the coupling strength η, the dynamics of the qubits also depends on the position q and the frequency detuning ∆ω S of the qubits.
From the spectral density (10), we see that the increasing of η only enhances the distribution weight for given ∆ω. It does not change the shape of the spectral density J(∆ω). However, varying the position of the qubits will greatly change the shape of the spectral density, and the corresponding weight as well, as clearly shown in figure 2, where outside the band, the spectral density will always take the value of zero. When the qubits are put in the position of q, the spectral density inside the band appears q peaks and q + 1 zero points. Interestingly, when q is odd, the center (∆ω = 0) is a peak, while when q is even, it becomes zero (i.e. J(∆ω) = 0). In table 1, we list the zero points for q ⩽ 4 inside the band of the spectral density. Below, we will show that when the frequency detuning of the qubits is set at the zero points of the spectral density, no matter it is inside or outside the band, the dynamics of the qubits will be qualitatively different from other nonzero points.

Exact results
Next, we will analyze the exact dynamics of the qubits by numerically solving the equation of motion (8) and analytically explaining the results.

Numerical results
After numerous calculations, we find that there are two types of dynamical behaviors of P(t) changing with the increasing of η. One type is the case that ∆ω S locates at the nonzero points of J(∆ω), e.g. ∆ω S = 0 as shown in figures 3 and 4. The other type corresponds to the case that ∆ω S situates at the zero points, including both inside and outside the band, of J(∆ω), e.g. ∆ω S = 1 as shown in figures 5 and 6.
On one hand, in the case of nonzero points of J(∆ω), with the increasing of η, when η surpasses a critical value η c , the dynamics of P(t) changes from the exponential (or exponential-like) decay Markovian (weak non-Markovian) behaviors to the oscillating (strong non-Markovian) behaviors. To quantitatively see how much information could be finally preserved in the qubits, we extract the information of the maximum value of P(t) that could be achieved after the qubits reach their steady state (or the long-time limit), denoted as P m s . That is, Figure 4(b) shows the maximum steady probability P m s versus the coupling strength η. With the increasing of η, P m s increases from zero (when η ⩽ η c ) to nonzero values (when η > η c ), even to a very high value (≃ 0.9) for large enough η. That is, strong couplings between the qubits and the waveguide induces the strong information exchange between them. This is a strong non-Markovian effect [11,16,19,54,[56][57][58][59] and the initial quantum information of the qubits could be highly preserved. This is a commonly phenomenon for other spectral densities as well, e.g. Ohmic-type spectral density [18].
On the other hand, for the zero points of J(∆ω), the dynamics of P(t) shows totally different behaviors for weak η compared to the case of the nonzero points. As shown in figures 5 and 6(a), for a weak coupling (e.g. η = 0.1), the dynamics of P(t) undergoes a dissipationless steady dynamics after an initial small decay. Thus the initial information of the qubits is highly preserved. With a small increase of η, e.g. η = 0.7, the steady value of P(t) decreases. Figure 6(b) displays the maximum steady probability P m s versus η. P m s first decays from a high value for small η, then surpassing a critical η c , it arises again to a high value again. This is totally different from the case with nonzero points of J(∆ω) as shown in figure 4(b).

Analytically explanations
Now let us roughly give a physically explanation about the reason why the dynamics is so different for different frequency detunings. From the dynamics of P(t) shown in figures 5 and 6(a), we see that dP(t) dt is almost smaller or equal to zero for η = 0.1. Thus this may be a Markovian or Markovian-like behavior. It is known that [11,56] for this model in the Born-Markovian limit, the constant decay rate γ BM ≃ J(∆ω S )/2. Hence, when ∆ω S is set at the points of J(∆ω) = 0, the decay rate γ BM ≃ 0 leading to a dissipationless Markovian dynamics. It seems that the system evolves almost closely. For a little lager η (e.g. η = 0.7), the Born-Markovian approximation is broken. So the dynamics of the system becomes a non-Markovian one, and the initial information of the quantum system may be partially lost. This is why the maximum steady probability P m s decreases first in figure 6(b). Moreover, with the increasing of η, the information exchanges strongly between the system and its environment, and the strong non-Markovian effect again takes advantage for information preservation of the system. Fortunately, we can give an analytical solution of the above two cases, and make a more accurate explanation in terms of the bound states. As shown in appendix B, the key integrodifferential equation (8) takes an analytical solution with the dissipationless steady value [19] Here {ω B i } are the poles of c q (t) at the complex number space after the Laplace transformation. They are induced by the discontinuity of the spectral density. The coefficients {Z i } give the corresponding amplitudes of the poles. Physically, these poles are just the eigen-energies, obeying the equation of (B7), of the entire system as shown in the appendix C. They correspond to the eigenstates (bound states) of the entire system in the single-excitation subspace |ψ i There we also prove that the amplitudes of the poles Z i = |d i ϕ | 2 . Thus the maximum steady probability is finally derived as, i.e.
Thus the ability of information preservation is indeed decided by the existence of the bound states, but the probability of information preservation is closely determined by the probability of the initial state of the system in the bound states of the entire system. For q = 1, there are two bound states lying outside the band edges. At the resonant point ∆ω S = 0, the two corresponding eigen-energies are ω B± = ω 0 ± 2η 2 ξ0 √ 2η 2 −1 existing when η > η c± = 1 [54]. They have two equal coefficients Z ± = 1 2 (1 − 1 2η 2 −1 ). So the steady amplitude is with Θ(x) the Heaviside step function, it equals 1 when x > 0, otherwise it is 0. Obviously, when η ⩽ 1, Z + = 0, the steady value of P(t) = 0 indicating that all the initial information of the qubits is lost. When η > 1, the steady amplitude a 1 (t) is an oscillating time function. Thus the maximum steady probability P m s = 4|Z + | 2 is obtained at the times 2η 2 ξ0 √ 2η 2 −1 t = 2kπ, k = 0, 1, 2, . . .. It is totally determined by the coefficient Z + who increases monotonously with the increasing of η. So the larger η is, the larger P m s is. In the strong coupling limit that η ≫ 1, P m s → 1. As shown in figure 4(b), this analytical result (the red line) coincides well with the numerical one (the blue circle). Thus, traditionally, for the nonzero points of J(∆ω), one could highly preserve the initial information by enhancing the coupling strength η [11,16,19,39,[56][57][58][59].
Interestingly, for the zero point with detuning ∆ω S = 1, the two bound states taking the energies 2η − 1)ξ 0 exist when η > η c+ = √ 2 and η > η c− = 0, respectively. The corresponding coefficients Z ± = 1 2 (1 ∓ 1 √ 2η∓1 ). Then the steady amplitude becomes When 0 < η ⩽ √ 2, only the bound state with energy ω B− exists. Thus the steady P(t) = |Z − | 2 is constant with Z − in inverse proportion to η. The smaller η is, the larger P(t) is. When η → 0, Z − → −1, and P(t) → 1. This is totally different from the case with ∆ω S = 0. When η > √ 2, the steady probability oscillates periodically due to the competition of the two bound states. Thus we can get the maximum steady probability P m s = (Z + + Z − ) 2 at 2ξ 0 ( 1 η−1 + η − 1)t = 2kπ, k = 0, 1, 2, . . .. The larger η is, the larger P m s is. In the limit that η ≫ 1, P m s → 1. From figure 6(b), we see that the analytical and numerical results coincide well with each other. Thus, when the detuning frequency locates at the zero points of J(∆ω), even in the weak coupling case, there exists one bound state leading to the dissipationless steady dynamics of the system. As a result, both weak η and strong η could bring high information preservation.
Furthermore, figure 7 clearly displays the maximum steady value P m s verse different ∆ω S in the cases of η = 0.1 and η = 3.0. We see that in the weak coupling case, only the zero points of J(∆ω), both inside and outside the band, takes nonzero probability (indeed almost to 1) of the information preservation. In this case, the Markovian effect with single bound state contributes more to high information preservation. On the contrary, for a large enough coupling strength, the ability of information preservation is almost irrespective of the detuning ∆ω S . In this case, the strong non-Markovian effect takes more advantage. That is, suitably setting the detuing ∆ω S , high information preservation could be obtained from both weak couplings and strong couplings.

Analyses for q = 2
Now let us further check the case for q = 2. We first focus on the special zero point of the spectral density ∆ω S = 0, which is just the resonant point. As shown in figure 8(a), for a small η, the qubits could highly preserve its initial entanglement. For a large enough η, the steady pattern of P(t) is a periodic oscillation with two different oscillating peaks. As shown in figure 8(b), with the increasing of η, the maximum steady P(t) decreases first, then exceeding a critical η c , it increases again. Below let us analytically check the results from the view of bound states.
At the zero point ∆ω S = 0, there are totally three bound states. Two of them lies outside the band edges, 2η)ξ 0 with the condition of η > η c± = 1/ √ 2. The additional bound state locates at the center of the band ω B0 = ω 0 with the condition η > η c0 = 0. The corresponding coefficients are Z ± = 1 2 2η 2 −1 2η 2 +1 and Z 0 = 1 2η 2 +1 . Thus the steady value is Obviously, when 0 < η ⩽ 1/ √ 2, the steady value a 2 (t) is only determined by the bound state with energy ω B0 , whose amplitude Z 0 decreases monotonously with the increasing of η. The smaller η is, the larger P(t) = |a 2 (t)| 2 is. When η → 0, P(t) → 1. When η > 1/ √ 2, all the three bound states ω B± and ω B0 contribute to the steady a 2 (t). The competition among them leads to two oscillating peaks in the steady dynamics for η = 1.5 and η = 3.0. In the limit that η ≫ 1, the three coefficients approximate to Z 0 → 0 and Z ± → 1 2 . Thus P(t) → | cos( √ 2ηξ 0 t)| 2 with P m s → 1 as well. This again proves that the existing of bound states induces the information preservation of the system, but high information preservation may happen both for very weak and very strong couplings, no matter the non-Markovianity of the dynamics.   For q = 2, we also display the maximum steady value P m s for different ∆ω S at η = 0.1 and η = 3.0, as shown in figure 9. Similar to the case of q = 1, for weak η, the system could highly preserve its initial entanglement at the zero points of J(∆ω), while for strong η, the information preservation is almost independent of ∆ω S .

Analyses for arbitrary q
Finally, in figure 10, we present the maximum steady probability P m s versus the position q for different coupling strength η at the resonant point ∆ω S = 0. As expected, for small η, one can get almost zero P m s at odd positions, and high P m s at even positions. This is because the resonant point ∆ω S = 0 is a zero point of J(∆ω) at even positions. Moreover, for even q, P m s decreases as q increases. This indicates that for weak η, q = 2 maybe the best position for information preservation. However, with the increasing of η, the values of P m s for the even and odd q are slowly converging on each other, and finally they are almost identical. Hence the weaker the coupling strength is, the stronger the dependence of the information preservation on the position is.

Conclusions
Now let us conclude our results. We study the information preservation of two qubits coupling to a photonic-crystal waveguide with semi-infinite cavity-array structure. Initially the qubits are maximally entangled. The spectral density is of finite-band type. For different positions of the qubits in the cavity array, the shape of the spectral density is greatly changed. This may induce more than three bound states of the entire system. From the exact dynamical equation, we clearly see that the dynamics of the qubits strongly depends both on their position and the frequency detuning from the cavity field. In terms of numerical analyses and analytical explanations, we study the probability of the qubits preserved in their initial entangled state and analyze its dependence on the position and the detuning.
We find that for weak couplings, the initial entangled information will be lost completely after the qubits reaches its steady state when the frequency detuning of the qubits locates at the nonzero points of the spectral density. While when the frequency detuning locates at the zero points of the spectral density, the initial entanglement could be highly preserved due to the appearance of a bound state of the entire system. The weak couplings correspond to a Markovian or weak non-Markovian dynamics of the qubits. For strong couplings, the bound states always induce high information preservation no matter the frequency detuning is at the zero or nonzero points of the spectral density. This is always corresponding to a strong non-Markovian dynamics. We also find that, for the weak couplings, the information preservation strongly depends on the position of the qubits. At resonant case, the even and odd positions show totally different behaviors due to they correspond to the zero and nonzero points of the spectral density, respectively. While for strong couplings, this dependence almost disappears.
Therefore, for weak couplings, the information preservation of the qubits depend strongly on both the position and the frequency detuning of the qubits. While for strong couplings, the information preservation is not only detuning-independent, but also position-independent. Physically, our analytical results show that the ability of information preservation is indeed determined by the existence of the bound states. However, the probability of information preservation is essential determined by the probability of the initial state of the qubits in the bound states.
In addition, the probability of information preservation in this work is the probability of the qubits in their initial entangled state. Hence our results could easily be verified in current experiments. Moreover, since weak couplings are usually more easily realized in experiments, our results provide a new and easy way for achieving high information preservation.

Discussions
Below we discuss the applications of the model under our consideration and the relevant experimental implementations.
In this work, we consider a system of two qubits coupled to a semi-infinite structure photonic-crystal waveguide. In fact, the semi-infinite structure of the waveguide could also be realized in other experimental platforms including optical fibers [68,69], superconducting microwave transmission lines [70,71], photonic nanowires [72] and plasmonic waveguide [73]. It has been widely utilized for the investigation of the geometric constraint effect of the electromagnetic field on the cavities or the qubits [74]. It hosts many interesting physical effects due to one end of the waveguide is terminated, such as the dynamical Casimir effect [75], giant Lamb shift [71,76,77], and suppression of spontaneous emission [71,[76][77][78][79][80]. Thus it attracts lots of applications in dynamical control of single photons [81], bright single-photon generation [82], single-photon detection [83], effective cavities with atomic mirrors [84], quantum information processing [85,86], and quantum entanglement generation [45].
It should be pointed out that the position dependence of the qubits in the photonic-crystal waveguide is testable experimentally. The position q of the qubits is equivalent to the distance x q = qR with R the length of the unit cell of each cavity. One way to measure the distance x q is directly from the scanning electron micrograph of the waveguide. For example, for photonic-crystal waveguide with cavity array structure [65], the length R = 1.5 µm for a photonic lattice constant of 245 nm. The cavity array has the total number of N = 50, which gives a total physical length of the waveguide L = 75 µm. Generally, the field of vision of a scanning electron microscope is about 100 µm [87]. Thus the whole waveguide could be scanned experimentally. In this way, the distance x q could be measured directly from the scanning electron micrograph.
Another way to detect the distance x q is in terms of the reflection spectrum of the input light. This is according to the physics that the qubits behave like a mirror, so an effective cavity is formed between the terminated end of the waveguide and the qubits [81,84]. The length of this effective cavity just equals to the distance x q . Then one can test the length of the effective cavity from the reflection spectrum of the input light using the optical spectrum analyzer. This is a standard approach to measure the cavity length [88]. That is, the valleys of the reflectivity is closely related to the phase φ of the round trip of the input light propagating inside the effective cavity. From the value of φ, one can easily get the distance x q from the relation of φ = 4π λ nx q where λ is wavelength of the input light and n is the refractive index of the waveguide. It is worth noting that, generally, the position of the qubits is not easy to be varied. However, in the platform of superconducting microwave transmission lines, which is terminated by a superconducting quantum interference device (SQUID), it was shown [89,90] that the distance x q could be changed effectively by varying the boundary condition at the terminal end of the transmission line through the magnetic flux threading the SQUID ring.

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

Appendix A. Derivation of equation (6)
Below we will derive the Lindblad form of the differential equation (6) of the density matrix ρ S (t) (5). First, take a derivative of equation (5) as Generally, it is not directly to get the Lindblad form of equation (6) from (A1). But we can conversely to prove that equation (6) is equivalent to (A1). By rewriting the spin operators as where |ϕ ′ ⟩ = 1 √ 2 (|1, 0⟩ − |0, 1⟩) S , one can easily rewrite the right hand side of equation (6) as That is, the Lindblad terms for the two qubits are the same with each other. Thus equation (6) equals to Compared the equation (A4) with equation (A1), one find that only if the rate in the Lindblad from satisfies the equations (A1) and (A4) are equivalent to each other. Hence one can say the derivative of (5) is just the differential equation (6), which is of the Lindblad from with the coefficient γ(t) taking the expression of (A5).

Appendix B. Analytical solution of the exact equation (8)
Below we will give the analytical solution of the exact equation (8) according to [19]. For clarity, let us start from the original equation of motion of c q (t) d dt c q (t) + iω S c q (t) +ˆt where the memory kernel of the environment takes the form Using the Laplace transformation C q (z) =´∞ 0 dt c q (t)e izt , one can obtain with Σ(z) the environment induced self-energy Then the analytical solution of c q (t) could be derived from the inverse Laplace transformation of C q (z) as Here the first term is the contribution from the frequencies {ω B i }, which are the poles of C q (z) induced by the discontinuity of the spectral density. This gives rise to dissipationless dynamics. The second term is a dissipation dynamics coming from the non-analyticity of the spectral density, which usually produces non-exponential decays characterized by the dissipation spectrum D(ω) [19]. Thus, in the long-time limit, the second term finally disappears, and only the first term retains. That is, Mathematically, the frequencies {ω B i } are determined by the equation of While the coefficient is the corresponding residue.
Therefore, the long-time solution (B6) tell us that, the maximum steady probability of the system in its initial state is generally decided the residues {Z i }.
We can easily get the corresponding eigen-energy obeying the equation with d ϕ is the probability amplitude of the initial state |ψ SE (0)⟩ = |ϕ⟩|{0}⟩ E in the bound state of the entire system.
Obviously, the equation of eigen-energy (C4) is exactly the same as the equation of the poles ω B i (B7). Thus the frequencies {ω B i } are just the eigen-energies of the entire system corresponding to the bound state (C1). If there are more than one solutions of E 1 , different E i 1 leads to different d i ϕ and different |ψ i 1 ⟩, then there will be more than one bound states. Thus the existence of ω B i is just equivalent to the formation of the bound state (C1) in the entire system.
Moreover, for the initial state |ψ SE (0)⟩, the amplitude Z i in (B8) is exactly the square of the absolute of the amplitude d ϕ in (C5), i.e. Z i = |d i ϕ | 2 . This is similar to the single qubit case [93]. In this way, the maximum steady probability (B9) becomes It is determined by the probability of the initial state of the system in the bound states of the entire system.