A Derivation of Fokker Planck equation

# Laser cooling of high temperature oscillator by multi-level system

## Abstract

We study the laser cooling of a mechanical oscillator through the coupling with a dissipative three-level system. Under a background temperature beyond the Lamb-Dicke regime, we extend the standard cooling analysis by separately studying the classical motion and the quantum dynamics of the oscillator. In Ladder-system cooling, the cooling rate degrades by orders of magnitude at large classical motion. This phenomenon causes a critical transition of the final temperature at a hot background. In stark contrast, electromagnetic-induced-transparency (EIT) cooling with a -system produces significant negative cooling rate at high motional excitation. At steady state, the oscillator could exhibit both cooling and lasing behaviours. We argue that a successful EIT cooling requires either a poor quality oscillator to suppress the lasing effect, or terminating the cooling process at a transient stage.

## I Introduction

A macroscopic mechanical oscillator with quantum-level motion is of great interest due to its wide range of applications in establishing quantum information channels and storage (1); (2); (3), testing the foundation of quantum mechanics (4); (5); (6), measuring forces as weak as the Casmir force (7), and many others (8). By using state-of-the-art refrigeration technology, a GHz range oscillator can be cooled to the single phonon level under a background temperature mK (9). In some applications, however, a cooled oscillator with lower frequency is more desirable due to its larger zero point motion and the slower control requirement. Even at cryogenic temperature, the thermal occupation of an oscillator with frequency 100MHz is much larger than unity. Motional excitation has to be removed by applying additional cooling processes, such as feedback cooling (10); (11), sideband cooling (12); (13); (14); (15), ultrafast pulsed laser cooling (16); (17), or laser cooling with a dissipative finite-level (qudit) system (18); (19); (20); (21); (22).

The latter method is particularly favourable for quantum information processors because the dissipative qudit can be implemented by the inherent quantum memory initialization mechanism. If the initial temperature of the oscillator is sufficiently low, i.e., within the Lamb-Dicke (LD) regime, the performance of laser cooling can be analysed by using the techniques developed in atomic systems (19). In this regime, a realistic mechanical oscillator could be cooled to near the motional ground state through the coupling with a quantum dot (19), a superconducting circuit (23); (20), or a diamond NV centre (24).

In practice, however, the initial equilibrium temperature of a mechanical oscillator is often far beyond the LD regime. The high temperature performance of a two-level-system (TLS) laser cooling have been studied by computing the transition matrix of the highly excited states (18), or considering an expansion of the derivatives of the quasi-probability distribution (21); (22). The cooling rate is generally found to drop significantly at high motional excitation. Nevertheless, the final excitation predicted by the LD regime analysis is valid even beyond the LD regime.

In some architecture, the oscillator is coupled with a system, such as a diamond NV centre, that intrinsically involves more than two internal levels. The performance of laser cooling could be improved by utilising the additional states. In this article, we specifically study the the laser cooling of a high temperature oscillator with a three-level system. The objective of our work is three-fold: first, we present a strategy to study the evolution of a highly excited oscillator by separating its classical motion from its quantum dynamics. Our method is an extension of the LD regime analysis, and assumes only the regularity of the quasi-probability distribution but not the magnitude of its derivatives. Second, TLS cooling is efficient only when the system parameters are appropriate. In the low temperature regime, desired effective parameters can be engineered by appropriately driving a multi-level system, such as a Ladder-system (25). Our work shows that in spite of the discrepancy of cooling rate in the high temperature regime, the final excitation produced by a Ladder-system is close to that of an effective TLS when cooling is efficient. Third, quantum coherence between multiple metastable states can introduce cooling effects that behave differently from TLS cooling. Particularly, a lower final excitation can be achieved by using electromagnetic-induced-transparency (EIT) cooling with a -system (26); (27); (28); (29); (30). In the high temperature regime, we find that EIT cooling generally exhibits both cooling and phonon-lasing effects. Modifications of the cooling process would be required for EIT cooling to be efficient.

