Accelerating relaxation through Liouvillian exceptional point

We investigate speeding up of relaxation of Markovian open quantum systems with the Liouvillian exceptional point (LEP), where the slowest decay mode degenerate with a faster decay mode. The degeneracy significantly increases the gap of the Liouvillian operator, which determines the timescale of such systems in converging to stationarity, and hence accelerates the relaxation process. We explore an experimentally relevant three level atomic system, whose eigenmatrices and eigenspectra are obtained completely analytically. This allows us to gain insights in the LEP and examine respective dynamics with details. We illustrate that the gap can be further widened through Floquet engineering, which further accelerates the relaxation process. Finally, we extend this approach to analyze laser cooling of trapped ions, where vibrations (phonons) couple to the electronic states. An optimal cooling condition is obtained analytically, which agrees with both existing experiments and numerical simulations. Our study provides analytical insights in understanding LEP, as well as in controlling and optimizing dissipative dynamics of atoms and trapped ions.


I. INTRODUCTION
Open quantum systems coupled to environments will relax toward a stationary state.The relaxation processes have rich properties from both dynamic and thermodynamic perspectives.Often an important question is to control the relaxation time [1][2][3][4], for instance, on a timescale as short as possible [see Fig. 1(a)].This problem is of great relevance to cases where one is concerned with properties of stationary states, such as ground state laser cooling [5][6][7][8][9], or aims to generate quantum states for quantum applications [10][11][12][13].
Starting from an arbitrary initial state, the relaxation timescale is largely characterized by the slowest decay mode of the Liouvillian operator.The gap is defined as modulus of the real part of its eigenvalue λ 1 [14][15][16], as depicted in Fig. 1(b).Therefore, relaxation speeding is achieved through increasing the gap.An alternative approach to speed the relaxation is offered by the so-called Mpemba effect [17][18][19][20][21][22], where an unitary operation on the initial pure state removes its overlap with the slowest decaying mode [1,2].This transformation can be exactly constructed provided that the initial state is a pure state.
In this paper, we show that, for an arbitrary initial state, if the slowest decay mode and its corresponding eigenvalue coalesce with a faster decay mode, one can maximize the gap and thus accelerates dynamics to reach stationary states.Our study exploits the novel nature of exceptional points (EPs), which are hallmarks of non-Hermitian systems [23][24][25][26][27][28][29][30].EPs are specific points in parameter space, where two or more eigenvalues of a non-Hermitian operator and their corresponding eigenvectors coalesce [31].The origin of non-Hermiticity is the coupling between the system and the environment.Liouvillian superoperator, which captures the time evolution of an open quantum system, is non-Hermitian.Therefore, it can exhibit EPs (referred to as Liouvillian EPs, or LEPs) [26,31].Properties of LEPs, with diverse unusual effects, have attracted considerable attention currently [32][33][34], such as dissipative phase transition [35][36][37], non-Hermitian skin effect [38], signatures of LEPs in the dynamics [39][40][41] and so on.
Here, we focus on analytical understanding and applications of LEPs in quantum control [42,43], and further explore how to accelerate the relaxation towards to stationarity in Markovian open quantum systems via LEPs.The basic mechanism underpinning our study is based on the fact that, at LEPs both the slowest decay mode and corresponding eigenvalue coalesce with a faster decay mode and the corresponding eigenvalue.We show that, when the stationary state of the system is unique and independent on system parameters, one can set the parameters at the LEP to speed up the relaxation process significantly.For certain quantum dynamic processes, such as ground state cooling, how to convergent to stationarity as quickly as possible is often a concern in the actual process of quantum manipulation.In addition, we find that relaxation processes can be further accelerated by periodically modulating the dissipation strength, i.e., the Floquet modulation can overcome the gap limit of the static case and realize faster relaxation [see Fig. 1(c For an open quantum system with an arbitrary initial state, the timescale of approaching to the final stationary state ρss is related to the slowest decay mode (with eigenvalue λ1) of Liouvillian superoperator.By tuning the parameter of the system near the LEPs, where both of the slowest decaying mode and the corresponding eigenvalue are merged with a faster decaying mode and the corresponding eigenvalue (for instance λ2), the system dynamics approaches the stationary state in a much faster way.(b) This feature is evident from the Liouvillian spectrum.The stationary state ρss is characterized by the largest eigenvalue λ0 = 0.The other eigenvalues, characterizing the decay modes, have non-positive real part and always appear as complex conjugates.The Liouvillian spectral gap (g = −Re[λ1]) determines the relaxation timescale, and can reach its maximum value at LEP. (c-d) The gap at LEP can be further increased by the Floquet method (red line).In contrast with the static case (blue line), the gap under time-periodic modulation can be significantly increased, which means that relaxation process will be accelerated by applying the Floquet method.

d)]
. We apply our approach to analyze ground state laser cooling based on sideband transitions and electromagnetically induced transparency (EIT).Optimal conditions are obtained analytically, which have been demonstrated in recent trapped ion experiments [7,9].Our study reveals the importance of LEPs in practical applications and provides insights in seeking optimal conditions in quantum control of open quantum systems.
This paper is organized as follows.In Sec.II, we introduce the master equation of the Markovian open quantum systems.A general framework that connects to its dynamics and eigenmatrices of the Liouvillian superoperator is provided.This provides an intuitive picture to understand the relaxation and the gap.In Sec.III, we study the dynamics of dissipative three-level system.Eigenmatrices and eigenvalues of the corresponding Liouvillian superoperator are obtained analytically.Based on the analytical calculation, we reveal that the relaxation towards to stationary state can be accelerated by exploiting static and Floquet-modulated LEPs.Next, in Sec.IV applications for ground state laser cooling are demonstrated.Two experimentally relevant scenarios, i.e. sideband cooling and EIT cooling, are examined.Optimal cooling conditions are obtained at corresponding LEPs.We conclude in Sec.V.

