1. Introduction

The STIRAP process allows efficient population transfer in three level systems using two delayed laser pulses [1]. We consider the usual three level Λ-system {∣1⟩, ∣2⟩, ∣3⟩}, of respective energies E ₁ < E ₂ < E ₃, with no coupling between ∣1⟩ and ∣3⟩. The population is initially in level ∣1⟩. Units are chosen such that ħ = 1.

The STIRAP process consists in applying the Stokes laser pulse (approximately tuned to the Bohr frequency E ₃ - E ₂) before the pump laser pulse (tuned to E ₂ - E ₁). (The opposite sequence does not lead to complete transfer). We consider here for simplicity lasers exactly tuned to the one-photon resonances. At the initial and final times (when the fields are off), the dressed states (or Floquet states) are in resonance, and hence degenerate. For any system, the key of the transfer for this process is (i) the initial and final liftings of degeneracy which give rise to a transfer state connecting level ∣1⟩ to ∣3⟩, (ii) the adiabatic following of the dynamics on the transfer state [2, 3, 4, 5, 6].

The usual rotating wave approximation (RWA) allows to treat the initial and final resonances as the lowest order of the stationary perturbative theory on Floquet states: it determines the lifting of the degeneracy. It gives the first order terms for the dressed eigenvalues and the zeroth order terms for the eigenvectors. (The dimensionless parameters of the expansion correspond to the ratios between Rabi and Bohr frequencies.)

If we consider ideal adiabatic evolution, the corrections of the dynamics during the process are given by the counter-rotating terms, which are often neglected when considering the STIRAP process. We can study the corrections perturbatively as long as they do not induce new resonances between dressed states. These new resonances, which would appear as avoided crossing at non-zero fields, are called nonlinear resonances. We are in particular interested in the case when one of the peak Rabi frequencies approaches the difference of the two frequencies.

We restrict ourselves to this case of absence of nonlinear resonances and study a systematic perturbative development to improve the quantitative description of the dynamics during the process. We also study the limitations of this perturbative development when we approach a non-linear resonance.

2. The full Hamiltonian

The free three-level system is decribed by a Hamiltonian H ₀ on the Hilbert space 𝛨 = ℂ³ spanned by the vector set {∣1⟩, ∣2⟩, ∣3⟩}. It is driven by the two smooth pulsed-shaped monochromatic fields, with the dipole moment μ,

where the time-dependent field envelopes, carrier frequencies and initial phases of the fields are respectively denoted α̱ = (α_p ,α_s ), ω̱ = (ω_p ,ω_s ) and θ̱ = (θ_p ,θ_s ). For each fixed value of the fields, we can solve the time-dependent Schrödinger equation by the multi-mode Floquet theory [7, 3], which includes photon exchanges between matter and light [8]. This gives rise to the quasi-energy operator

It is defined in the enlarged space 𝛫 = 𝛨 ⊗ ℒ₂(dθ_p /2π) ⊗ ℒ₂(dθ_s /2π) where each ℒ₂(dθ_i /2π) is a space of square integrable functions of an angle θ_i , corresponding to a monochromatic photon field.

The eigenelements can be indexed with two indices: one, denoted n, refers to levels of the (dressed) molecule, and another one, denoted ḵ = (k_p , k_s ), for the relative photon number in each mode. The eigenvalues, which are two-mode periodic (quasiperiodic), are denoted λ_{n,k _} = λ_n,0 + ḵ∙ω and the eigenvectors ∣n,ḵ⟩.

Since the envelopes of the pulses vary slowly, we expand the solution of the time-dependent Schrödinger equation in the Floquet basis and apply adiabatic principles. If we consider as a first approximation “exact” adiabatic following of the transfer state, the time evolution can be written in terms of the eigenelements of K ^α_. In the following, we develop a systematic method to determine perturbatively the eigenelements of K ^α_. We consider here for simplicity that the field peak amplitudes are both α _max and equal couplings μ ₁₂ = μ ₂₃ (μ ₁₃ = 0).

3. The perturbative analysis

3.1 Preparing the Hamiltonian: The Rotating Wave Transformation