Our article is outlined as follows. In Sec. II, we briefly review the theoretical analysis of laser cooling in the LD regime for future comparison. In Sec. III.1, we derive an effective Fokker-Planck equation to describe the dynamics of the Glauber P distribution during laser cooling under a high temperature background. In Secs. IV and V, we apply our tools to analysis the cooling performance of a Ladder- and a -system respectively. We conclude the article in Sec. VI. In this article, we attribute the subscript () to the quantities belonging to the qudit (oscillator mode), and we have chosen units such that .

## Ii Cooling in Lamb-Dicke regime

In this section, we briefly review the analysis of laser cooling in the LD regime. More details can be found in Refs. (31); (23); (21). We consider a mechanical oscillator is coupled with a system with multiple internal levels. By selecting appropriate system parameters, we can focus on the evolution of only one oscillation mode and a -level system (qudit). The dynamics of the total mode-qudit state, , is governed by the master equation,

 ˙ρ=Lqρ+Lintρ+L0aρ+LDaρ . (1)

The first part accounts for the unperturbed evolution of the qudit, , where is the Hamiltonian of the qudit. accounts for the unperturbed evolution of the oscillation mode with the annihilation operator and the oscillation frequency . accounts for the interaction between the mode and the qudit. The lowest order mode-qudit interaction is usually linear to the mode operators, i.e.,

 Hint=λV(^a+^a†) , (2)

where is an operator acting on only the qudit levels. The effective Lamb-Dicke parameter is defined as the ratio of the interaction strength to the mode frequency, i.e., (19). and account for the open-system dynamics of the qudit and oscillation mode respectively. In particular,

 LDaρ=γ(Nth+1)D[^a]ρ+γNthD[^a†]ρ , (3)

where for some operator ; the bare damping rate is defined from the Q factor, i.e., , where is the number of cycle to lose of motional energy to a zero temperature background. is the mean phonon number of the mode if the oscillator is in equilibrium with a background with temperature .

The LD regime considers the scenario of a weak interaction strength, i.e., , and the equilibrium phonon number satisfying the criterion

 η√Nth≪1 . (4)

In this regime, the mode-qudit interaction described by is always a perturbative effect. In the master equation Eq. (1), the unperturbed dynamics is governed by ; the interaction is a first order effect; and the mode dissipation is assumed to be in the second order.

At the leading order, standard laser cooling analysis considers the system could be represented by a separable state . The qudit steady state can be obtained self-consistently by solving

 ˙ρss=Lqρss−i[λV(αss+α∗ss),ρss]=0 . (5)

The steady state displacement, , induced by the qudit steady state follows the relation

 (iν+γ/2)αss+iλ⟨V⟩ss=0 , (6)

where . The steady state contribution can be subtracted off by considering the displaced mode , where .

By adiabatically eliminating the qudit contribution, the quantum dynamics of the mode state follows

 ˙~ρa = L0a~ρa+LDa~ρa (7) Missing or unrecognized delimiter for \Big

where , and . Due to limited time correlation, the upper bound of the integration is usually replaced by . After applying Markov and rotating wave approximations, Eq. (7) becomes

 ˙~ρa = L0a~ρa+LDa~ρa−iλ2Im(S(ν)+S(−ν))[^a†^a,~ρa] (8) Missing or unrecognized delimiter for \Big

The spectral function is given by

 S(ν)=∫∞0Tr{δVUq(t−t′)δVρss}eiν(t−t′)dt′ , (9)

where is the evolution operator with respect to , i.e., . The spectral function depends on only the qudit state, and can be calculated by quantum regression theorem (32).

By using Eq. (8), the mean phonon number of the mode, , follows

 ⟨˙~n⟩=−Γc⟨~n⟩+γN , (10)

where is the cooling rate, and is the heating rate. When the mode reaches a steady state, i.e., , the final excitation of LD regime laser cooling is hence . Here we have assumed that the steady state displacement is removed by coherent operation after cooling, otherwise the final excitation is modified by a high order factor as .

## Iii Cooling in high temperature regime

