Liouvillian of the Open STIRAP Problem

Abstract

With the corresponding Liouvillian as a starting point, we demonstrate two seemingly new phenomena of the STIRAP problem when subjected to irreversible losses. It is argued that both of these can be understood from an underlying Zeno effect, and in particular both can be viewed as if the environment assists the STIRAP population transfer. The first of these is found for relative strong dephasing, and, in the language of the Liouvillian, it is explained from the explicit form of the matrix generating the time-evolution; the coherence terms of the state decay off, which prohibits further population transfer. For pure dissipation, another Zeno effect is found, where the presence of a non-zero Liouvillian gap protects the system’s (adiabatic) state from non-adiabatic excitations. In contrast to full Zeno freezing of the evolution, which is often found in many problems without explicit time-dependence, here, the freezing takes place in the adiabatic basis such that the system still evolves but adiabatically.

1. Introduction

Coherent control has become an essential part of many branches in quantum physics, ranging from atomic and molecular thermal gases to solid state devises and ultracold atomic condensates [1,2,3,4]. Utilising adiabatic driving is a method to circumvent any errors arising from timing of applied pulses—as long as the process is adiabatic, the target state is reached. When a direct coupling between the initial and target state is forbidden by selection rules, an intermediate state is often used in a so called Raman coupling scheme. As often is the case, this additional state is an excited state that suffers from spontaneous emission. This source of dissipation/decoherence can be bypassed by using instead the so called STIRAP—stimulated Raman adiabatic passage method [1,5,6]. In addition, the STIRAP makes use of an intermediate state, but by properly choosing the two Raman coupling pulses, the coherent transfer of population can be made perfect without ever populating the excited intermediate state. As such, deficiencies due to spontaneous emission of that state are greatly suppressed [7,8]. However, dephasing of the two lower states may well occur, for example due to elastic collisions between particles [9]. The process is then no longer coherent, which indeed affects the population transfer. While the evolution is no longer unitary, it is still linear and generated by a Liouvillian operator. Naturally, the properties of the Liouvillian will determine the evolution of the system’s state.

In AMO experiments, the coupling to any environment can usually be made rather weak and especially when working in the optical regime the system dynamics can be well approximated by a Markovian master equation. The general form of such an equation is given by the Lindblad type [10]

$${\partial }_t\widehat{\rho }=\widehat{\mathcal{L}}\left[\widehat{\rho }\right]=i[\widehat{\rho },\widehat{H}]+\mathcal{D}\left[\widehat{\rho }\right],$$

(1)

with the disipator

$${\displaystyle \mathcal{D}\left[\widehat{\rho }\right]=\sum _i{\gamma }_i\left(2{\widehat{L}}_i\widehat{\rho }{\widehat{L}}_i^{†}-{\widehat{L}}_i^{†}{\widehat{L}}_i\widehat{\rho }-\widehat{\rho }{\widehat{L}}_i^{†}{\widehat{L}}_i\right).}$$

(2)

Here, we have defined the Liouvillian $\widehat{\mathcal{L}},$ which generates the full evolution of the state, the unitary time-evolution governed by the system Hamiltonian $\widehat{H}$, and the part incorporates the effects of the environment. The effective system–environment couplings are ${\gamma }_i$ with corresponding Lindblad jump operators ${\widehat{L}}_i$. The Liouvillian on the form above preserves the physical properties of the state, that is $\widehat{\rho }\left(t\right)$ is normalized and positive semi-definite for all times. That implies also that the evolution is in agreement with the fluctuation-dissipation theorem stating that any sort of dissipation inevitably causes some sort of decoherence. Neglecting such fluctuations, one of the first studies on the open STIRAP consisted of adding a complex component $-i\gamma |2\rangle \langle 2|$ to the Hamiltonian [7]. In this simplified picture, it was found, via adiabatic elimination, that the transition probability $P_i$ for $i=3$ falls off exponentially in terms of the dissipation rate $\gamma$ for weak dissipation $g_0\gg \gamma$, and as ${\gamma }^{-2}$ in the opposite limit $\gamma \gg g_0$. Ref. [8] considered instead the full Lindblad master Equation (1) describing spontaneous emission of the excited intermediate level. For a large regime of decay rates and being deep in the adiabatic regime, their results showed that the success rate was in principle unaltered by the openness of the problem. In a more recent work [11], the open STIRAP was reexamined by deriving a microscopic master equation that takes into account the explicit time-dependence of the Hamiltonian (In the standard derivation of the Lindblad master equation, the Hamiltonian is assumed time-independent. The derivation cannot be performed following the same procedure once the Hamiltonian becomes time-dependent.). In the overdamped regime, it was found that an increased population transfer is to be found, which results from a Zeno-type process. Dephasing between the two states $|1\rangle$ and $|3\rangle$ is expected to have a greater effect on the STIRAP process. In particular, dephasing among the lower states will effectively couple the different adiabatic states, and hence result in deterioration of the population transfer success rate. Influence of dephasing on the STIRAP problem has been analysed in Ref. [9] in the weak system–environment coupling limit.