We start from the full Floquet Hamiltonian (2). It can be expressed as a 3 by 3 matrix (in the basis of H ₀), whose elements are θ-dependent. To calculate the exact eigenelements of K, we have to diagonalize the full Hamiltonian. That can be done numerically in a truncated Fourier decomposition for each frequency (this comes down to a discretization of the variables θ̱). The idea is to extract from the full Hamiltonian the dominant θ̱-independent terms in a perturbative series.

Because of the initial and final degeneracies, perturbative series cannot converge without a preliminary treatment of K. This treatment is the usual Rotating Wave Transformation (RWT) represented by the diagonal matrix:

It is denoted RWT as oposed to RWA because the counter-rotating terms are not discarded. We obtain (setting E ₂ = 0)

with

The usual RWA consists in neglecting the θ̱-dependent operator V ₁, i. e. the counter-rotating terms. We remark that the RWA is equivalent to the application (in one Floquet block) of quasi-degenerate stationary perturbation theory on the Floquet Hamiltonian to lowest order, i.e. just to take the good linear combinations in the degenerate subspace. The first term of Eq. (5) contains the counter-rotating terms of the pump laser on the 1–2 transition and of the Stokes laser on the 2–3 transition. The other terms correspond to the interactions of the pump laser on the 2–3 transition and of the Stokes laser on the 1–2 transition.

We next have to consider the diagonalization of the θ̱-independent part of the Hamiltonian (4)

where K ⁽⁰⁾ is the diagonalized usual STIRAP Hamiltonian

with the eigenvalues (including all the Brillouin zones) ${\mathrm{\lambda }}_{n,k}^{\left(0\right)}$ _{_} = ḵ ∙ ω̱ + ${\mathrm{\lambda }}_n^{\left(0\right)}$ , for n = {1, 2, 3},

The orthogonal matrix T ₀ contains the normalized eigenvectors of H ⁽⁰⁾ as column vectors.

We have thus written the transformed operator as K̃ = K ⁽⁰⁾ + εV ^(l) with $K^{\left(0\right)}=i\underset{̅}{\omega }·\frac{\partial }{\partial \underset{̅}{\theta }}+D^{\left(0\right)}$ , D ⁽⁰⁾ being diagonal and εV ^(l) = $T_0^{-1}$ V ₁ T ₀. We have introduced the formal parameter ε in order to treat eV ^(l) perturbatively.

3.2 The perturbative algorithm

We start with a quasienergy Hamiltonian K written (exactly) as

where ε is a small parameter. D ⁽⁰⁾ is diagonal and independent of θ̱.

We construct a unitary transformation exp(εW), with W† = -W antihermitian, such that

where D ^(l) is a diagonal part, of order e and independent of θ̱, and V ⁽²⁾ is a remaining correction of order ε ² (or higher). The unitary transformation reduces the size of the perturbation from order e to order e2. This method is known under different names, like “contact transformation”, KAM transformation, or van Vleck method [9, 10, 11]. Iterating this procedure is an alternative to expansions in power series (see e.g. [12]) which yields improved convergence [11]. In the present context, we will only do one step, which yields eigenvalues that contain all the corrections up to order ε ² and eigenvectors up to order ε. Maybe more importantly, this method allows one to distinguish in a systematic way the dominant contributions of the perturbation.

Inserting the unitary transformation in (10), expanding the exponential and identifying the terms of order ε, we obtain the equations that determine the unknown W and D ⁽¹⁾:

Expressing these equations in terms of the matrix elements with respect to the basis {∣m⟩} of eigenvectors of K ⁽⁰⁾ (we use a unique integer index m for simplicity), the solution of (11) can be written as

where we have denoted the eigenvalues of K ⁽⁰⁾ as ${\mathrm{\lambda }}_m^{\left(0\right)}$. The choice of W is not unique: one could add to it in (11) an arbitrary operator A that commutes with K ⁽⁰⁾. We choose A = 0.

In the present case, we have D ^(l) = 0, since we have already absorded the diagonal part into K ⁽⁰⁾.

The first three terms of the remaining correction of order ε can be written as:

4. The first corrections to the usual STIRAP

We start with the (full) prepared Hamiltonian (6). We apply once the previous scheme to first detect the dominant corrections in the usual regime of small Rabi frequencies:

We next treat the corrections pertubatively.

4.1 Dominant corrections