In the rest of this article, we retain the assumption of weak interaction strength, i.e., , but consider background temperatures that may exceed the LD regime, i.e., . In this situation, the interaction in Eq. (2) cannot be treated perturbatively. We observe that in laser cooling literatures, the mode state can always be viewed as an ensemble of coherent state with different displacement, i.e.,

 ρa(t)=∫P(α,α∗,t)|α⟩⟨α|d2α , (11)

where the Glauber P function is positive and smooth. For this kind of states, the total state evolution can be obtained by studying the collective dynamics of each differently displaced sub-ensembles.

Our strategy based on this intuition is illustrated in Fig. 1. A high temperature ensemble is equivalent to a collection of sub-ensembles with large displacement, i.e.,

 ρ(t)=∫pc(α,α∗,t)D(α)~ρα(t)D†(α)d2α , (12)

where is the classical probability of finding a sub-ensemble that is displaced by , and is the mode-qudit state of that sub-ensemble. The dynamics of each sub-ensemble obeys Eq. (1). The mode-qudit interaction would generate a displacement-dependent force that in-turn affects the classical motion of the displacement. If the mode-qudit interaction is perturbative after subtracting the classical contribution, its effect on the mode state can be analysed by modifying the approach in the LD regime. Combining both effects for all sub-ensembles, we deduce a Fokker-Planck equation to describe the evolution of the Glauber P function of the total mode state.

A sufficient validity of our formalism is the possibility to express the total state as Eq. (12) and each sub-ensemble exhibits only perturbative mode-qudit interaction, i.e., . This condition is self-consistently satisfied in our regime of interest where and the Glauber P function is a probability of coherent states.

### iii.1 Classical Displacement Damping

Here we just outline the key steps of our derivation; a detailed treatment is presented in Appendix A. We first note that Eq. (12) is over-determined that we have the freedom to choose the implicit time dependence of . The choice determines the evolution of and ; a proper choice of could simplify the problem. Here we require to satisfy

 ˙α(t)+(iν+γ/2)α(t)+iλ⟨Vα(t)⟩=0 , (13)

such that the relative displacement of the sub-ensemble vanishes, i.e., .

With this choice, can be recognized as the classical displacement of the sub-ensemble. Eq. (13) is hence the classical equation of motion of the sub-ensemble mode state. Particularly, accounts for the damping effect due to the mode-qudit interaction. At the leading order of , the sub-ensemble qudit state, , follows

 ˙~ρα,q=Lq~ρα,q−iλ(α(t)+α∗(t))[V,~ρα,q] . (14)

When comparing to the equilibrium condition in the LD regime, Eqs. (5) and (6), here the classical motion induces a time-dependent potential on the qudit. Solving Eqs. (13) and (14) self-consistently is possible, but we further simplify the problem by considering that the displacement is dominated by the unperturbed harmonic motion, i.e., , where is defined as the degree of motional excitation. Then the motion-induced potential becomes periodic with the frequency .

In analogy to the equilibrium assumption in the LD regime, we consider the qudit state converges to the dynamic steady state (33), which is defined as

 ρss(t)=∞∑n=−∞ρ(n)sseinνt , %where˙ρ(n)ss=0 . (15)

By inserting Eq. (15) into Eq. (14), all coefficients can be obtained with the Floquet method (34).

With the dynamic steady state, we can obtain the dynamic steady value of the qudit-induced damping effect in Eq. (13), i.e., and . Particularly, the effect is dominated by the near-resonant term . As we will see, this is also the main contribution of the cooling rate for the total mode state.

### iii.2 Quantum Dynamics Around Large Displacement

After subtracting the contribution of the classical motion, the mode-qudit interaction is perturbative, and so the quantum dynamics of the sub-ensemble can be analysed by modifying the techniques in the LD regime. Our intuition is that if the mode state is initially localised within a small region in the phase space, for a sufficiently short time the coupling with the qudit and the environment will not project the state far from the region. In analogy to Eqs. (8) and (9), we find that the sub-ensemble mode state, , evolves as

 ˙~ρα,a=L0a~ρα,a+LDa~ρα,a (16) −λ2(Sα(ν,t)[^a+^a†,^a~ρα,a]+Sα(−ν,t)[^a+^a†,^a†~ρα,a] −S∗α(−ν,t)[^a+^a†,~ρα,a^a]−S∗α(ν,t)[^a+^a†,~ρα,a^a†]) ,