II. LIOUVILLIAN GAP, DYNAMICS AND LEP
We consider an open quantum system evolving under Markovian dynamics, governed by master equation ρ(t) = Lρ(t), where the generator L, normally called Li-ouvillian superoperator, has the form [44,45] Here, ρ(t) is the state of the system at time t, H is the system Hamiltonian, and J α are quantum jump operators which provide coupling of the system to the environment.Since the Liouvillian L acts linearly on ρ(t), one can obtain information about the relaxation in terms of its eigenmatrices R i and the corresponding complex eigenvalues λ i via the relation LR i = λ i R i .Note that, duo to the Hermiticity of L, if λ i is complex, λ * i must also be an eigenvalue of L [1,26,36,46].Therefore, the eigenvalues are symmetrically distributed with respect to the real axis as shown in Fig. 1(b).
The stationary state of the system under consideration is then given by the density matrix ρ ss such that Lρ ss = 0, i.e., ρ ss = R 0 , which corresponds to the zero eigenvalue λ 0 = 0 and is independent of the initial state.If the eigenvalues are ordered by decreasing their real parts, it is known that the negative real parts of the eigenvalues [47], Re[λ i>0 ] < 0, determine the relaxation rates of the system towards the nonequilibrium stationary state, and the corresponding eigenmatrices R i>0 are called decay modes [15,48].While the imaginary parts describe the oscillatory processes which may take place.We can then write the time dependence of the density operator from an initial state ρ in as where a i = Tr[L i ρ in ] are coefficients of the initial state decomposition into the eigenmatrices of Here R i and L i are referred as right and left eigenmatrices (eigenmodes), respectively, and can be normalized by Tr[L i R j ] = δ ij .The trace preservation of the dynamics implies that Tr[ρ(t)] = Tr[ρ ss ] = 1 = Tr[L 0 R 0 ], and thus L 0 is the identity (L 0 = I).It also implies that Tr[R i⩾1 ] = 0, which means other right eigenmatrices do not correspond to quantum states.A particular interesting case is when eigenvalue λ i is real, where the corresponding eigenmatrix can be diagonalized [36].We can rewrite it as superposition of eigenstates from the diagonalization [36] with where |ψ i n ⟩ are eigenvectors of R i with eigenvalues p i n .With this definition, R ± i are arranged to proper density matrices.If λ i is complex, one can define a pair of eigenmatrices R i +R † i and i(R i −R † i ), then their corresponding eigenvalues are real (i.e.real and imaginary parts of λ i ).This allows to diagonlize their new eigenmatrices.
A fundamental role in the system dynamics is played by λ 1 (R 1 ), which possesses the slowest decay rate on the condition a 1 ̸ = 0. Then the Liouvillian gap, defined by [14,49] is thus an important quantity determining the timescale of the final relaxation to the stationary state.If the slowest decay mode λ 1 (R 1 ) coalesce with a faster decay mode when we set the system parameters at the so-called LEP, where λ 1 (R 1 ) = λ 2 (R 2 ), the gap will have extreme values g LEP .Consequently for long times one has where In such a case, the state would relax at the fastest rate with timescale 1/g LEP for arbitrary initial state.As we will demonstrate latter, this is useful in some quantum applications, where long relaxation timescales become impractical or even harmful to the coherence.Therefore how to quickly approach the steady state becomes necessary in these applications .