In this work, instead of directly integrating the Lindblad master Equation (1), which has been the standard route in the past, we start from the properties of the Liouvillian $\widehat{\mathcal{L}}$ in order to analyse the general open STIRAP problem. We present new results that in a natural way can be explained in terms of the Liouvillian. In particular, by extracting the (complex) spectrum of $\widehat{\mathcal{L}}$, the different relevant time-scales can more easily be identified. For example, one may expect that, if the STIRAP process is slow, the environment will have a greater influence on the system dynamics, and possibly deteriorate the success rate. Simultaneously, a too fast STIRAP will imply non-adiabatic excitations taking you out from the desired instantaneous eigenstate. This expected behaviour is indeed found for weak couplings to the environment, and the optimal time for the process should be such that the inherent unitary and the external irreversible time scales agree. For strong dephasing, a surprising evolution is demonstrated. Here, it might actually be preferable to consider a slow process despite the coupling to the environment. It is explained from a Zeno freezing effect of the population transfer between the different diabatic states. Another Zeno manifestation is found for dissipation of the excited state. This is understood from the non-vanishing Liouvillian gap that tends to project the state back onto the desirable adiabatic state. There is a trade-off though, a too strong coupling to the environment implies a qualitatively different regime where the quantum Zeno mechanism projects you down to a diabatic state instead of the adiabatic state.

The paper is structured as follows. In the next section, we start by recapitulating the general idea behind STIRAP. The following subsection introduces the Bloch representation of the Lindblad master equation and we briefly discuss some needed properties of the Liouvillian matrix. We continue in Section 3 with the actual numerical results. In particular, the spectrum of the Liouvillian matrix is analysed in some detail. Furthermore, the full open STIRAP problem is simulated numerically, which confirms our predictions drawn from the structure of the Liouvillian matrix. Finally, we conclude in Section 4 with a summary and some remarks. In the Appendices, we especially point out some remarks about the Liouvillian matrix.

2. The Open STIRAP Model

The STIRAP model is a paradigm, both for demonstrating adiabatic ideas in general and as an example of a time-dependent three-level system. In this respect, it is a perfect starting point when analysing adiabatic processess from the perspective of the Liouvillian $\widehat{\mathcal{L}}\left(t\right)$. Indeed, the corresponding Liouvillian encapsulates many novel features, such as exceptional points (EP’s) [12,13]. The Liouvillians for open two-level systems, such as the Landau–Zener problem [14,15], seem to be lacking many of these interesting aspects.

STIRAP for Closed Systems

The standard STIRAP setup is the $\Lambda$ one depicted in Figure 1; two stable states, $|1\rangle$ and $|3\rangle$ respectively, are laser coupled by a pump $G_1\left(t\right)$ and a Stokes field $G_2\left(t\right)$ to an excited intermediate state $|2\rangle$. Throughout, we assume a two-photon resonance transition, i.e., the photon frequency difference $\hslash ({\omega }_2-{\omega }_1)$ (with ${\omega }_{1,2}$ the frequency of the pump and Stoke lasers respectively) matches the energy difference between the bare states $|1\rangle$ and $|3\rangle$. For the main part of this analytical subsection, we do not restrict the analysis to single photon resonance transitions, and thereby the introduction of a detuning $\Delta$ [1]. It is found, however, that the qualitative results will not depend on $\Delta$ and, for all numerical simulations, we thereby let $\Delta =0$ for simplicity (see, for example, Ref. [16]). Before analysing the corresponding Liouvillian, we next summarise the basics of the STIRAP mechanism.

Within the rotating wave approximation and using the bare (diabatic) basis introduced above, the Schrödinger equation becomes [1] ($\hslash =1$ throughout)

$$i{\partial }_t|\mathsf{\Psi }\left(t\right)\rangle ={\widehat{H}}_{\mathrm{d}}\left(t\right)|\mathsf{\Psi }\left(t\right)\rangle =\left[\begin{array}{ccc}0& G_1\left(t\right)& 0\\ G_1\left(t\right)& \Delta & G_2\left(t\right)\\ 0& G_2\left(t\right)& 0\end{array}\right]|\mathsf{\Psi }\left(t\right)\rangle .$$

(3)

