title: Performance of superadiabatic stimulated Raman adiabatic passage in the presence of dissipation and Ornstein-Uhlenbeck dephasing
authors: Kostas Blekos Dionisis Stefanatos Emmanuel Paspalakis
text: 
  
  
  
  
  
  
  
  
  
  dionisis@post.harvard.edu
  
  
  
  Materials Science Department, School of Natural Sciences, University of Patras, Patras 26504, Greece
  
  
  
  
               
  
  
  
  In this paper we evaluate the performance of two superadiabatic stimulated Raman adiabatic passage (STIRAP) protocols derived from Gaussian and sin-cos pulses, under dissipation and Ornstein-Uhlenbeck noise in the energy levels. We find that for small amplitudes of Stokes and pump pulses, the population transfer is mainly achieved directly through the counterdiabatic pulse, while for large amplitudes the conventional STIRAP path dominates. This kind of “hedging" leads to a remarkable robustness against dissipation in the lossy intermediate state. For small pulse amplitudes and increasing noise correlation time the performance is decreased, since the dominant counterdiabatic pulse is affected more, while for large pulse amplitudes, where the STIRAP path dominates, the efficiency is degraded more for intermediate correlation times (compared to the pulse duration). For the Gaussian superadiabatic STIRAP protocol we also investigate the effect of delay between pump and Stokes pulses and find that under the presence of noise the performance is improved for increasing delay. We conclude that the Gaussian protocol with suitably chosen delay and the sin-cos protocol perform quite well even under severe noise conditions. The present work is expected to have a broad spectrum of applications, since STIRAP has a crucial role in modern quantum technology.
  
  
  
  
                               
  
                                
  Performance of superadiabatic stimulated Raman adiabatic passage in the presence of dissipation and Ornstein-Uhlenbeck dephasing
      Emmanuel Paspalakis
      December 30, 2023
  ================================================================================================================================
  
  
  
  
  
  
  § INTRODUCTION
  
  
  One of the most successful quantum control methods for population transfer between the levels of a quantum system is Stimulated Raman adiabatic passage (STIRAP) <cit.>. The prototype STIRAP system consists of three energy levels in the Λ-configuration. In order to transfer population from state |1⟩ to state |3⟩, passing through the intermediate level |2⟩, two laser pulses are applied in counterintuitive order, the Stokes pulse coupling states |2⟩-|3⟩ and the pump pulse coupling states |1⟩-|2⟩. A coherent superposition is formed by states |1⟩ and |3⟩, which adiabatically evolves from state |1⟩ initially to state |3⟩ finally, while the lossy intermediate state |2⟩ is barely populated. STIRAP finds a wide range of applications in modern quantum science, from optical wavequides <cit.> and matter waves <cit.> to nitrogen-vacancy centers in diamond <cit.> and superconducting quantum circuits <cit.>, for more details see the recently published roadmap <cit.>. The most important advantage of STIRAP is its robustness against moderate variations of experimental parameters. Its major drawback is that adiabatic transfer requires long times, leading to reduced efficiency when undesirable interactions with the environment, for example decoherence and dissipation, are present.
  
  During the last decade, a series of methods characterized as shortcuts to adiabaticity have been developed with the aim to improve the performance of quantum adiabatic evolution by reducing the required duration <cit.>. The common basic idea of these techniques is to drive faster the quantum system at the same final state as the slow adiabatic evolution. This goal is achieved either bypassing the intermediate adiabatic states, or by introducing an extra term in the Hamiltonian to suppress the diabatic transitions and evolve the system along the adiabatic path of the original Hamiltonian. The latter approach is called superadiabatic, assisted adiabatic passage or transitionless tracking algorithm. Both methods are widely exploited in modern quantum technologies, as is thoroughly discussed in the recently published review <cit.>, and have also been used to increase STIRAP efficiency, see for example Refs. <cit.> and <cit.> for the first and second method, respectively. 
  
  The influence of noise on the efficiency of STIRAP has been the subject of several works. The effect of dephasing caused by classical Ornstein-Ulhenbeck noise <cit.> in the level energies was studied in Ref. <cit.>, while the dephasing due to quantum baths was evaluated in Refs. <cit.>. In the context of superconducting artificial atoms, the influence of broadband colored noise on population transfer was considered in Refs. <cit.>. Along with these studies which focus on conventional STIRAP, there are also works which examine the performance of STIRAP shortcuts in the presence of noise. For the case where no additional counterdiabatic field is used, i.e. the classical STIRAP framework is preserved, there are some recent works in the physical context of nitrogen-vacancy centers in diamond <cit.>. Specifically, in Ref. <cit.> the efficiency of STIRAP shortcuts developed in Refs. <cit.> was investigated in the presence of dissipation and spectral diffusion, while in Ref. <cit.> the performance of a Floquet-engineered shortcut under colored noise was evaluated. In our recent work <cit.> we studied the efficiency of the shortcuts derived in <cit.>, in the presence of classical Ornstein-Ulhenbeck noise in the energy levels, as in Ref. <cit.>. There are also a few works <cit.> which study the influence of noise on superadiabatic STIRAP (SA-STIRAP, in the spirit of Giannelli and Arimondo <cit.>), where an extra counterdiabatic term is exploited in the Hamiltonian. Specifically, the authors of Ref. <cit.> consider the effect of non-Gaussian distribution of energy fluctuations, while those of Ref. <cit.> study the effect of ground state dephasing.
  
  In this article we evaluate the performance, in the presence of noise and dissipation, of two SA-STIRAP protocols derived from Gaussian and sin-cos pulses <cit.>.  We use classical Ornstein-Ulhenbeck noise processes with exponential correlation functions <cit.> in the energy levels, as in Refs. <cit.>. This type of noise provides a relatively simple means to study the effect of colored noise, with non-zero correlation time. It has been used to model the fluctuations in the energy levels for molecules in a liquid <cit.>. It may also be present in the phases of the applied laser fields, appearing in the rotating wave approximation as a corresponding noise term in the energy levels <cit.>. Also note that the sensitivity of shortcuts to adiabaticity to Ornstein-Ulhenbeck noise has also been studied in the context of fast shuttling of an atom using a moving optical lattice <cit.>. From numerical simulations, we find that for small amplitudes of Stokes and pump pulses the population transfer is mainly achieved directly from level |1⟩ to level |3⟩ through the counterdiabatic pulse, while for large amplitudes the STIRAP path |1⟩→|2⟩→|3⟩ dominates. The presence of two alternative paths leading to the target state results in a remarkable robustness against dissipation in the lossy intermediate level. For small pulse amplitudes and increasing noise correlation time the performance is decreased, since the dominant counterdiabatic pulse is affected more. For large pulse amplitudes the efficiency is degraded more for intermediate correlation times (compared to the pulse duration), while it is better for small or large noise correlation times. This behavior is similar to that observed for conventional STIRAP <cit.>, something expected since for large amplitudes STIRAP is the dominant mechanism for population transfer. For Gaussian SA-STIRAP we also investigate the effect of delay between pump and Stokes pulses and find that under the presence of noise the performance is improved as the delay increases (at least up to the considered values). We conclude that the Gaussian SA-STIRAP protocol with suitably chosen delay and the sin-cos SA-STIRAP protocol perform quite well even under severe noise conditions.
  
  The structure of the article is as follows. In the next section we derive for completeness the Gaussian and sin-cos SA-STIRAP protocols in the absence of noise. In section <ref> we study their performance when dissipation and noise are present. Section <ref> concludes this paper.
  
  
  
  § SA-STIRAP IN THE ABSENCE OF DISSIPATION AND DEPHASING
  
  
  
  
  
  
  
  
  
  
  The reference Hamiltonian for STIRAP in both one-photon and two-photon resonance is
  
      H_0(t)=ħ/2(
      [      0 Ω_p(t)      0; Ω_p(t)    -iΓ Ω_s(t);      0 Ω_s(t)      0 ]),
  
  where Ω_p(t), Ω_s(t) are the Rabi frequencies for the pump and Stokes lasers, respectively, and Γ is the dissipation rate from level |2⟩. In this section we derive SA-STIRAP for the ideal case Γ=0
  
  
  
  
  while in the next section we study the effect of nonzero dissipation and dephasing.
  
  
  
  If we define the time-dependent amplitude Ω(t) and mixing angle θ(t) through the relations
  
      Ω(t)=√(Ω_p^2(t)+Ω_s^2(t)),  tanθ(t)=Ω_p(t)/Ω_s(t),
  
  then the instantaneous eigenstates of H_0(t) are
  
      |ϕ_0(t)⟩=
      ([  cosθ;     0; -sinθ ]),
        
      |ϕ_±(t)⟩=1/√(2)([ sinθ;  ± 1; cosθ ]),
  
  with corresponding eigenvalues
  
      E_0(t)=0,   E_±(t)=±ħΩ(t)/2  .
  
  In conventional STIRAP the angle θ changes slowly (adiabatically) from the initial value θ(t_i)=0 at t=t_i to the final value θ(t_f)=π/2 at t=t_f, while the the system follows the dark adiabatic state |ϕ_0(t)⟩=cosθ(t)|1⟩-sinθ(t)|3⟩, from the initial state |1⟩ to the final state |3⟩.
  
  The adiabatic approximation fails for fast variations of angle θ, in which case some population remains in levels |1⟩ and |2⟩ at the final time. A method to accomplish the desired population transfer to level |3⟩ in arbitrarily short times is to add in Hamiltonian H_0 an extra counter-diabatic Hamiltonian H_cd, which cancels the diabatic terms arising when H_0 is transformed to the time-dependent adiabatic basis,
  
      H_cd(t)   =   iħ∑_n=0,±[|ϕ̇_n(t)⟩⟨ϕ_n(t)|
            -⟨ϕ_n(t)|ϕ̇_n(t)⟩|ϕ_n(t)⟩⟨ϕ_n(t)|].
  
  Using Eq. (<ref>) we find for the STIRAP system
  
      H_cd(t)=ħ/2(
      [     0     0  iΩ_d;     0     0     0; -iΩ_d     0     0 ]),
  
  where
  
      Ω_d(t)=2θ̇(t)
  
  is a π-pulse connecting directly states |1⟩ and |3⟩, since ∫_t_i^t_fΩ_d(t)dt=∫_t_i^t_f 2θ̇(t)dt=π.
  Under the total Hamiltonian H=H_0+H_cd the system can track with perfect fidelity and for arbitrarily short times the dark state |ϕ_0(t)⟩ of the reference Hamiltonian H_0, and for this reason the method is called transitionless tracking algorithm.
  
  Note that usually the direct transition between levels |1⟩ and |3⟩ is electric dipole forbidden. One physical system where such transition is possible is described in Refs. <cit.>. It is a Λ-system between J=1 Zeeman sublevels and an excited J=0 state. The coupling between the |1⟩ and |3⟩ states may be accomplished by a magnetic dipole interaction between the atomic or molecular angular momentum 𝐉 and an external magnetic field.  Another example, also discussed in Ref. <cit.>, contains a ladder system which includes both single-photon transitions for the two electric dipole allowed transitions and a two-photon transition for the electric dipole forbidden transition. Different systems where all transitions in the three-level system are electric dipole allowed exist in asymmetric quantum systems with broken inversion symmetry, e.g. asymmetric molecules <cit.>. We note that the applicability of the method is not restricted to atoms or molecules, but can also be exploited for population transfer in artificial atoms, like superconducting quantum structures and semiconductor nanostructures <cit.>. For example, in the recent experiment <cit.>, SA-STIRAP was used in the first three states of a superconducting transmon circuit, to transfer population between the ground state and the second excited state, by combining two single-photon transitions and a two-photon transition. The advantage of the SA-STIRAP method, compared to a direct π-pulse between states |1⟩ and |3⟩ or classical STIRAP with only pump and Stokes fields, is its increased fidelity and robustness, since it combines the useful characteristics of both methods. The increased fidelity and robustness is particularly significant for quantum technology applications.
  
  In the next section we will evaluate the performance of two SA-STIRAP protocols under noise. For the first protocol the pump and Stokes pulses are Gaussian <cit.>
  
      Ω_p(t)=Ω_0 e^-(t-τ/T)^2,  Ω_s(t)=Ω_0 e^-(t+τ/T)^2,
  
  with corresponding counterdiabatic pulse
  
      Ω_d(t)=4τ/T^2(4τ t/T^2).
  
  Note that 2τ is the delay between the Gaussian pulses and T is the pulse width.
  For the second protocol the pump and Stokes pulses have the sin-cos shape <cit.>
  
      Ω_p(t)=Ω_0 sin (π t/2T ),  Ω_s(t)=Ω_0 cos (π t/2T ),
  
  while the corresponding counterdiabatic pulse is constant,
  
      Ω_d(t)=π/T.
  
  
  
  
  
  
  § PERFORMANCE OF SA-STIRAP PROTOCOLS UNDER DISSIPATION AND ORNSTEIN-UHLENBECK DEPHASING
  
  
  
  
  We study the effect of noise by adding to the total Hamiltonian a noise term,
  
      H(t)=H_0(t)+H_cd(t)+H_ϵ(t),
  
  where
  
      H_ϵ(t)=ħ/2(
      [ ϵ_1(t)      0      0;      0 ϵ_2(t)      0;      0      0 ϵ_3(t) ])
  
  and ϵ_i(t), i=1,2,3, are independent Ornstein-Uhlenbeck noise process. The latter are defined by the stochastic differential equations <cit.>
  
      ϵ̇_i=-1/τ_cϵ_i+1/τ_cg_i,
  
  where g_i(t), i=1,2,3, are independent Gaussian white noises with zero mean and correlations
  
      ⟨ g_i(t)g_j(t+τ)⟩=2σ^2τ_cδ_ijδ(τ).
  
  The stochastic processes ϵ_i are also Gaussian with zero mean, correlations
  
      R_ij(τ)=⟨ϵ_i(t)ϵ_j(t+τ)⟩=σ^2δ_ije^-|τ|/τ_c,
  
  and steady state probability distributions
  
      P_i(ϵ)=1/√(2πσ^2)e^-ϵ^2/2σ^2,
  
  where note that parameter τ_c corresponds to the noise correlation time. The power spectral density of each Ornstein-Uhlenbeck process is
  
      S_ii(ω)=∫_-∞^∞R_ii(τ)e^-iωτdτ=2σ^2τ_c/1+(ωτ_c)^2
  
  and the expectation value of the total power is
  
      ⟨ϵ_i^2(t)⟩=R_ii(0)=1/2π∫_-∞^∞S_ii(ω)dω=σ^2.
  
  
  We implement the Ornstein-Uhlenbeck processes using the algorithm described in Ref. <cit.>, with specific standard deviation σ=5/T, corresponding to fixed noise power, and three correlation times τ_c/T=0.008, 0.08, 0.8. More details about the numerical implementation can be found in the Appendix. 
  In the following Figs. <ref>-<ref> we present simulation results for both STIRAP and SA-STIRAP using various pulses. In all these figures we plot the transfer fidelity F (average final population of level |3⟩ over 200 stochastic runs) versus the peak amplitude Ω_0 of pump and Stokes pulses, for various values of delay τ and dissipation Γ. We start with Fig. <ref>, where we use Gaussian pulses with delay parameter τ/T=1/4. In Figs. <ref>, <ref>, <ref>, <ref> we display results for the conventional STIRAP and four dissipation values, Γ=0, 1/T, 4/T, 10/T, respectively. Solid blue line corresponds to the case where there is no dephasing noise, while the other three lines correspond to different noise correlation times: τ_c/T=0.008, (orange dashed line), τ_c/T=0.08 (green dotted line), and τ_c/T=0.8 (red dashed-dotted line). For the cases where dephasing noise is present, we also show the 98% confidence interval, i.e. the range of values where lies the fidelity of 98% of the stochastic runs. In Figs. <ref>, <ref>, <ref>, <ref> we display similar results but for the corresponding SA-STIRAP protocol (with the additional counterdiabatic pulse), also for Γ=0, 1/T, 4/T, 10/T, respectively. In Figs. <ref>, <ref>, <ref> we show analogous results for Gaussian protocols with delay parameters τ/T=1/3, 1/2, 3/4, respectively, while in Fig. <ref> for the sin-cos protocol.
  
  Several interesting observations can be made from these figures. First note that in the absence of noise or for small noise correlation time the STIRAP fidelity shows an oscillatory behavior for smaller delay values τ/T, Figs. <ref> and <ref>, which disappears as τ/T increases, Figs. <ref> and <ref>. As discussed in Ref. <cit.>, these fidelity maxima occur when the Rabi oscillations between the levels are synchronized with the interaction time. Also observe that, although the first maximum in Figs. <ref> and <ref> appears at values Ω_0 lower than those needed in Figs. <ref> and <ref>, high levels of fidelity are maintained only for a narrow window around this maximum. An oscillatory fidelity is also observed for STIRAP with sin-cos pulses in the absence of dephasing noise or for noise with small τ_c, see Fig. <ref>. The STIRAP efficiency is reduced in the presence of dephasing noise, as well as for nonzero dissipation, requiring larger values of Ω_0, see the first column of Figs. <ref>-<ref>.
  
  We now move to discuss the performance of SA-STIRAP protocols. The most important observation when comparing the first (STIRAP) and second (SA-STIRAP) columns of Figs. <ref> to <ref> is that the presence of the counterdiabatic pulse leads to nozero fidelity even for small values of Ω_0. In the absence of dephasing noise a perfect transfer is achieved for every value of Ω_0, while the efficiency is reduced when noise is present. From these plots becomes obvious that for SA-STIRAP and in the case of small Ω_0 the population is mainly transferred directly from state |1⟩ to |3⟩ using the counterdiabatic pulse, while for larger Ω_0 dominates the STIRAP path passing through the intermediate state |2⟩. This is also demonstrated in Fig. <ref>, where the final total population P in all three levels is plotted versus Ω_0, for Gaussian SA-STIRAP with delay τ/T=1/2 and dissipation rate Γ=4/T. Observe that, as Ω_0 increases, there appear population losses (P<1), indicating that part of the population passes through the lossy level |2⟩. For lower Ω_0 the transfer to level |3⟩ cannot be completed and most of the  population in the intermediate state |2⟩ is dissipated. As Ω_0 increases the STIRAP path is established and the transfer to state |3⟩ is completed, as indicated by the convergence of total population P to the fidelity F when comparing Figs. <ref> and <ref>.
  
  We emphasize that SA-STIRAP is different than the shortcut method used in our recent work <cit.> since, in addition to the pump and Stokes pulses, it exploits the extra field Ω_d connecting directly states |1⟩ and |3⟩. The advantage of using the additional field in the presence of dissipation is demonstrated in Fig. <ref>, where note that for simplicity we ignore dephasing noise. In Fig. <ref> we display the shortcut pump and Stokes pulses, derived from Gaussian profiles, that we used in Ref. <cit.> (see Fig. 2 there), while in Fig. <ref> we plot the corresponding time evolution of the third level population p_3, for different values of dissipation rate Γ=0, 1/T, 4/T, 10/T (from top to bottom). Observe that, as the dissipation increases, the transfer fidelity at the final time drops considerably. The reason is that the transfer is accomplished through the lossy intermediate state |2⟩. In Fig. <ref> we plot the Gaussian SA-STIRAP pump and Stokes pulses, as well as the counterdiabatic pulse Ω_d (red solid line), for Ω_0=5/T and τ/T=3/4, while in Fig. <ref> the corresponding time evolution of population p_3, for the same dissipation values as before, Γ=0, 1/T, 4/T, 10/T. Observe that now the effect of dissipation can be hardly distinguished. The reason is that the transfer is mostly accomplished directly from |1⟩ to |3⟩ by the extra field Ω_d. For comparison we mention that the common area of SA-STIRAP pump and Stokes pulses is 8.8558 units, while those of the shortcut pulses are 10.3563 and 12.3597 units. By inspecting the second columns of Figs. <ref>-<ref>, it becomes evident that SA-STIRAP shows a remarkable robustness against dissipation Γ, even in the presence of dephasing noise.
  
  
  Another interesting observation in the case of Gaussian SA-STIRAP is that for fixed delay τ between the Stokes and pump pulses, the fidelity for small Ω_0 drops with increasing noise correlation time τ_c, as is evident in the second columns of Figs. <ref>-<ref>. This behavior can be explained as follows. Recall that in the range of small Ω_0 the desired transfer is achieved through the counterdiabatic pulse given in Eq. (<ref>) which, for fixed delay τ, is specified. The Fourier transform of this pulse is
  
      F_d(ω)=∫_-∞^∞Ω_d(t)e^-iω tdt=π(π T^2/8τω),
  
  from which it is obvious that the pulse spectrum is mainly concentrated between 0 and 8τ/(π T^2). Now observe from Eq. (<ref>) that, as the correlation time τ_c increases, the noise power is concentrated in lower frequencies around zero, affecting thus more the baseband counterdiabatic pulse. This is also demonstrated in Fig. <ref>. The sin-cos protocol also exhibits a similar behavior for increasing τ_c, see the second column of Fig. <ref>. A closely related observation is that for fixed τ_c, Γ and increasing τ, i.e. when examining the same type of curve in subfigures belonging to the same row of the second column across Figs. <ref> to <ref> (f.e. green dotted curve corresponding to τ_c/T=0.08 across Figs. <ref> to <ref>, where Γ=0), the fidelity for small Ω_0 increases. The explanation is that, as the delay τ increases, the spectrum of the counterdiabatic pulse broadens and thus is less affected by the baseband noise, as shown in Fig. <ref>.
  
  For large Ω_0 it is evident from Figs. <ref>-<ref> that the fidelity for dephasing noise with large correlation time τ_c/T=0.8 (red dashed-dotted line) is similar or in most cases better than that for noise with intermediate correlation time τ_c/T=0.08 (green dotted line). This behavior has been observed for STIRAP in Ref. <cit.> and is confirmed here also for SA-STIRAP since, as we pointed out above, for large Ω_0 the STIRAP path dominates. We next try to understand this behavior using Fig. <ref>, where we employ Gaussian STIRAP with delay τ/T=1/2, dissipation Γ=0 and no dephasing noise, thus the reference for fidelity is Fig. <ref>. In Fig. <ref> we plot the time evolution of level |2⟩ population p_2 and in Fig. <ref> the trajectory on the Bloch sphere, for large Ω_0=60/T.  The oscillations of p_2, most of which occur with a period of about 0.2T, correspond to the revolution of Bloch vector (red line) around the instantaneous total field which, as the pump and Stokes pulses evolve, rotates from the north pole to the equator (on the plane of the meridian shown). Dephasing noise with small correlation time τ_c/T=0.008 is averaged during an oscillation period and its net effect is small, as pointed out in <cit.>. For large τ_c/T=0.8, the noise remains constant for large parts of the evolution, and its influence is reduced since Ω_0=60/T>σ=5/T.  It is the noise with intermediate correlation time τ_c/T=0.08, comparable with the period of oscillations, which affects the most the desired transfer, and this is clearly demonstrated in Fig. <ref> for Ω_0=60/T. In Figs. <ref>, <ref>, we plot the population p_2 and the trajectory on the Bloch sphere, respectively, for lower Ω_0=15/T. Observe that here the period of oscillations is comparable to τ_c/T=0.8, consequently the noise with this correlation time is now more effective, as depicted in Fig. <ref> for Ω_0=15/T.
  
  Closing, we would like to emphasize that, by inspecting fidelity in Figs. <ref> to <ref>, we conclude that the best performance under severe noise conditions is observed for Gaussian SA-STIRAP with delay τ/T=3/4 and the sin-cos SA-STIRAP protocol.
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  § CONCLUSION
  
  
  
  
  We investigated the performance of two SA-STIRAP schemes obtained from Gaussian and sin-cos pulses, in the presence of dissipation and dephasing by exponentially correlated noise. For small amplitudes of Stokes and pump pulses we found that the population transfer is mainly accomplished directly by the counterdiabatic pulse, while for large amplitudes conventional STIRAP dominates. This kind of “hedging" leads to a remarkable robustness against dissipation in the lossy intermediate level. For small pulse amplitudes and increasing noise correlation time we found a decreasing performance, since the dominant counterdiabatic pulse is affected more, while for large pulse amplitudes, where the STIRAP path dominates, the efficiency is degraded more for intermediate correlation times (compared to the pulse duration). We also studied for the Gaussian SA-STIRAP the effect of delay between pump and Stokes pulses and found that in the presence of noise the performance is improved for increasing delay. Our conclusion is that the Gaussian SA-STIRAP protocol with suitably chosen delay and the sin-cos SA-STIRAP protocol perform quite well even under severe noise conditions. The current work is expected to find application in the implementation of emerging quantum technologies, since STIRAP is a widely used method in this research area.
  
  *
  
  
  
  § NOTES ON NUMERICAL METHOD
  
  
  
  
  The simplest type of colored noise is that with exponential correlation function.
  It requires only one more parameter than white noise while being a more
  realistic noise source in the context of laser noise.
  In general, a stochastic differential equation
  with an exponentially correlated colored noise term
  can be modeled as
  
      ẋ_i    = f(x_i) + ϵ_i 
      ϵ̇_i    = -1/τ_cϵ_i +1/τ_cg_i
  
  where g_i is simple zero-mean Gaussian white noise.
  
  
  
  
  
  
  To integrate this stochastic system
  an explicit fixed time step method is preferable
  to control the stochastic dynamics.
  We use an explicit 4th order Runge-Kutta method
  to integrate the first, deterministic part of the system
  (equation <ref>)
  coupled with
  the algorithm by Fox et al. <cit.>
  to integrate for the noise (Eq. <ref>).
  A sufficiently small time step must be used to ensure
  both that the noise terms are properly integrated
  and the deterministic part of the system is accurately
  simulated.
  In general we have used a time step Δ t that is at most
  1/10 the value of correlation time τ_c.
  
  Algorithm <ref> shows the steps
  of the integration of Eq. <ref>.
  First, an initial noise value ϵ_i|_t=0 is generated
  and the intensity-related parameter E is calculated for the
  integration interval Δ t and correlation time τ_c.
  At each time step of the integration,
  a new ϵ_i value is calculated,
  based on the previous value of ϵ_i
  and a Gaussian parameter h of mean zero and standard deviation σ;
  then, ϵ_i are a series of values that
  realize a  zero-mean Ornstein-Uhlenbeck process with correlation time τ_c,
  standard deviation σ, obeying the correlation relation
  ⟨ϵ_i(t)ϵ_i(t+τ) ⟩ = σ^2 e^-|τ|/τ_c,
  as discussed in the previous sections.
  
  In Figs. <ref> and <ref>
  we compare the numerical results to the
  theoretically desired properties of the noise.
  The numerical data are averaged over 10 iterations
  of the integration procedure.
  Fig. <ref> shows
  the relative frequency of the generated ϵ_i values.
  Fig. <ref>
  shows the mean correlation vs the lag-time τ.
  It is obvious that
  the numerical results are in excellent agreement
  with the theoretically desired properties.
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  99
  
  Bergmann98
  K. Bergmann, H. Theuer, and B. W. Shore, Rev. Mod. Phys. 70, 1003 (1998).
  
  Kobrak98
  M. N. Kobrak and S. A. Rice, Phys. Rev. A 57, 2885 (1998).
  
  Vitanov01
  P. Král, I. Thanopulos, and M. Shapiro, Rev. Mod. Phys. 79, 53 (2007).
  
  Vitanov17
  N. V. Vitanov, A. A. Rangelov, B. W. Shore, and K. Bergmann, Rev. Mod. Phys. 89, 015006 (2017).
  
  Bergmann19
  K. Bergmann et al., J. Phys. B: At. Mol. Opt. Phys. 52, 202001 (2019).
  
  
  
  
  Dreisow09
  F. Dreisow, A. Szameit, M. Heinrich M, R. Keil, S. Nolte, A. Tünnermann, and S. Longhi, Opt. Lett. 34 2405 (2009).
  
  Menchon16
  R. Menchon-Enrich, A. Benseny, V. Ahufinger, A. D. Greentree, T. Busch and J. Mompart, Rep. Prog. Phys. 79 074401 (2016).
  
  Golter14
  D. A. Golter and H. L. Wang, Phys. Rev. Lett. 112, 116403 (2014).
  
  Kumar16
  K. S. Kumar, A. Vepsäläinen A, S. Danilin, and G. S. Paraoanu, Nat. Commun. 7 10628 (2016).
  
  
  
  Odelin19
  D. Guéry-Odelin, A. Ruschhaupt, A. Kiely, E. Torrontegui, S. Martínez-Garaot, and J. G. Muga, Rev. Mod. Phys. 91, 045001 (2019).
  
  Unayan97
  R. G. Unanyan, L. P. Yatsenko, K. Bergmann, and B. W. Shore, Opt. Commun. 139, 48 (1997).
  
  Demirplak03
  M. Demirplak and S. A. Rice, J. Phys. Chem. A 107, 9937 (2003).
  
  Berry09
  M. V. Berry, J. Phys. A: Math. Theor. 42, 365303 (2009).
  
  Motzoi09
  F. Motzoi, J. M. Gambetta, P. Rebentrost, and F. K. Wilhelm, Phys. Rev. Lett. 103, 110501 (2009).
  
  Chen10a
  X. Chen, A. Ruschhaupt, S. Schmidt, A. del Campo, D. Guéry-Odelin, and J. G. Muga, Phys. Rev. Lett. 104, 063002 (2010).
  
  Masuda10
  S. Masuda and K. Nakamura Proc. R. Soc. A 466, 1135 (2010).
  
  Deffner14
  S. Deffner, C. Jarzynski C and A. del Campo, Phys. Rev. X 4, 021013 (2014).
  
  Claeys19
  P. W. Claeys, M. Pandey, D. Sels, and A. Polkovnikov, Phys. Rev. Lett. 123, 090602 (2019).
  
  
  
  Li16
  Y.-C. Li and Xi Chen, Phys. Rev. A 94, 063411 (2016).
  
  Clerk16
  A. Baksic, H. Ribeiro, and A.A. Clerk, Phys. Rev. Lett. 116, 230503 (2016).
  
  Du16
  Y-X. Du, Z.-T. Liang, Y.-C. Li, X.-X. Yue, Q.-X. Lv, W. Huang, X. Chen, Hui Yan, and S.-L. Zhu, Nat. Commun. 7, 12479 (2016).
  
  Zhou17
  B. B. Zhou, A. Baksic, H. Ribeiro, C. G. Yale, F. J. Heremans, P. C. Jerger, A. Auer, G. Burkard, A. A. Clerk, and D. D. Awschalom, Nat. Phys. 13 330 (2017).
  
  Mortensen18
  H. L. Mortensen, J. J. Sørensen, K. Mølmer, and J. F. Sherson, New J. Phys. 20 025009 (2018).
  
  
  
  
  Demirplak05
  M. Demirplak and S. A. Rice, J. Phys. Chem. B 2005, 109, 6838 (2005).
  
  Chen10b
  X. Chen, I. Lizuain, A. Ruschhaupt, D. Guéry-Odelin, and J. G. Muga, Phys. Rev. Lett. 105, 123003 (2010).
  
  Giannelli14
  L. Giannelli and E. Arimondo, Phys. Rev. A 89, 033419 (2014).
  
  Torosov14
  B. T. Torosov, G. Della Valle, and S. Longhi, Phys. Rev. A 89, 063412 (2014).
  
  Masuda15b
  S. Masuda and S. A. Rice, J. Phys. Chem. A 2015, 119, 3479 (2015).
  
  Masuda15
  S. Masuda and S. A. Rice, J. Chem. Phys. 142, 244303 (2015).
  
  Vepsalainen19
  A. Vepsäläinen A, S. Danilin, and G. S. Paraoanu, Science Adv. 5, eaau5999 (2019).
  
  Impens19
  F. Impens and D. Guéry-Odelin, Scient. Rep. 9, 4048 (2019).
  
  
  
  
  
  
  
  
  
  Fox88
  R. F. Fox, I. R. Gatland, R. Roy and G. Vemuri, Phys. Rev. A 38 5938 (1988).
  
  Demirplak02
  M. Demirplak and S. A. Rice, J. Chem. Phys. 116, 8028 (2002).
  
  Shi03
  Q. Shi and E. Geva, J. Chem. Phys. 119, 11773 (2003).
  
  Vitanov04
  P. A. Ivanov, N. V. Vitanov, and K. Bergmann, Phys. Rev. A 70, 063409 (2004).
  
  Zeng19
  Y. X. Zeng, T. Gebremariam, M. S. Ding, and C. Li, Ann. Phys. (Berlin) 531, 1800234 (2019).
  
  
  
  
  Falci13
  G. Falci, A. La Cognata, M. Berritta, A. D'Arrigo, E. Paladino, and B. Spagnolo, Phys. Rev. B 87, 214515 (2013).
  
  DiStefano16
  P. G. Di Stefano, E. Paladino, T. J. Pope, and G. Falci, Phys. Rev. A 93, 051801(R) (2016).
  
  
  Boyers19
  E. Boyers, M. Pandey, D. K. Campbell, A. Polkovnikov, D. Sels, and A. O. Sushkov, Phys. Rev. A 100, 012341 (2019).
  
  Ibanez12
  S. Ibáñ́ez, X. Chen, E. Torrontegui, J. G. Muga, and A. Ruschhaupt, Phys. Rev. Lett. 109, 100403 (2012).
  
  Stefanatos20
  D. Stefanatos, K. Blekos, and E. Paspalakis, Appl. Sci. 10, 1580 (2020).
  
  Issoufa14
  Y. H. Issoufa and A. Messikh, Phys. Rev. A 90, 055402 (2014).
  
  Laine96
  T. A. Laine and S. Stenholm, Phys. Rev. A 53, 2501 (1996).
  
  Chen12
  X. Chen and J. G. Muga, Phys. Rev. A 86, 033405 (2012).
  
  Lu20
  X.-J. Lu, A. Ruschhaupt, S. Martínez-Garaot, and J. G. Muga, Entropy 22, 262 (2020).
  
  Kral03a P. Král, I. Thanopulos, M. Shapiro, and D. Cohen,
  Phys. Rev. Lett. 90, 033001 (2003).
  
  Kral01a P. Král and M. Shapiro, Phys. Rev. Lett. 87, 183002 (2001).
  
  Nori05a Y.-X. Liu, J. Q. You, L. F. Wei, C. P. Sun, and F. Nori,  Phys. Rev. Lett. 95, 087001 (2005).
  
  
  
  