K ⁽⁰⁾ is defined by Eq. (7) and we have εV ₍₁₎ = $T_0^{-1}$ V ₁ T ₀.

By construction, we have the following two-mode Fourier developments

with the set ḵ = {(-2, 0); (2, 0); (0, -2); (0, 2); (-1, -1); (1,1); (-1,1); (1, - 1)} and, from the definition (12) of W,

for n= 1, 2, 3 and the eigenvalues ${\mathrm{\lambda }}_n^{\left(0\right)}$ defined in (8).

Taking into account the hypothesis (14), it appears clearly that the denominators appearing in W carry the dominant contribution for the set ḵ = {-ḵ̂; k̂} = {(-1,1); (1, -1)}. More precisely, these denominators become small when

Hence the second order gives the dominant contribution for the part of V ₁ corresponding to the modes {-ḵ̂;ḵ̂, i.e. for the last term of (5) [13].

4.2 Treatment of the corrections without nonlinear resonances

Keeping the dominant modes {-ḵ̂;ḵ̂} we obtain for the second order correction (the first commutator of (13)):

The second order corrections of the eigenvalues are given by the diagonal part of V ⁽²⁾:

with

and

The first order eigenvectors (θ̱-dependent) of $R_0^{-1}$ KR ₀ are given by:

This scheme is correct if the left-hand side of (17) does not approach too closely to ∣ω_p - ω_s ∣, otherwise the corresponding denominators become very small (and even zero) and induces the divergence of the perturbative scheme: this produces nonlinear resonances, that have to be tretated specifically with a second local RWT.

4.3 Population transfer in the adiabatic regime

It has been shown that at the first order the middle eigenvalue ${\mathrm{\lambda }}_2^{\left(0\right)}$ is always connected to level 1 at the beginning and to level 3 at the end of the process [3]. The second order eigenvalue (20) also connects 1 to 3, in the regime of absence of nonlinear resonances. In the adiabatic regime, this eigenvalue characterizes the transfer state.

Fig. 1a displays, for δ = 2 and α _max = 1 the second order eigenvalue curves (19) and (20), in comparison with the true quasienergies (obtained numerically): They are in quite good agreement. On Fig. 1b, the differences are plotted. We have also plotted the differences taking into account the diagonal part of the fourth order of V ⁽²⁾ (18). The accuracy is improved.

 

Figure 1. For δ = 2 and squared trig function pulse (of length 1 and delay 0.33): a) Exact (full lines) and second order (dashed lines) eigenvalue curves; b) Differences between the exact eigevalues and: the fourth order ones (full lines), the second order ones (dashed lines), and the ones from adiabatic elimination (dotted lines).

Download Full Size | PDF

We remark that a full description of the dynamics requires, besides the corrections due to counterrotating terms discussed here, corrections due to deviations from the adiabatic limit (nonadiabatic corrections). These deviations can be analyzed in terms of superadiabatic expansions [14, 15, 16, 17]. We notice that the corrections we obtain from the counterrotating terms in the present regime are larger than the nonadiabatic corrections obtained in Ref. [16]. For the present example, the corrections due to counterrotating terms do not affect the connectivity, but they can modify population of level 2 during the process.

5. Comparison with adiabatic elimination of dressed states

We compare the previous eigenvalues with the ones obtained from the Hamiltonian simplified with adiabatic elimination under the hypothesis $\sqrt{{\alpha }_p^2+{\alpha }_s^2}\ll \mid \delta \mid$:

More precisely, we obtain it following Ref. [18] for adiabatic elimination that is applied on the quasienergy operator projected in the relevant basis

The result is the STIRAP process with time-dependent Stark shifts (on-diagonal elements).

On Fig. 1b, we show the comparison of the eigenvalues obtained with different approximations.

6. Conclusion

In summary, we have discussed a systematic method to do perturbation analysis in the Floquet representation, based on an iterative scheme. We have calculated an explicit formula for corrections to second order of the eigenvalues. The comparison with the exact eigenvalues (computed numerically) shows a good agreement, provided that the peak intensities are sufficiently small to avoid nonlinear resonances. The results allow one to conclude that in this regime, the complete transfer of population is still possible. However, the transfer state contains a component on level 2 during the process. This may cause a partial loss of population, if level 2 is lossy.