where the time dependent spectral function is

 Sα(±ν,t) (17) = ∫t0Tr{δVα(t)Uq(t−t′)δVα(t′)ρss(t′)}e±iν(t−t′)dt′ .

We show in Appendix A.2 the techniques for obtaining the dynamic steady value of .

### iii.3 Fokker-Planck Equation

The time evolution of the total mode state involves two parts: the redistribution of the classical displacement probability in Eq. (12), and the quantum evolution of the sub-ensemble mode states. After combining both effects through the steps in Appendix A.3, the total mode state evolution can be described by a Fokker-Planck equation of the Glauber P function,

 ˙P = ∂∂α((iνα+γ2α+iλ⟨Vα(t)⟩)P)−λ2∂2∂α2Sα(−ν,t)P (18) +12∂2∂α∂α∗(γNth+2λ2ReSα(−ν,t))P+h.c.

Eq. (18) can be simplified by making the rotating wave approximation with respect to . Furthermore, if the oscillator state is initially symmetric against the phase of , and the initial excitation is much larger than , the P function would vary with respect to only . The Fokker-Planck equation for this component of P function is

 ˙P(r,t)=1r∂∂r(r22Γc(r)P(r)+r4∂∂rγN(r)P(r)) , (19)

where the collective cooling rate is , and the collective heating rate is ; the dynamic steady spectral functions are defined as and .

When , the steady state P function is given by

 Peq(r)=AγN(r)exp(−∫r02r′Γc(r′)γN(r′)dr′) , (20)

where the normalization constant guarantees . The final excitation is hence . We note that if the steady state displacement is neglected, the final excitation has to be modified as .

We apply the tools developed in the previous section to study Ladder-system cooling in the high temperature regime. Specifically, we consider the following Hamiltonian and dissipative dynamics:

 Hq = Δ1σee+(Δ1+Δ2)σdd (21) +Ω12(σeg+σge)+Ω22(σed+σde) , V = σee , (22) LDqρ = Γ1D[σgd]ρ+Γ2D[σed]ρ , (23)

where for qudit states and . The layout of the qudit levels is shown schematically in Fig. 2(a). For the purpose of cooling, the drive between and is red-detuned from the transition frequency, i.e., . In each cycle of Ladder-system cooling, a phonon is absorbed during the transition from to . Then the excitation at is restored to through an excitation to the fast-decaying state .

The cooling rates produced by two typical Ladder-systems are shown in Fig. 3. In general, the cooling rate is high when the degree of motional excitation, , is small. As increases to the order of , the cooling rate is reduced by orders of magnitude, and eventually dominated by the bare damping rate .

The final excitation of Ladder-system cooling is shown in Fig. 4. For a wide range of , which could exceed the LD criterion in Eq. (4), the final excitation matches the prediction of LD regime analysis, . Above some critical background temperature, cooling becomes inefficient and the final excitation rises until reaching . As a side note, the critical transition is not determined by the initial excitation of the mode, because the steady state given by Eq. (20) is unique.

The criticality of the final excitation is also observed in TLS cooling (21). We attribute the criticality to the reduction of cooling rate at large . To illustrate our idea, we study a toy model that the cooling rate is a step function, i.e., and for some ‘jump’ point . is by construction the upper-bound of of a Ladder-system. If a final excitation transition appears in the toy model, so does it happens in the realistic case. For simplicity, we assume the heating rate is constant as , which is a good approximation for a high temperature background or a mode with a low Q factor.

The analytical expression of the toy model final excitation is

 nf=n2LD+e−r2c/nLD(n+(n++r2c)−nLD(nLD+r2c))nLD% +e−r2c/nLDn+ , (24)