III. ANALYTICAL LEP THEORY OF DISSIPATIVE THREE-LEVEL SYSTEM
A. The model Consider a simple dissipative three-level system of Fig.
and a jump operator J = √ γ|a⟩⟨c|.Qualitatively, there will be competition between the reversible coherent coupling between |b⟩ ↔ |c⟩ at frequency Ω and the population loss of |c⟩ at a rate γ.As shown in Fig. 2(a), if Ω ≫ γ, Rabi oscillations occur before |c⟩ eventually decays to |a⟩ and dynamics of the system exhibits damped oscillations, which corresponds to under-damped (UD) dynamics.At very long times, the probability in the state |b⟩ tends towards zero since the system ends up in level |a⟩.If, on the other hand, Ω ≪ γ, we expect an over-damped (OD) evolution for the probability of state |b⟩, which tends exponentially towards zero.The level |b⟩ is irreversibly damped via its coupling with the strongly relaxing level |c⟩, which appears then as an environment for the level |b⟩.
This regime provides a tunable dissipation channel [10][11][12]41] and recently, it is widely used to simulate paritytime(PT )-symmetric Hamiltonians with postselection of the jump results [30,32,50].It can also describe the dynamics of the simplest situation of spin-spring system with relaxation processes [51].For instance the damped vacuum Rabi oscillation by the state definition |a⟩ = |g, 0⟩, |b⟩ = |e, 0⟩ and |c⟩ = |g, 1⟩, where states |a⟩ and |b⟩ are coherently coupled by the Jaynes-Cummings Hamiltonian, while |b⟩ decays towards |c⟩ at the rate γ.Further more, as we will show below that, this simple model is the core physical principle of phonon ground state cooling [5,52].