This equation defines the time-dependent Hamiltonian ${\widehat{H}}_{\mathrm{d}}$, and furthermore we have assumed real couplings $G_{1,2}\left(t\right)$ without loss of generality. The adiabatic basis $\left\{|{\phi }_{+}\left(t\right)\rangle ,|{\phi }_0\left(t\right)\rangle ,|{\phi }_{-}\left(t\right)\rangle \right\}$ is given by the instantaneous (adiabatic) eigenstates of ${\widehat{H}}_{\mathrm{d}}\left(t\right)$. Written as a unitary matrix given in the diabatic basis, these states are

$$\widehat{U}\left(t\right)=\left[\begin{array}{ccc}sin\varphi sin\theta & cos\theta & cos\varphi sin\theta \\ cos\varphi & 0& -sin\theta \\ sin\varphi cos\theta & -sin\theta & cos\varphi cos\theta \end{array}\right],$$

(4)

with the parametrisation $tan\theta =G_1\left(t\right)/G_2\left(t\right)$, $tan2\varphi =2G_0\left(t\right)/\Delta$, and with $G_0^2\left(t\right)=G_1^2\left(t\right)+G_2^2\left(t\right)$. The corresponding instantaneous (adiabatic) eigenstates are $E_{\pm }\left(t\right)=\left(\Delta \pm \sqrt{{\Delta }^2+4G_0^2\left(t\right)}\right)/2$ and $E_0\left(t\right)=0$. Thus, we have ${\widehat{H}}_{\mathrm{ad}}\left(t\right)\equiv \widehat{U}{\widehat{H}}_{\mathrm{d}}\left(t\right){\widehat{U}}^{-1}=\mathrm{diag}\left(E_{+}\left(t\right),E_0\left(t\right),E_{-}\left(t\right)\right)$, with ${\widehat{H}}_{\mathrm{ad}}\left(t\right)$ the adiabatic Hamiltonian. The Schrödinger equation in the adiabatic basis reads

$$i{\partial }_t|{\mathsf{\Psi }}_{\mathrm{ad}}\left(t\right)\rangle =\left[{\widehat{H}}_{\mathrm{ad}}\left(t\right)-i{\widehat{U}}^{-1}{\partial }_t\widehat{U}\right]|{\mathsf{\Psi }}_{\mathrm{ad}}\left(t\right)\rangle .$$

(5)

The adiabatic approximation consists of dropping the “gauge potential” $\widehat{A}=-i\widehat{U}{\partial }_t{\widehat{U}}^{-1}$ that comprises the non-adiabatic couplings [1,17,18]. Thus, within this approximation, the adiabatic states evolve as $|{\phi }_i\left(t\right)\rangle \to exp\left(-i{\int }_0^t{E}_i\left(t^{\prime }\right)d{t}^{\prime }\right)|{\phi }_i\left(t\right)\rangle$ with $i=\pm ,0$ (Here, we put $\widehat{A}\equiv 0$, but it is noted that this gauge potential is responsible for the Berry phase when encircling a closed loop in parameter space). The adiabatic state (also referred to as dark state) $|{\phi }_0\left(t\right)\rangle$ is particularly attractive for practical purposes as it does not contain the bare excited state $|2\rangle$, which is typically subject to spontaneous emission. Now, if the couplings are chosen

$$\begin{array}{c}{\displaystyle \underset{t\to -\infty }{lim}\frac{G_1\left(t\right)}{G_2\left(t\right)}=0,\theta \to 0,}\\ {\displaystyle \underset{t\to +\infty }{lim}\frac{G_2\left(t\right)}{G_1\left(t\right)}=0,\theta \to \frac{\pi }{2},}\end{array}$$

(6)

it follows that, provided that the evolution is adiabatic, the adiabatic dark state obeys