where ; for a sufficiently large . As shown in Fig. 4, the final excitation of the toy model exhibits a transition. Specifically, if , the second parts in both the denominator and the numerator are suppressed by an exponentially small factor . Therefore, the final excitation can be approximated by , which matches the prediction of the LD regime analysis. On the other hand, if , the second parts become dominant and the final excitation is . In this regime cooling becomes inefficient.

The critical value of in Ladder-system or TLS cooling is closely related to the system parameters. Nevertheless, we can estimate the value of by considering at which the LD model becomes invalid. We consider the interaction picture with respect to the motion-induced potential; the state in this picture is transformed as . After applying the Jacobi-Anger expansion (35), Eq. (14) becomes

 ˙ρ′α,q=−i[H′q,ρ′α,q]+LDqρ′α,q , (25)

where

 H′q = Δ1σee+(Δ1+Δ2)σdd (26) +Ω12∞∑n=−∞Jn(2ηr)(einνtσeg+e−inνtσge) +Ω22∞∑n=−∞Jn(2ηr)(einνtσed+e−inνtσde) ;

and are Bessel functions of the first kind. The transformation changes only the unperturbed Hamiltonian of the qudit, but not the qudit-induced damping effect in Eq. (13), i.e., .

When comparing with Eq. (21), the drives in the original picture are effectively suppressed by a factor of . At the same time, additional sideband drives are induced with a magnitude proportional to . We can claim that the cooling mechanism in the LD regime is no longer a proper description of the system dynamics when , at where the magnitude of the original drives is significantly reduced and the sideband drives become non-negligible. Combining with the intuition obtained from the toy model, the final excitation matches the prediction of the LD regime analysis only if,

 η√nLD≪1 , (27)

which is a weaker condition than Eq. (4). Eq. (27) can also be understood as a self-consistent condition: at the steady state the mode-qudit interaction has to be a perturbative factor so that the LD regime analysis is valid.

Finally, we discuss the quality of the effective TLS engineering. In the LD regime, an effective TLS can be engineered by appropriately driving a Ladder-system and adiabatically eliminating the fast-decaying state (25). This matches our observation when is small: as shown in Fig. 3, the cooling rates of both Ladder-system are close to that of the TLS with the same effective system parameters. On the other hand, the Ladder-system cooling rate is much smaller in magnitude than the TLS cooling rate when is large. We attribute this effect to the shift of state level by the classical motion-induced potential, so that the transition becomes off-resonant. This suppresses the restoration of the qudit state to , and hence reduces the cooling rate. We also note that, at any a Ladder-system with a larger acts as a better approximation to the effective TLS. This is because the off-resonance of transition is less important if the line-width of is large.

The cooling rate discrepancy can only be experienced by a sub-ensemble if it is at large . Therefore, the final excitation can be affected only if the total mode state consists of a significant probability of highly excited sub-ensembles, in other words, only if cooling is inefficient. The idea can be observed in Fig. 4 that, the final excitation of Ladder-system and TLS cooling mismatch only when transits from to . Nevertheless, our main interest is the cases that cooling is efficient. Then the final excitation is mainly determined by the cooling rate at the small regime. In this regime, the cooling rate, and hence the final excitation, produced by the Ladder-systems and TLS are close. Our conclusion is that even though the background temperature is beyond the LD regime, if TLS cooling is found to be efficient with certain system parameters, the effective parameters can be engineered by driving a Ladder-system without compromising the cooling performance.

## V Λ system

The Hamiltonian and the dissipative dynamics of a -system are

 Hq = −Δ1σgg−Δ2σee (28) +Ω12(σgd+σdg)+Ω22(σed+σde) , V = −σgg+σee , (29) LDqρ = Γ1D[σgd]ρ+Γ2D[σed]ρ . (30)

The layout of the qudit levels is shown schematically in Fig. 2(b). -system has been extensively studied not only because it appears in various quantum systems, but it also exhibits EIT that an excitation from is dependent of the drive between other states. EIT has a variety of applications ranging from storing light to improving photon-photon interaction (36). Particularly it can be applied in laser cooling (26); (27); (28); (29); (30).