B. The Liouvillian spectra and LEPs
For this level system, the stationary state of the system is always |a⟩ no matter how the initial state and parameters of the system change.If the goal is to prepare or use this state for related applications, it is unnecessary and even harmful to wait for a long relaxation timescales.In this instance, the most practical construction is to set the optimal parameters to ensure the approaching the stationary state on a timescale which is as short as possible.It can be obtained quantitatively by solving the spectrum of Liouvillian superoperator L, as shown in Fig. 2(c), and it is and it highlights three regimes for the dynamics as a function of Ω/γ as shown in Fig. 2(c).When Ω < γ/2, L exhibits a real spectrum, which means that all the excited eigenmodes exponentially decay with time, which corresponds to the OD regime.For Ω > γ/2, on the other hand, L exhibits a complex spectrum and the system still exhibits Rabi oscillations.For arbitrary initial state, it eventually is damped out with the effectively decay rate which is determined by the gap g = γ/4.This is UD regime.When Ω = γ/2, two pairs of eigenvalues and eigenvectors of Liouvillian coalesce simultaneously (see Appendix A for the exact form of the eigensystem of Liouvillian L), giving rise to two second-order and a third-order LEPs.It corresponds to a critical damping making the boundary between the OD and the UD regime [26,40].In particular, at LEP, the gap reaches the maximum value g max = γ/4 = Ω/2, corresponding to that the dynamics at the LEP situation is fastest.
In order to investigate the physical connotations of the two types of LEPs, we reduce the master equation to the non-zero matrix elements of ρ, The coherent terms couple the populations ρ bb and ρ cc to the coherence ρ bc and ρ cb , but have no contribution to the dynamics of coherence ρ bc + ρ cb , because ρbc + ρcb = −γ/2(ρ bc + ρ cb ), which is only exponentially damped dynamics with decay rate γ/2.This means that we can not characterize LEP by observe the dynamics of (ρ bc + ρ cb ).Meanwhile, as shown in Eq. ( 9)(a), the damping term, proportional to γ, does not affect ρ bb .It contributes to the decay of ρ cc and to the corresponding increase of ρ aa .The competition between the coherent coupling and the dampling term of two states {|b⟩, |c⟩} induce the third-order LEP (with the average decay rate γ/2), witch is the phase transition point of passive PT Hamiltonian [30].And their contributions to |a⟩ give rise the second-order LEP (with average decay rate γ/4 and half rotating frequency of the third-order LEP).There are only three independent variables in Eqs. ( 9): x = ρ bb , y = ρ cc and z = −i(ρ bc − ρ cb ), which describe the dynamics of the subsystem {|b⟩, |c⟩}.With these new notations, we find the dynamics, Eigenvalues of the 3 × 3 matrix in Eq. ( 10) are −(γ ∓ κ)/2 (= λ 5,6 ) and −γ/2 (= λ 7 ).The eigenvalues can be real or complex, leading to the two different regimes qualitatively analysed above and inducing the the thirdorder LEP at κ = 0.Besides that, we also derive their dynamical evolution analytically We show the dynamics of x(t), y(t), z(t) in Fig. 3.When the decay rate is weak (γ < 2Ω), the evolution is described by damped oscillation with decay rates γ/2 (subsystem) and γ/4 (full system), respectively.Obviously, increasing γ, the evolution approaching to the stationary state will become faster, which corresponds to the quantum anti-Zeno effect [53].When the decay rate γ > 2Ω, all the eigenvalues are real and the dynamics exhibits an irreversible damping.In the limit of strong decay, γ ≫ 2Ω, the relaxation time scale is determined by so that the system will experience the metastable process for a long time scale when the system appears stationary, before eventually relaxing to ρ ss = |a⟩.This means that the larger γ is, the slower the system relaxes, which is a manifestation of the quantum Zeno [50,[54][55][56][57]. Our results show that, for a dissipative system, quantum Zeno and anti-Zeno effects correspond to the dynamical phenomena with strong and weak dissipation strengths, respectively.The LEP is thus the boundary between the quantum Zeno and anti-Zeno regimes and bridges the two previously independent effects [58].
As we mentioned before, this dynamics also leads to an effective decay from the state |b⟩ to state |a⟩, with the effective decay rate The same result can be found in [41,59] by employing perturbation theory and adiabatic elimination of states |c⟩ for a weakly driven between |b⟩ ↔ |c⟩.The above analysis shows that, our dissipative three-level model can be used to engineer decay processes between state |b⟩ and |a⟩ just by tuning the Rabi frequency Ω.
C. Engineering the relaxation dynamics

Control Liouvillian dynamics through the initial state
As we discussed in last subsection, there exist two timescales of the relaxation process depending on the space spanned by the initial state.If the initial state is an arbitrary state in space {|a⟩, {|b⟩, |c⟩}, the relaxation timescale approaching to the stationary state |a⟩ is determined by λ 1 = −(γ − κ)/4, and the fastest dynamical relaxation happens at LEP (γ = 2Ω) (see Fig. 4(a)).On the other hand, if the initial state ρ in ⊆ {|b⟩, |c⟩}, and the stationary state ρss = |a⟩ for an initial random state in the full space with γ = Ω (dashed green line), γ = 2Ω (solid red line), and γ = 3Ω (dot-dashed blue line), respectively.In this original case, the approach to stationary is governed by the eigenvalue λ1, and LEP leads to an exponentially faster convergence to the steady state with the rate gLEP = γ/4 = Ω/2.(b) Distance between the time-evolved state ρ(t) and the stationary state ρss.We compare the case of an initial random state in the full space (red line) with the time evolution ensuing the initial state in the subspace of {|b⟩, |c⟩} (green line).While in the original case, the approach to stationary is governed by the eigenvalue λ1(dashed red line), the special set of initial state leads to an exponentially faster convergence to the steady state with the rate given by λ5 (dashed green line).(c) Population dynamics versus evolution time for an initial random state in the full space (green line:Tr[ρ(t)|b⟩⟨b|], red line:Tr[ρ(t)|c⟩⟨c|]), and the time scale is governed by the eigenvalue λ5 (dashed line).(d) Observable dynamics versus evolution time for an initial random state in the full space (green line: ⟨σx⟩ = Tr[ρ(t)σx], red line: ⟨σy⟩ = Tr[ρ(t)σy], dashed blue line: ⟨σz⟩ = Tr[ρ(t)σz], and the time scale are different: for σx it is governed by the eigenvalue λ7 (dashed green line) and for σy,z they are governed by the eigenvalue λ5 (dashed red line).All the y axises are in logarithmic scale and the parameters for (b-d) are γ/Ω = 3.We have to mention here that the real dynamics and exponential decay function do not coincide at short times.This is because at short time, the decay rate is determined by all decay modes while at long time it decays with time exponentially.
as shown in Fig. 4(b), the relaxation timescale is determined by λ 5 = −(γ−κ)/2.This means that we can speed up relaxation in the convergence to stationarity by engineering the initial state, which is, the so-called Mpemba effect [1,2,17].We can understand this by looking the coefficients a i of the initial state decomposition into the left eigenmatrices L i .It can be shown that the coefficients of subspace {|b⟩, |c⟩} decomposition into L 1∼4 are all vanished, i.e., a 1∼4 = Tr[L 1∼4 ρ in ] = 0 (see Appendix A for further details).In this case therefore the asymptotic decay rate is −Re[λ 5 ] = (γ − κ)/2, which can get g LEP = γ/2.In Fig. 4(b), we compare the timescales for different initial states.It presents that if the initial state is in the full space, the approach to the stationary state is governed by the eigenvalue λ 1 (red dashed line), while the initial state in the subspace leads to an exponentially faster relaxation to the stationary state with the rate given by λ 5 = −(γ − κ)/2 (green dashed line).
In addition to the dependence of the relaxation rate on the initial state, we also find that the observable vales have significant effects on the relaxation rate (see Fig. 4(c,d)).For instance, because Tr[|i⟩⟨i|R 1∼4 ] = 0(i = a, b, c), the dynamics of the state populations approaching to stationarity is governed by the eigenvalue λ 5 (Fig. 4(c)).Moreover, the eigenmatrices R 5,6 describe the decay of σ y,z in the subspace {|b⟩, |c⟩} with rates Re[λ 5,6 ] = (γ ∓ κ)/2, while R 8 (λ 8 ) is associated with σ x (see Appendix A for further details).Whereas there occurs only damped dynamics in the subspace spanned by operators σ x , the oscillatory evolution at frequency |Im[λ 5 ]| in the subspace spanned by vectors σ y,z allows to identify the third-order LEP [26,33,36].Considering that the final state ρ ss = |a⟩ is independent on the parameters of the system, we can speed up the relaxation process by combining the acceleration effect of LEP and the initial state generation.

Tuning the Liouvillian gap through Floquet modulation
We find that the Liouvillian gap can be further increased under time-periodic (Floquet) dissipation with dissipation rate γ given by Here n ∈ Z, T = 2π/ω is the period of the Liouvillian i.e., L(t + T ) = L(t) with ω the modulation frequency, and τ is the off-duty time interval with no decay in each cycle.The density matrix at any time t is determined by the time-evolution operator P(t) = Texp( t 0 )L(t ′ )dt ′ , wihere T is the time-ordering operator.
In analogy to the case of a non-Hermitian Hamiltonian system [56,60], we can now formally define a Floquet generator for our case, an effective time-independent generator L F such that P(T ) = exp(L F T ) [61][62][63][64] transition.Here, λ P ± denote two eigenvalues with bifurcation structure, and µ = 0 marks the two eigenvalues are complex conjugate and the system is in the UD regime, while µ > 0 are in the OD regime (see Fig. 5(a)).We can see the Floquet method enriches the phase diagram.In contrast with the static dissipation (ω = 0), where the phase transition and LEPs appears at Ω/γ = 2, phase transitions under time-periodic dissipation depends on the modulation frequency ω and can occur at vanishing small dissipation strength.
Beyond that, we are more interested in the effect of modulation on the energy gap.As shown in Fig. 5(b-c), Floquet method increase the gap and the maximal gap appears at the LEP which is a different point with the static case.In the case of static dissipation (ω = 0), the g max = g LEP = Ω/2 when γ/Ω = 2.The gap under time-periodic dissipation depends on the modulation frequency ω and can even be significantly increased to bigger than Ω (see Fig. 5(c)).Fig. 5(d) plots the dynamics of the population of state |b⟩.We compare the two different cases, ω = 0 (dashed lines) and ω = Ω (the solid lines).These results illustrate that the LEP gap can surpass its static limit through Floquet engineering and thus further accelerates the relaxation.

IV. APPLICATIONS IN GROUND STATE COOLING OF TRAPPED IONS
In the following, we demonstrate the power of this approach with a practical application, i.e., the ground state cooling of trapped ions.Through analytical and numerical analysis, we will illustrate that optimal cooling conditions in the sideband and EIT approaches can be obtained, which agree with existing experiments.Our LEP gap condition provides a new perspective on optimal cooling conditions and may simulate more studies for a wide range of quantum engineering applications.

A. Sideband cooling
As shown in Fig. 6, we consider laser-ion interactions in the Lamb-Dicke limit.Dynamics is governed by Hamiltonian and jump operator J = √ γ|g⟩⟨e|.γ is the linewidth of the state |e⟩, which is coupled to state |g⟩ by a cooling laser field of frequency ω l , Rabi frequency Ω g , and detuning ∆ = ω ge − ω l , where ω ge is the frequency of the bare atomic transition |e⟩ ↔ |g⟩.ν is the trap frequency and a(a † ) is annihilation (creation) operator of phonons.Ω = ηΩ g describes the effective coupling between the phonon and internal state and η is the Lamb-Dicke parameter.When ∆ ≃ ν, the red sideband transition is nearly resonant, and the non-resonant transitions, i.e., the carrier transition and blue sideband transition, will induce the ac Stark shift to |e⟩(|g⟩) by δ(−δ), respectively.This is a good approximation that just considering shift caused by carrier transition, and under this approximation we get δ = ( Ω 2 g + ∆ 2 − ∆)/2.In order to obtain the optimal cooling condition, we reduce the overall dynamics to a low-dimensional subsystem {|g⟩|1⟩, |g⟩|0⟩, |e⟩|1⟩, |e⟩|0⟩} to obtain analytical results (see Appendix B for details).It is very helpful for us to understand the whole cooling process.Based on the perturbative calculations for this finite systems, we get that λ 1(3) = −(γ ∓ κ ′ )/4 + i(∆ + 2δ + ν)/2 and the gap In Fig. 7 (ad), we plot the analytical results of the spectra.It is obvious that the LEPs can only occur under the condition ∆ + 2δ − ν = 0.With δ = ( Ω 2 g + ∆ 2 − ∆)/2, we obtain the condition to generate LEP, Under this condition, the eigenvalues λ 1(3) become λ 1(3) = −(γ ∓ κ)/4 + iν, whose real parts are the same with three-level dissipative system shown in Fig. 2(b) and the imaginary parts connote rotating.The physical mechanism underlying condition (18) is that the detuning ∆ needs to be adjusted according to the ac Stark shift of the atomic levels to ensure that the red sideband transition is exactly on resonance.Under this premise, the level structure shown in Fig. 6 can be considered as a simple three-level dissipative system discussed in section III.When Ω ⩾ γ/2, Re[λ 1 ] = Re[λ 3 ] = −γ/4, and we get the maximum value g max = γ/4.Particularly, when In Fig. 7 (e-f), we compare the gap given by Eq. ( 17) (the solied lines) with the numerical results calculated from the full master equation (the dashed lines).Although our analytical results about spectrum and gap are obtained from the subsystem of the sideband cooling, they match very well with the numerical results calculated from the full-system.As shown in Fig. 8, the real parts of the eigenvalues λ i⩾1 , which give the relaxation rates of all the decay modes of the system, mainly can be divided into several characteristic intervals.It approximately could be {−γ/4, −γ/2, −3γ/4, −γ} under the condition (18) (see Fig. 8 (a-b)), otherwise, it becomes to be {0, −γ/2, −γ} when Ω/γ > 1/2 (see Fig. 8 (c)).As shown in Fig. 8 (c-d), different with the 3level dissipation system, it features a so-called metastable regime either Ω/γ > 1/2 or Ω/γ < 1/2, which occurs when low lying eigenvalues become separated from the rest of the spectrum [14] .The imaginary parts of the eigenvalues, which give the rotating rates of the decay modes approximately equal to Im[λ i⩾1 ] ≈ (∆+2δ+nν)/2 are mainly determined by phonon energy.
In Fig. 9(a), we plot the gap g as functions of Ω/γ and ∆/ν by using the full master equation.The point A corresponds to the LEP, and the dashed white line is the condition of g = g max that combines red sideband transition resonant condition in Eq. ( 18) as we discussed before.The numerical results and the analytical results are matched very well.Fig. 9(b) shows the the dynamics of the full system for some set of parameters (the points A, B, C, D in Fig. 9(a)).It indicates that the gap g provides a good description of the cooling time.And at the LEP, the system not only reaches stationary state at a significantly faster pace, but also obtains a lower phonon number (see Fig. 9(c-d)).Therefore, we believe that the best cooling effect can be obtained by the system parameters at the LEP.

B. EIT cooling
For the EIT cooling method discussed in [5,7,9], we find that the optimal parameter selection could be explained by the gap at LEP.As shown in Fig. 10 (  If we turn the detuning frequency ∆ g of Ω g closing to ∆ r , then |+⟩ = sin ϕ|e⟩ + cos ϕ|r⟩ (tan ϕ = Ω r /( Ω 2 r + ∆ 2 r + ∆ r )) is chosen to replace |e⟩ in the sideband cooling model (see Fig. 10 (c)).Here we just replace ), ( 19) Then the optimal condition of g max shown in Eq. ( 18) will become to If we ignore the higher order terms O(η 2 ), Eq. ( 22) can be rewritten as It is in accordance with the generalized cooling condition given in experiment works of Ref. [7,9], thus our LEP method offers a fresh understanding of optical cooling condition.
V. DISCUSSION In summary, we have studied how to engineer relaxation dynamics of Markovian open quantum systems with an arbitrary initial state.Our results have shown, for an arbitrary initial state, a speed-up relaxation can be achieved by setting the parameters of the system at LEP, where the slowest decay mode degenerates with a faster decay mode.In addition, our LEP-based accelerated approach can also be applied to accelerate the relaxation to stationarity in Floquet dissipative quantum dynamics.We have shown the relaxation processes can be dramatically faster than the static case by periodically modulating the dissipation strength.Finally, We have demonstrated the applications of our method for speeding up cooling processes in ground state cooling of trapped ions.In a broader view, our ideas may be still instructive for optimal parameter options to accelerate the cooling process even for simultaneous cooling of multiple phonon modes in a ion crystal.Therefor, our method would be in general and would facilitate the optical parameters setting in the experiments with open quantum many-body systems.Together with well-developed techniques of engineering quantum states, our work provides a powerful tool for exploring and utilizing true quantum LEP effects as examples of engineered relaxation dynamics [66].
Appendix A: Liouvillian spectrum of a dissipative three-level system We consider a simple dissipative three-level system of Fig. 2(b), and the dynamics is described by a Lindblad master equation with the Hamiltonian H = Ω/2(|c⟩⟨b| + |b⟩⟨c|) and a jump operator J = √ γ|a⟩⟨c|.
To study the Liouvilian spectra and LEPs, we first represent the Liouvillian superoperator L in a matrix form by recasting the above master equation as a matrix differential equation for the vectorized state of the density operator ρ.With the definitions that |a⟩ = (1, 0, 0) T , |b⟩ = (0, 1, 0) T , |c⟩ = (0, 0, 1) T , the Liouvillian superoperator is given by The eigenvalues of L are with κ = γ 2 − 4Ω 2 .Both the right and left eigenmatrices of the Liouvillian superoperators can be constructed to be Hermitian, they are  We are interested with the low lying eigenvalues, especially λ 1 which determines the spectral gap g = 1 4 (γ − κ ′ ).

Figure 1 .
Figure 1.(a)For an open quantum system with an arbitrary initial state, the timescale of approaching to the final stationary state ρss is related to the slowest decay mode (with eigenvalue λ1) of Liouvillian superoperator.By tuning the parameter of the system near the LEPs, where both of the slowest decaying mode and the corresponding eigenvalue are merged with a faster decaying mode and the corresponding eigenvalue (for instance λ2), the system dynamics approaches the stationary state in a much faster way.(b) This feature is evident from the Liouvillian spectrum.The stationary state ρss is characterized by the largest eigenvalue λ0 = 0.The other eigenvalues, characterizing the decay modes, have non-positive real part and always appear as complex conjugates.The Liouvillian spectral gap (g = −Re[λ1]) determines the relaxation timescale, and can reach its maximum value at LEP. (c-d) The gap at LEP can be further increased by the Floquet method (red line).In contrast with the static case (blue line), the gap under time-periodic modulation can be significantly increased, which means that relaxation process will be accelerated by applying the Floquet method.

Figure 4 .
Figure 4. (a)  Distance between the time-evolved state ρ(t) and the stationary state ρss = |a⟩ for an initial random state in the full space with γ = Ω (dashed green line), γ = 2Ω (solid red line), and γ = 3Ω (dot-dashed blue line), respectively.In this original case, the approach to stationary is governed by the eigenvalue λ1, and LEP leads to an exponentially faster convergence to the steady state with the rate gLEP = γ/4 = Ω/2.(b) Distance between the time-evolved state ρ(t) and the stationary state ρss.We compare the case of an initial random state in the full space (red line) with the time evolution ensuing the initial state in the subspace of {|b⟩, |c⟩} (green line).While in the original case, the approach to stationary is governed by the eigenvalue λ1(dashed red line), the special set of initial state leads to an exponentially faster convergence to the steady state with the rate given by λ5 (dashed green line).(c) Population dynamics versus evolution time for an initial random state in the full space (green line:Tr[ρ(t)|b⟩⟨b|], red line:Tr[ρ(t)|c⟩⟨c|]), and the time scale is governed by the eigenvalue λ5 (dashed line).(d) Observable dynamics versus evolution time for an initial random state in the full space (green line: ⟨σx⟩ = Tr[ρ(t)σx], red line: ⟨σy⟩ = Tr[ρ(t)σy], dashed blue line: ⟨σz⟩ = Tr[ρ(t)σz], and the time scale are different: for σx it is governed by the eigenvalue λ7 (dashed green line) and for σy,z they are governed by the eigenvalue λ5 (dashed red line).All the y axises are in logarithmic scale and the parameters for (b-d) are γ/Ω = 3.We have to mention here that the real dynamics and exponential decay function do not coincide at short times.This is because at short time, the decay rate is determined by all decay modes while at long time it decays with time exponentially.

Figure 6 .Figure 7 .
Figure 6.Schematic of the sideband cooling process.The cooling laser with frequency ω l drives the transition |e⟩ ↔ |g⟩ with Rabi frequency Ωg and detuning ∆ = ωeg − ω l , which leads to the ac Stark shift δ.Ω is the effective coupling strength between the red sideband transition |g⟩|n⟩ ↔ |e⟩|n − 1⟩ with n the phonon number and η the Lamb-Dicke parameter.

Figure 10 .
Figure 10.(a) Levels and transitions of the EIT cooling scheme (found in many species used for ion trapping).(b) The dressed levels what the cooling laser Ωg couples.(c) When the cooling laser is near resonant with the red sideband transition of the dressed state |+⟩, then this EIT cooling model can be equivalent to the model we discussed in example II.