$$|{\phi }_0\left(t\right)\rangle =\left\{\begin{array}{ccc}|1\rangle ,\hfill & \hfill & t=-\infty ,\hfill \\ |3\rangle ,\hfill & \hfill & t=+\infty .\hfill \end{array}\right.$$

(7)

This defines the STIRAP in our $\Lambda$-configuration; if we prepare our state in $|1\rangle$ and adiabatically turn on the couplings according to Equation (6), we will steer the state into $|3\rangle$ without ever populating the mediating state $|2\rangle$. One simple choice of couplings obeying this condition is for two symmetric Gaussians:

$$\begin{array}{c}{\displaystyle G_1\left(t\right)=g_{0}exp\left[-\frac{{(t-a\tau )}^2}{2{\left(a\sigma \right)}^2}\right],}\\ {\displaystyle G_2\left(t\right)=g_{0}exp\left[-\frac{{(t+a\tau )}^2}{2{\left(a\sigma \right)}^2}\right].}\end{array}$$

(8)

Here, $g_0$ is the pulse amplitude, $a\sigma$ the pulse width, and $2a\tau$ the distance between the pulses. These are shown in Figure 2, and they are also the pulses used throughout this manuscript. What is especially worth pointing out is the counterintuitive order of the pulses; pulse 2, the pump, which couples the initially empty states, is turned on before pulse 1, the Stokes. We have introduced a to serve as our single “adiabaticity parameter” (qualitatively, it is the pulse area setting the level of adiabaticity [5], and hence one could imagine varying other parameters, like the pulse amplitude, instead). Thus, we will keep $g_0$, $\tau$ and $\sigma$ fixed in all numerical simulations (more precisely, $g_0=1$ and $\tau =\sigma =10$) and instead vary a alone.

To get a better insight into the parameters rendering adiabatic evolution, we make use of the criteria for adiabaticity [18]

$$\frac{|\langle {\phi }_{+}\left(t\right)|({\partial }_t{\widehat{H}}_{\mathrm{d}})|{\phi }_0\left(t\right)\rangle |}{E_{+}^2\left(t\right)}\ll 1,$$

(9)

where we have used the symmetry of the $|{\phi }_{\pm }\left(t\right)\rangle$ adiabatic states relative to $|{\phi }_0\left(t\right)\rangle$, and that $E_0\left(t\right)=0$. This results in

$$\mathcal{A}\left(t\right)\equiv \frac{G_1\left(t\right)G_2\left(t\right)}{G_0^3\left(t\right)}\frac{a\tau }{{\left(a\sigma \right)}^2}\ll 1.$$

(10)

For adiabatic evolution, condition (10) should be fulfilled for all times $t\in [t_i,t_f]$ between initial and final times. For not too large initial and final times, $\mathcal{A}\left(t\right)$ peaks at $t=0$, and we directly notice that since $\mathcal{A}\left(t\right)\sim 1/g_0$ and $\mathcal{A}\left(t\right)\sim 1/a$ a large amplitude $g_0$ and/or a large a favour adiabaticity. As we already pointed out, we fix $g_0=1$ and instead vary a in order to analyse the influence of non-adiabatic excitations. As we will see, for the open STIRAP problem, a is in general not a good adiabaticity parameter.

3. Results and Discussion

There are in particular two types of dissipation/decoherence processes occurring naturally in realistic settings, spontaneous emission of the excited state $|2\rangle$ and dephasing. Previous works have considered direct numerical integration of the master equation or adiabatic elimination schemes [7,8,9]. Here, we examine the open STIRAP problem by rewriting the Lindblad master equation as a linear differential equation for the Bloch vector of the state. In two dimensions, the resulting equations are the famous Bloch equations [19]. In higher dimension, however, much less is known about these “generalised Bloch equations”, and the interpretation and analysis of the equations become also much more complex as will be discussed below.

The system evolution is modelled within the Lindblad master equation formalism, Equation (1). We will work in the diabatic basis such that the Hamiltonian is taken as ${\widehat{H}}_{\mathrm{d}}$ of Equation (3), and we consider two different situations of loss channels:

(1)

Case (a). Dephasing of the lower states $|1\rangle$ and $|3\rangle$ implemented by the Lindblad jump operator

$$\widehat{L}=|1\rangle \langle 1|-|3\rangle \langle 3|.$$

(11)

We disregard any dephasing between the other levels as these are typically of less importance [9]. More precisely, once the dephasing arising from the jump operator (11) is taken into account, the additional dephasing occurring between the other levels do not qualitatively alter the results. This is especially true when the intermediate state is only slightly populated.

(2)

Case (b). Spontaneous emission of the excited state $|2\rangle$ to the states $|1\rangle$ and $|3\rangle$. The corresponding jump operators are

$${\widehat{L}}_1=|1\rangle \langle 2|,{\widehat{L}}_2=|3\rangle \langle 2|.$$

(12)

Note that we do not study the situation of losses to a fourth level $|4\rangle$ as was the scenario of Refs. [7,8]. Furthermore, we assume the same decay rates, ${\kappa }_1={\kappa }_2=\kappa$.

These two cases are schematically presented in Figure 3.

3.1. The Liouvillian Matrix for the $\Lambda$ System

Any state in any dimension D can be written on the form [20,21]

$$\widehat{\rho }=\frac{1}{D}\left(\mathbb{I}+\sqrt{\frac{D(D-1)}{2}}\mathbf{R}·\lambda \right),$$

(13)

where $\mathbf{R}=(r_1,r_2,\cdots ,r_{D^2-1})$ is the generalised Bloch vector and $\lambda =({\lambda }_1,{\lambda }_2,\cdots ,{\lambda }_{D^2-1})$ a vector with the generalised Gell–Mann matrices as elements [22] (see Appendix A). For pure states, one has $\left|\mathbf{R}\right|=R=1$, and naturally for the maximally mixed state $\mathbf{R}=\mathbf{0}$. Solving the master equation for $\widehat{\rho }$ now transforms into the problem of solving an equation for the Bloch vector on the form [23,24,25]

$${\partial }_t\mathbf{R}=\mathbf{M}\mathbf{R}+\mathbf{b}.$$

(14)

Here, $\mathbf{M}$ is a $(D^2-1)\times (D^2-1)$ matrix generating the time-evolution of the Bloch vector, and will be hereafter called the Liouvillian matrix, while $\mathbf{b}$ we denote Liouvillian pump for reasons to be explained. In the Appendix B, we give several general properties of $\mathbf{M}$. What is especially important is that it is not an Hermitian matrix, and thereby its eigenvalues ${\mu }_i$ are typically complex. The eigenvalues must, however, appear in complex conjugate pairs since the trace of $\mathbf{M}$ is real. Moreover, it follows that, whenever an eigenvalue is complex, the corresponding eigenvector must also be complex. As a result, since the Bloch vector is real for any physical state, most eigenvectors do actually not represent physical states $\widehat{\rho }$. Another peculiarity, when dealing with non-Hermitian matrices, is that different eigenvectors need not be mutually orthogonal. In particular, at the EPs when two (or more) eigenvalues ${\mu }_i$ become degenerate the corresponding eigenvectors are parallel [12,13]. One quantity that is of special importance is the eigenvalue with the smallest (but non-zero) absolute value of its real part, i.e.,

$$\tilde{\Delta }=\underset{i}{min}\left[\mathrm{Re}(-{\mu }_i)\right],$$

(15)

which sets an upper bound for the time-scale for reaching the steady state, and it is often referred to as the Liouvillian gap [26,27].

The Liouvillian matrices corresponding to the two cases (11) and (12), together with some of their properties, are presented in the Appendix C. As the couplings $G_1\left(t\right)$ and $G_2\left(t\right)$ now are time-dependent different scales become relevant. We will come back to this in the following subsection for the general case. In this subsection, however, we fix the values of the couplings, i.e., we consider the standard Raman $\Lambda$ model. In particular, we consider $G_1=G_2=1$ and $\Delta =0$. We have verified that varying these particular values does not alter our conclusions. With these parameters fixed, we ask how the eigenvalues of $\mathbf{M}$ vary with the system–environment couplings $\gamma$ and $\kappa$. While keeping the couplings constant presents an oversimplified picture (loosely speaking assuming that the short time-scale is that of environment induced relaxation), it does give valuable insight. In this respect, this subsection is to be seen as providing an intuition for the dynamics given some non-Hermitian matrix $\mathbf{M}$.

We may note that for $\gamma =0$ or $\kappa =0$, and since the Liouvillian matrices are skew-symmetric, we find two zero eigenvalues. The corresponding eigenstates are the steady states of the model. From the expression (A7) of the Appendix, it is easy to identify the Bloch vector

$${\mathbf{R}}_0={\left[\begin{array}{cccccccc}0& 0& \sqrt{3}{c}^2& -\sqrt{3}sc& 0& 0& 0& (1-3s^2)\end{array}\right]}^t/2$$

(16)

as one of them, with $s=sin\theta$ and $c=cos\theta$. This Bloch vector corresponds to the state with zero eigenvalue of the Raman model, which for the STIRAP system would be the adiabatic state $|{\phi }_0\left(t\right)\rangle$ introduced in the previous section. As $\gamma$ or $\kappa$ become non-zero, the two zero eigenvalues split but stay real (and negative) [28]. The Bloch vector (16) is no longer representing a steady state. With a non-zero Liouvillian gap, one expects relaxation to the maximally mixed state, which would be the unique steady state. This is also what one finds for case $\left(a\right)$ of dephasing. For case $\left(b\right)$, the Liouvillian pump $\mathbf{b}$ is non-zero and this implies that the Bloch vector is still representing steady state. One may say that the pump counteracts the relaxation such that the steady state becomes a non-trivial state. This observation is important for understanding the STIRAP problem exposed to dissipation as we discuss in the following subsection.

In Figure 4 and Figure 5, we display the real and imaginary parts of the eigenvalues of the Liouvillian matrices for case (a) and (b), respectively. In both cases, we have that the Liouvillian gap $\tilde{\Delta }$ is non-zero whenever the coupling to the environment is present. For case (a), where $\mathbf{b}=\mathbf{0}$, as we already mentioned this results in that the steady state is the maximally mixed state with $\mathbf{R}=\mathbf{0}$. For case (b), the real parts of the eigenvalues are also all non-positive, but as we pointed out the non-vanishing pump makes it possible that the steady state contains quantum coherences. The complex conjugation pairing of the eigenvalues is evident in the plots of the imaginary parts. The disappearance of imaginary parts at the EPs implies a splitting of the real parts, reminiscent of a bifurcation. At these points, the eigenvectors become real, which, however, is not enough to warrant that they represent physical states. In order for them to serve as physical states, their lengths must be relatively short. Indeed, given that the Bloch vector is real, it is always possible to construct a physical state $\widehat{\rho }$ from it given that one shrinks its length sufficiently much [29]. In case (a), for large enough $\gamma$s, all eigenstates become purely real after the three EPs. In case (b), however, a “reversed” EP takes place where two purely real eigenstates become imaginary upon increasing $\kappa$.

The EPs apparent in both spectra, Figure 4 and Figure 5, clearly demonstrate an example of non-analytic behaviour. A natural question arises whether their presence can result in visible effects in the system evolution, and in particular non-analytical behaviours. Indeed, the connection between EP’s and non-equilibrium phase transitions has been discussed in the past [30,31]. One may envision signatures of these in our model, for example the time-evolution of some initial state could qualitatively change if one varies $\gamma$ or $\kappa$; whenever a splitting of the real parts occurs, the exponents rendering the decay change that could alter the system’s evolution. We explore this by calculating the population imbalance $Z=\mathrm{Tr}\left[\left(\right|1\rangle \langle 1|-|3\rangle \langle 3\left|\right)\widehat{\rho }\left(t\right)\right]$ for an initial random pure state. The result for case (a) is shown in Figure 6. For small $\gamma$ (<$G_1,G_2$) and sufficiently short times, the evolution is dominated by unitary time-evolution, i.e., the population imbalance displays Rabi oscillations between the two states. When $\gamma$ is increased, the inverse ${\gamma }^{-1}$ determines the fast time-scale and the relaxation occurs before any Rabi oscillations have time to manifest. The presence of non-analytic features of the eigenvalues are, however, not reflected in the time-evolution of the physical states of Figure 6.

It can be shown, however, that, for non-physical states, corresponding to complex eigenstates ${\mathbf{R}}_i\left(t\right)$ non-analytic time-evolution is found upon quenching through an EP by considering a time-dependent $\gamma \left(t\right)$ [23]. In this respect, we think of the instantaneous eigenstates ${\mathbf{R}}_i\left(t\right)$ as the adiabatic states of $\mathbf{M}\left(t\right)$, and at an EP the characteristic time-scale diverges such that adiabaticity necessarily breaks down [32].

3.2. Dynamics—Numerical Results

The previous subsection analysed general effects stemming from the environment for a frozen Raman Hamiltonian, and not the actual interplay between external influence and the inherent STIRAP dynamics. This is the main focus of the present subsection.

If the Liouvillian gap $\tilde{\Delta }$ is large in comparison to the inverse time that sets the inherent time scale, Hamiltonian adiabatic evolution is not guaranteed by a large adiabaticity parameter a. Generally speaking, a large parameter a favours internal adiabatic evolution, but it implies an extended coupling to the environment, which, in return, tends to take the system out of its instantaneous adiabatic Hamiltonian eigenstate [23,24]. Consequently, one expects that there is an optimal $a_{opt}$ such that the intrinsic unitary evolution is close to adiabatic and at the same time excitations due to the environment are not too definite. This scenario should apply to case (a) representing dephasing, while for case (b) there should not be a trade-off between the two processes, i.e., a slow passage does not automatically lead to environment induced excitations since the state remains in desired adiabatic state. Recently, a similar scenario was also discussed in quantum many-body systems [33]. However, the parameter dependence of the optimal time $a_{opt}$ is not known. It is especially interesting to compare the optimal $a_{opt}$ with the rate $\gamma$. Note from Figure 4 that, for $\gamma <2,$ the real parts of the eigenvalues ${\mu }_i$ scale linearly with $\gamma$, meaning that there is a linear dependence between the characteristic relaxation time scale and the inverse of the decay rate. It is reasonable to assume that $a_{opt}$ thereby must decrease for growing $\gamma$. We will see that this is only partly true, namely for small $\gamma$.

The numerical results for the full time-dependent STIRAP problem for case (a) are displayed in Figure 7. As a measure of the efficiency of STIRAP, we calculate the population

$$P_i=\langle i|\widehat{\rho }\left(t_f\right)|i\rangle ,i=1,2,3,$$

(17)

at the final time $t_f$. For perfect transfer, all population ends up in the state $|3\rangle$, i.e., $P_3=1$. The integration interval $t\in [t_i,t_f]$ ($t_f=-t_i=100$) is chosen long enough such that convergence has been reached. The uppermost curve in plot (a) represents the closed STIRAP and we find that a serves as an adiabaticity parameter, i.e., there is a one-to-one relation between its value and the success rate of the STIRAP. In the other examples displayed in the figure, $\gamma \ne 0$ and the maximum population transfer is not found for $a\to \infty$, but, for a finite a, which defines $a_{opt}$. For the decay rates used in the figure, the convergence to $P_3=1/3$ occurs for moderate a values. In fact, for any non-zero $\gamma$, we have that $P_1=P_2=P_3=1/3$ in the limit of large a. In this limit, the state is the maximally mixed one $\widehat{\rho }=\mathbb{I}/3$, which can be understood since it is also the instantaneous steady state of the corresponding Liouvillian. Nevertheless, when $G_1=G_2=0,$ any diagonal density matrix, in the diabatic basis, is a steady state, and hence, since the STIRAP pulses go to zero for large times, one could expect an asymptotic state different from the maximally mixed one. However, for large as, the pulses are “sufficiently” non-zero to warrant relaxation towards the steady state ${\widehat{\rho }}_{ss}=\mathbb{I}/3$. The $\gamma$-dependence of the optimal $a_{opt}$ is shown in Figure 7b. For larger values of $\gamma$ beyond 2, there is no longer a clear maximum any longer. The plot clearly indicates the divergence for $\gamma \to 0$, and as $\gamma$ becomes nonzero we see a decrease in $a_{opt}$ as expected—the processess should not be too slow since then the environmental fluctuations hinder perfect transfer. In this regime, we find $a_{opt}\tau \sim {\gamma }^{-1},$ which is expected as it says that the STIRAP time-scale should be of the same order as the relaxation time-scale. However, for larger $\gamma$, beyond 0.4, the optimal time-scale starts to grow again! One roughly finds

$$a_{opt}\tau \approx c_0+\frac{1}{\gamma }+c_1\gamma ,$$

(18)

for some constants $c_0$ and $c_1$ (that depend on the remaining system parameters). The linear $\gamma$-dependence setting in for larger decay rates is not intuitive. This suggests a slower STIRAP process when the coupling to the environment is getting larger. Of course, the mechanism behind this must be different from the one explained above leading to a ${\gamma }^{-1}$-dependence. If we study the actual form of the Liouvillian matrix (A7), we see that a non-zero $\gamma$ has the effect of decreasing all the Bloch vector components $r_i$ apart from $r_3$ and $r_8$ which represent the diagonal elements of the density operator. Thus, as expected, the dephasing diminish the of diagonal coherence terms. If this reduction is definite, the dynamics enters into the Zeno regime—the absence of coherences block the population transfer. It seems that the explanation to this unintuitive result is a Zeno effect. For large $\gamma$, it might be favorable to consider slow processes in order to prohibit further population transfer, which would lower the final population of the desired state. We should point out, however, that this regime occurs for decay rates that are on the same order as the couplings $g_0$, and one may question the justification of describing the system by a Markovian Lindblad master Equation (1). Experimentally, it might be possible to simulate this regime and remain Markovian by utilizing so called engineered dissipation [34], which means that the Lindblad jump operators are monitored by the experimenter by cleverly couple stable states to dissipative ones.

Case (b) is conceptually different from case (a). The STIRAP adiabatic state $|{\phi }_0\left(t\right)\rangle$ is a zero eigenvalue eigenstate of the Lindblad jump operators ${\widehat{L}}_{1,2}$ of Equation (12), meaning that it is an instantaneous steady state and moreover transparent to the environment. Perfect population transfer is thereby expected in the deep adiabatic regime. This is also found numerically, as demonstrated in Figure 8. Contrary to the case of dephasing, here the population of the target state is always a monotonously increasing function of a regardless of the value of $\kappa$, or, put in other words, a is a proper parameter to characterise efficiency of the process.

In principle, adiabaticity and efficiency need not be equivalent for open systems. In fact, adiabaticity can be introduced in various ways for open quantum systems. This issue has been raised in a series of papers [23,24,35,36]. One definition is to define it by saying that the population transfer between different adiabatic states of $\mathbf{M}$ vanishes (A4). At an EP, adiabaticity is then doomed to break down. Following this, the adiabatic criterium (9) can be modified to open systems in a rather direct way [23]. As is clear, the openness may greatly affect the criteria and typically the environment induces additional excitations [37,38]. In this respect, something interesting occurs for small a when the STIRAP is not fully adiabatic; the spontaneous emission increases the efficiency of the population transfer. This derives from the presence of a Liouvillian gap $\tilde{\Delta }$, which implies that the adiabatic state $|{\phi }_0\left(t\right)\rangle$ is partly protected from non-adiabatic excitations. One way of seeing it is that a non-zero $\kappa$ helps the process to become more adiabatic [35]. Thus, contrary to the standard situation where the fluctuations from the environment induces excitations, in the non-adiabatic regime, the environment prevents the system from taking it out from the instantaneous Hamiltonian adiabatic state. The phenomenon is related, but still different from the Zeno-effect discussed in Ref. [11]. There, in the overdamped regime, the environment prohibits the system from leaving the initial state $|1\rangle$. In the present situation, we are far from overdamped (except the dashed green line of Figure 8) and the environment induced relaxation takes you back onto the adiabatic state and not onto the initial state. Remember that the adiabatic state is also an instantaneous steady state and as soon as the Hamiltonian drives you out of this state, the environment pushes you back towards it. It should be clear that the same mechanism does not work for case (a) since the desired adiabatic state is not an instantaneous steady state of the full model. There, we found instead a different manifestation of the Zeno physics, namely a freezing of population transfer between the different diabatic states.

As a final remark, the jump operators (12) describe the decay into an incoherent mixture of $|1\rangle$ and $|3\rangle$. One may ask whether a “coherent” decay [39], like, for example, as represented by the Lindblad jump operator $\widehat{L}=|1\rangle \langle 2|+|3\rangle \langle 2|$, would affect the result of Figure 8. We tried for a couple of different such jump operators and found that the coherence does not play a quantitative role; any such decay protects the adiabatic state during fast processes.

4. Conclusions

Starting from the Liouvillian, we have analysed the open STIRAP problem. The situations of either dephasing or dissipation were considered separately and in both cases we found surprising results. An important difference between the cases is that for dissipation the Lindblad jump operator is Hermitian and one consequence is that, for the generalized Bloch Equation (14), one inevitably has a non-vanishing pump term $\mathbf{b}$ [32]. This, in return, implies that the steady state is non-trivial, in contrast to the case of dephasing where the instantaneous steady state is the maximally mixed one. In particular, the desired adiabatic state is also the instantaneous steady state for the case of dissipation, and the slower the process becomes, the more efficient population transfer to the target state. This is in agreement with the definition of adiabaticity for open quantum systems in terms of steady states [35], which states that, if the system is initially in a steady state, and the process is infinitely slow, then the system remains in the same instantaneous steady state for all times. For dissipation, the interesting new result is found for processes that are not perfectly adiabatic, i.e., during the evolution, non-adiabatic excitations take the system out from the instantaneous steady state. Here, it is found that the environment can “protect” the steady state in the sense that excitations out of it are counteracted by forcing the state to relax back to the instantaneous steady state. This is a sort of Zeno effect, and the system is locked to be in the adiabatic state.

In the situation of dephasing, we found that the evolution can be separated into two regimes depending on the strength of the decay rate $\gamma$. For weak $\gamma$, the environment indeed causes fluctuations that decrease the transfer efficiency, and one should therefore not make the process too slow. Instead, one should choose it as a trade-off between minimizing both the Hamiltonian non-adiabatic excitations and the environmental induced excitations. This result has been mentioned in the past [9,32,33,40], but here we suggest that the optimal time-scale for the process is of the same order as the inverse decay rate. What was surprising, however, was that, for larger decay rates $\gamma$, the numerics showed that a slower process might actually be advantageous. This is in conflict with the general knowledge that the environment will induce additional unwanted excitations. We argued that this is again a result of a Zeno-type effect. If the process is slow, all coherence terms of the state are quenched and the population transfer between the diabatic states is considerably slowed down. Of course, the Zeno freezing should not be perfect early on, since then there will be no population transfer at all. Instead, the “freezing” should occur at the time when the population in the target state is the largest.

To the best of our knowledge, both of these Zeno effects are new. In particular, they are different from the Zeno effect analyzed in Ref. [11]. In that work, the dynamics was completely frozen by the strong coupling to the environment, while, in our work, for dissipation, the dynamics are only partly frozen; the system is still evolving, but it follows an adiabatic state. Thus, the Zeno freezing we consider is not perfect in the sense that evolution still occurs.

An experimental system that could be suitable for exploring our theoretical predictions is trapped ions [41], where STIRAP is a standard tool for state preparation (see, for example, [42]). In these systems, both the initialisation and detection is almost perfect. By engineering the dissipation/decoherence, which has also been demonstrated experimentally [34], one should be able to construct the desired jump operators and the corresponding decay rates $\gamma$ and $\kappa$.