The operational condition of EIT cooling is remarkably different from TLS cooling: the metastable states are driven in blue-detuning, instead of red-detuning, from the resonant transition frequency, i.e., . Besides, efficient cooling requires quantum coherence between and . EIT cooling has several advantages over TLS cooling schemes. For examples, a lower final excitation can be attained due to the darkness of the steady state (26); (28), and multiple oscillator modes can be cooled simultaneously (29).

The EIT cooling rate with a typical -system is plotted in Fig. 5 for different . Similar to the Ladder-system and TLS cooling in the small regime, the EIT cooling rate decreases as increases. The crucial difference appears in the large regime. Instead of reducing to negligible in magnitude, the qudit contribution of the cooling rate, , remains negative for a wide range of . This is because the motion-induced potential shifts the energy levels, so the Raman transition between and is no longer in resonance. Then both and are driven in blue-detuned frequency, which is a process that would increase the motional excitation.

The regime of negative cooling rate causes dramatic effect on the steady state. In Fig. 6, we show two typical behaviours of final excitation produced by EIT cooling. Below some threshold Q factor, EIT cooling behaves similarly as Ladder-system and TLS cooling: the final excitation is close to for small , and cooling becomes inefficient when is large. On the other hand, when the Q factor is above some threshold, the final excitation is much higher than, and only slightly depending on, the background temperature. Our results contradict to the common belief that a higher Q factor would lead to a lower final excitation.

The criticality of Q factor can be understood from the two types of steady state P function that can be produced by EIT cooling. As shown in Fig. 7, for a sufficiently small Q factor the negativity of is compensated, so is positive for every . Therefore, the mode is always cooled and the steady state P function is a single peak around .

On the other hand, a large Q factor results in negative at some . The motion of the sub-ensemble will be amplified if at its displacement is negative. For this range of Q factor, in addition to the usual cooling peak at , the steady state P function consists of a large peak centred at where the cooling rate vanishes, i.e., . For a larger Q factor, the large peak contributes more significantly and eventually dominates the behaviour of the steady state. In analogy to the equilibrium photon state of a laser (37), here the large peak is a signature of the phonon lasing effect due to the blue sideband drives from and .

In conventional laser cooling, the cooling drives are applied to the qudit until the mode has reached the steady state. Following this convention, EIT cooling can reduce the motional excitation of a mode only if its Q factor is sufficiently small that the lasing effect is suppressed. Nevertheless, EIT cooling is still applicable to cool a mode with large Q factor, if the operation is halted at a transient stage before the steady state is reached.

To illustrate the idea, we simulate the mode state time evolution by numerically integrating Eq. (19). As shown in Fig. 8, the mean excitation is reduced significantly after EIT cooling is applied, and stays at the minimum for a fairly long time before it rises. The intuition here is that is positive around , at where cooling is efficient. If the initial excitation of the mode is sufficiently low, with a high probability a sub-ensemble is found within the efficient cooling regime. The emergence of phonon lasing effect requires the total mode state to diffuse, due to the heating rate , from the regime of positive to negative . If the diffusion process requires a longer time than cooling, at some transient stage the mode is efficiently cooled.

We note that the minimum transient excitation achieved by this method depends on the initial excitation. In contrast, the final temperature of conventional cooling is independent of the initial excitation, due to the uniqueness of the steady state. Our results suggest that pre-cooling is important for EIT cooling to be efficient.

## Vi Conclusion

In this article, we study the laser cooling of a mechanical oscillator with a dissipative three level system. Under a background temperature exceeding the Lamb-Dicke regime, we propose a formalism to separate the large classical motion from the quantum dynamics of the mode, and calculate the total cooling and heating rate in a self-consistent way.

In Ladder-system cooling, the cooling rate generally reduces by orders of magnitude when the displacement reaches . We show that the reduction of cooling rate causes a critical transition of the steady state excitation, which indicates that cooling is no longer efficient when the background temperature is beyond some critical value. We provide a simple criterion to estimate the regime of efficient cooling. In this regime, a Ladder-system can be used to engineer the system parameters of an effective two-level system without compromising the cooling performance.

When EIT cooling is operated with a -system, its high temperature performance is remarkably different from the Ladder-system and the TLS cooling. For a wide range of displacement, the EIT cooling rate could be negative, so that the steady state simultaneously exhibits both cooling and lasing effects. We suggest that successful EIT cooling would require either a small Q factor to suppress the lasing effect, or terminating the operation at a transient stage. In the latter case, the minimum transient temperature is dependent of the initial excitation of the mode.

We end this article with a discussion about possible strategies for better cooling performance. In TLS-like cooling scheme, the remaining question is to minimise the steady state excitation at the LD regime. This could be done by exploring the optimal system parameters and the internal level layout, e.g. Ladder-system or V-system. On the other hand, optimising EIT cooling would be more complicated because the minimum excitation state may be transient. Possible strategies includes applying TLS-like cooling to pre-cool the mode before applying EIT cooling, or dynamically varying the system parameters so that the majority of sub-ensembles experience positive cooling rate throughout the process. Overall, our work provides intuitions and tools for designing and analysing laser cooling schemes for high temperature oscillators.

H.-K. L. would like to acknowledge the support from the Croucher Foundation. This work was supported by an Alexander von Humboldt Professorship.

## Appendix A Derivation of Fokker Planck equation

This appendix provides details of the tedious steps that enter the derivation of the Fokker-Planck equation, Eq. (19), for the radial P function.

### a.1 Master equation in harmonic reference frame

We first derive Eq. (16) that describes the quantum dynamics of a sub-ensemble with classical displacement . If the displacement satisfies Eq. (13), then the sub-ensemble mode-qudit state obeys

 ˙~ρα=Lq~ρα+L1(t)~ρα+~Lint(t)~ρα+L0a~ρα+LDa~ρα , (31)

which is the same as Eq. (1) except the addition of a time dependent potential induced by the classical displacement, i.e., , and the mode-qudit interaction is replaced by , where . If the unperturbed evolution of the qudit is much faster than that of the mode and the mode-qudit interaction, then at the leading order the qudit is disentangled from the mode state and its dynamics follows Eq. (14).

We then assume the classical motion is approximately harmonic, . The assumption is valid in our case because the bare damping rate and the mode-qudit interaction strength are much smaller than . Here we have neglected the phase of the displacement for simplicity; it can be trivially added back by redefining the time reference of the oscillating phase.

Under the harmonic assumption, the motion-induced potential becomes periodic with . We then conduct Floquet analysis to obtain the dynamic steady state of the qudit in Eq. (15) (34). Following the same idea as the LD regime analysis, our aim is to adiabatically eliminate the contribution of the qudit and to construct an equation for the mode state only.

Specifically, we define the projection operators

 PX(t)≡Trq{X(t)}⊗ρss(t) , (32)

and for any density operator in the mode-qudit Hilbert space. After applying the projection onto Eq. (31), we have

 Trq{˙~ρ}⊗ρss(t)=(L0a+LDa)P~ρ+PLint(t)Q~ρ (33) ˙Q~ρ=(Lq+L1(t)+L0a+LDa)Q~ρ+QLint(t)P~ρ , (34)

where we have used the identities , , , , and .

In analogy to Eq. (7), each sub-ensemble mode state follows

 ˙~ρα,a(t) = L0a~ρα,a(t)+LDa~ρα,a(t) (35) +Trq{∫t0~L% int(t)Uq(t−t′)Ua(t−t′) ×~Lint(t′)(~ρα,a(t′)⊗ρss(t′))dt′} ,

where the evolution operators are defined as for satisfying Eq. (14); as is assumed to be at the second order of .

Before moving forward, we note several important differences between Eqs. (7) and (35). First, the time independent steady quantities in the LD regime are substituted by the dynamic steady quantities that are varying as time. Therefore the time ordering in the integral becomes important. Second, the upper bound of the integration cannot be replaced by . As we will see, the dependence in the integration bound is necessary for obtaining the correct evolution of the spectral function.

In analogy to the derivation of Eq. (8), we consider the evolution of the mode is dominated by the free evolution, i.e., , then Eq. (35) can be written as Eq. (16).

Solving Eq. (16) requires the dynamic steady values of the qudit state and the spectral function . Here we discuss the techniques for computing these values.

The dynamics of the qudit state following Eq. (14) can be solved by considering the evolution of the generalized bloch vector, , where is a complete set of operator in the qudit state space. In our three-level system analysis, we pick . The Bloch vector version of Eq. (14) is given by

 ˙⟨σ⟩=(M−iλ(α(t)+α∗(t))V)⋅⟨σ⟩ , (36)

where the time independent matrix elements are defined as , and .

After making the harmonic approximation on the displacement, finding the dynamic steady state is equivalent to finding the solution of Eq. (36) that satisfies

 ⟨σ(t)⟩=∞∑n=−∞→σ(n)einνt , (37)

where .

Eq. (36) is not easy to solve because is singular. The singularity originates from the conservation of the trace, . Such condition can be removed from the equation by transforming the Bloch vector as , and then considering only the second to ninth entries. We specifically consider is a 9x9 matrix that corresponds to the transformation , , , and otherwise. After the transformation, we have

 ⟨˙~σ⟩=(~M−iλr(e−iνt+eiνt)~V)⋅⟨~σ⟩+u , (38)

where , , and . The truncation matrix is defined as , where is a 8x1 zero vector; is a x identity matrix.

Substituting the dynamic steady Bloch vector in Eq. (37) into Eq. (38) and matching the components with the same frequency, we can see the solution obeys

 inν→σ(n)=~M⋅→σ(n)−iλr~V⋅(→σ(n+1)+→σ(n−1))+δn,0u , (39)

where is the Kronecker delta. By assuming vanishes for a large enough , the dynamic steady , and hence , can be obtained efficiently by the continued fraction method (21).

The spectral function can be solved in a similar way. We first note that the spectral function in Eq. (17) is identical to

 Sα(±ν,t)=∫t0Tr{VUq(t−t′)δVα(t′)ρss(t′)}e±iν(t−t′)dt′ . (40)

We then define the spectral vector as

 S(±ν,t)=∫t0Tr{σUq(t−t′)δVα(t′)ρss(t′)}e±iν(t−t′)dt′ . (41)

If the spectral vector is obtained, the spectral function can be obtained by picking the specific components, i.e., for Ladder-system, and for -system.

Differentiating Eq. (41) with respect to and applying the quantum regression theorem (32), we have

 ˙S(±ν,t) = (±iνI9+M−iλ(α(t)+α∗(t))V)⋅S(±ν,t) (42) +Tr{σδVα(t)ρss(t)} .

We note that the last term, which is essential to recover the LD regime result at the limit , would be missing if the integral upper bound in Eq. (17) is taken as .

Under the harmonic approximation of the displacement, the dynamic steady spectral function can be obtained from the solution of Eq. (42) that satisfies

 S(±ν,t)=∞∑n=−∞→Sn(±ν)einνt, (43)

and . Nevertheless, solving Eq. (42) also suffers from the problem of the singularity of . We again extract the singular component by the transformation . In the transformed basis, Eq. (44) can be written as

 ˙~S(±ν,t) = Missing or unrecognized delimiter for \big (44) ⋅~S(±ν,t)+Tr{~σδVα(t)ρss(t)} .

Substituting Eq. (43) into Eq. (44), the dynamic steady solution obeys

 inν→Sn(±ν) = (±iνI8+~M)⋅→Sn(±ν) (45) −iλr~V⋅(→Sn+1(±ν)+→Sn−1(±ν))+→vn ,

where . By assuming vanishes for a sufficiently large , the solution can be obtained by the continued fraction method (21).

Before moving forward, we note that the above method can reproduce the LD regime steady state and spectral function by setting and considering only the non-oscillating component, i.e., .

### a.3 Fokker-Planck equation

Differentiating Eq. (12) and tracing out the qudit contribution, we get

 ˙ρa(t) = ∫˙pc(α,α∗,t)D(α)~ρα