A simple, analytic solution of the semiclassical Rabi model in the red-detuned regime

Physics Letters A 383 (2019) 1997–2003

Contents lists available at ScienceDirect

Physics Letters A www.elsevier.com/locate/pla

A simple, analytic solution of the semiclassical Rabi model in the red-detuned regime L.O. Castaños Instituto Tecnológico y de Estudios Superiores de Monterrey, Campus Santa Fe, Avenida de los Poetas 100, Santa Fe, La Loma, 01389 Ciudad de México, Mexico

a r t i c l e

i n f o

Article history: Received 14 February 2019 Received in revised form 15 March 2019 Accepted 27 March 2019 Available online 2 April 2019 Communicated by M.G.A. Paris

a b s t r a c t We consider the semiclassical Rabi model in the large, red-detuned regime. Using the method of multiple-scales we obtain a simple, analytic, and approximate solution that describes the evolution of the system accurately for long times and for arbitrary values of the qubit-field coupling. It is used to characterize the probability to find the qubit in the excited state and the trajectory of the associated Bloch vector. Finally, we present physical situations where the results can be applied. © 2019 Elsevier B.V. All rights reserved.

Keywords: Semiclassical Rabi model Dispersive deep strong-coupling Two-level system Semiclassical radiation theory

1. Introduction The semiclassical Rabi model (SRM) and its extensions [1–6] constitute a fundamental model used to describe the interaction of radiation with matter [7]. The SRM considers a qubit (i.e., a two level system) with transition frequency ω0 interacting with a single-mode, classical field with frequency ω1 . When ω1  ω0 and the qubit-field coupling is much smaller than ω0 , the model is simplified using the rotating-wave-approximation (RWA) [8], which consists in neglecting terms whose effect averages to zero because they oscillate rapidly. In this simpler form, the model has been used extensively to describe accurately physical systems such as those found in the area of cavity quantum electrodynamics (QED) [9,10]. Nevertheless, there are other physical systems where the RWA is not applicable because either the qubit-field coupling is comparable or larger than ω0 or |ω1 − ω0 | can be large [1,11–21]. This occurs, for example, in the area of circuit QED [1,11–13] and, in some cases, in atomic and molecular physics [14–17] and experiments with cold atoms, see [18–20] and references therein. As it so happens, the SRM is difficult to solve exactly, so several approximate treatments have been developed [15–17,22–25]. In particular, one treatment [17] considers the resonant and blue-detuned cases and is based on the multiple-scales method, which has the advantage of providing simple, analytic, and approximate solutions that describe the evolution of the system accurately for arbitrary qubitfield couplings and that can be interpreted and manipulated easily.

E-mail addresses: [email protected], [email protected]. https://doi.org/10.1016/j.physleta.2019.03.039 0375-9601/© 2019 Elsevier B.V. All rights reserved.

The purpose of this article is to extend the treatment presented in [17] to include the red-detuned regime, that is, the case where the transition frequency of the qubit is larger than the frequency of the field. The article is organized as follows. Section 2 introduces the model under consideration, while Section 3 presents the multiplescales solution, establishes its accuracy, and considers the probability to find the qubit in the excited state and the trajectory of the associated Bloch vector. In addition, it is discussed to which physical systems the results are applicable. Finally, a summary and the conclusions are given in Section 4. 2. The model Consider a qubit (that is, a two level system) interacting with a classical, single-mode field. An orthonormal basis for the state space H A of the system is {|2, |1}, where |1 is the ground state of the qubit and |2 is its excited state. The Hamiltonian of the system is

H (t ) =

h¯ ω0 2

σˆ z − h¯ 0 cos(ω1 t )σˆ x .

(1)

ω0 > 0 is the angular transition frequency of the qubit, ω1 > 0 is the angular frequency of the field, and 0 > 0 is an Here

angular frequency, called the Rabi frequency, that measures the strength of the qubit-field coupling. Also, σˆ x , σˆ y , and σˆ z are the Pauli operators given by

1998

L.O. Castaños / Physics Letters A 383 (2019) 1997–2003

d

σˆ z = |22| − |11|, σˆ + = |21|, σˆ − = |12|, σˆ x = σˆ − + σˆ + ,

σˆ y = i (σˆ − − σˆ + ).

(2)

In all that follows ρ (t ) denotes the density operator of the system. Its evolution is governed by von Neumann’s equation

ih¯

d dt

ρ (t ) = [H (t ), ρ (t )] .

(3)

Instead of working directly with the density operator, it is more convenient to introduce the Bloch vector r(t ) given by

r(t ) = (α1 (t ), α2 (t ), α3 (t )) T ,

(4)

with T the transpose and

α2 (t ) = 2Im [1|ρ (t )|2] , α3 (t ) = 2|ρ (t )|2 − 1|ρ (t )|1.

(5)

Recall that |r(t )| ≤ 1 for all times. Here and in the following |v| is the Euclidean or 2-norm of a vector v. We note that one can reconstruct the density operator from the Bloch vector, since it is straightforward to show that

1

1

2

2



α1 (t )σˆ x + α2 (t )σˆ y + α3 (t )σˆ z ,

(6)

with I the identity operator in H A . Using von Neumann’s equation in (3), one can obtain an evolution equation for the Bloch vector. It is given by

d dt





r(t ) = −20 cos(ω1 t )ˆx + ω0 zˆ × r(t ),

(7)

where xˆ , yˆ , and zˆ are unit vectors in the positive direction of the x-, y-, and z-axes. Observe that this is the equation of a point r(t ) that is rotating  around the vector [−20 cos(ω1 t )ˆx + ω0 zˆ ] with angular speed

2 0

ω

+ 420 cos2 (

ω1 t ) [26].

The equation in (7) is quite difficult to solve. In the next section we present a simple, approximate solution using the multiplescales method. 3. The multiple-scales solution The resonant and blue-detuned ω1 > ω0 cases were presented in [17]. In this article we treat the large, red-detuned case, that is, in the rest of the article we assume that

ω0  ω1 .

(8)

Given the relationship in (8), it is convenient to measure time in units of 1/ω0 . Hence, we introduce the non dimensional quantities





(10)

Again, this is the equation of a point that is rotating around

˜ 0 cos(ω˜ 1 τ )ˆx + zˆ ] with the angular speed [−

˜ 2 cos2 (ω˜ 1 τ ). 1+ 0

˜ 0 , 1/ω˜ 1 , and 1. One can recognize three time scales in (10): 1/ ˜ 1 . Then, one has two competing From (8) one knows that 1  1/ω ˜ 0 cos(ω˜ 1 τ )ˆx + zˆ ] processes: a fast rotation around the vector [− and a slow variation with respect to time of this vector. It follows that the multiple-scales method [27] can be used to obtain an accurate approximation that holds for long times. From (8) one knows that a slow time-scale is

˜ 1τ . t2 = ω

(11)

Using the multiple-scales method [27] it follows that an appropriate fast time scale is given by (see the Appendix for the details of its deduction)

α1 (t ) = 2Re [1|ρ (t )|2] ,

ρ (t ) = I +

dτ



˜ 0 cos(ω˜ 1 τ )ˆx + zˆ × r0 (τ ). r0 (τ ) = −



τ , τ = ω0t , r0 (τ ) = r ω0 ˜ 0 = 20 , ω˜ 1 = ω1 .  ω0 ω0

t1 =



˜2 1+ 0

τ 0



  ˜2     0  ˜ 1τ  . dτ 1− sin2 ω 2 ˜ 1 + 0

Notice that the integral in (12) is an elliptic integral of the second kind [26]. To get a better idea of the fast time scale t 1 , one can consider ˜ 0  1, then one can make a linear approxtwo limiting cases. If  ˜ 0 to get imation in 

t1 = τ .

(13)

˜ 0  1, so 1 is the fast time-scale. In this case  ˜ 0  1, then one can make a linear apOn the other hand, if  ˜ proximation in 1/0 to get

    ˜ 0  2ω˜ 1 τ     ˜ 1 τ π   2ω  + sin ω˜ 1 τ − sin , t1 = ω˜ 1 π π 2 

t1 

2

π

˜ 0τ . 

(15)

˜ 0. ˜ 0  1, so the fastest time scale is 1/ In this case  Using the multiple-scales method [27] and the time-scales defined in (11) and (12) one can obtain a first term approximation r00 (τ ) for r0 (τ ) (see the Appendix for the details of the calculation): r00 (τ ) = U0 (t 1 , t 2 )U1 [θ(t 2 )] r(0).

(16)

Here, U0 (t 1 , t 2 ) is a 3 × 3 rotation matrix (a real, orthogonal matrix with determinant equal to 1) given by

⎛

−U 21

U 11 ⎝ U 21 U0 (t 1 , t 2 ) = y (t 2 ) U 31

Notice that τ is a non dimensional time and that r0 (τ ) is the Bloch ˜ 0 and vector as a function of the non dimensional time τ . Also,  ω˜ 1 express the size of the qubit-field coupling 0 and field frequency ω1 in terms of ω0 . Using (7) and (9) one can obtain an evolution equation for r0 (τ ):

(14)

with a the greatest integer less than or equal to a. To obtain a ˜ 1 τ /π  1 more transparent expression, one can consider times 2ω ˜ 1 τ /π  2ω˜ 1 τ /π in so that one can make the approximation 2ω (14) to get

1

(9)

(12)

U 22 U 32

⎞

U 31 −U 32 ⎠ , U 33

(17)

with

˜ 20 cos2 (t 2 ), U 33 = 1 +  ˜ 20 cos2 (t 2 )cos(t 1 ), y (t 2 ) = 1 +  ˜ 20 cos2 (t 2 ), U 22 = cos(t 1 ) y (t 2 ), U 11 = cos(t 1 ) +  

U 21 = sin(t 1 )

y (t 2 ),

˜ 0 cos(t 2 )sin2 U 31 = −2



˜ 0 cos(t 2 )U 21 , U 32 = − t1 2



.

(18)

L.O. Castaños / Physics Letters A 383 (2019) 1997–2003

1999

Also, U1 [θ(t 2 )] is a rotation matrix around the y-axis:

⎛

cos [θ(t 2 )] 0 U1 [θ(t 2 )] = ⎝ −sin [θ(t 2 )]

⎞

0 sin [θ(t 2 )] ⎠, 1 0 0 cos [θ(t 2 )]

(19)

with

cos [θ(t 2 )] =

˜ 2 cos(t 2 ) 1+ 0

√



˜2 y (t 2 ) 1 +  0

,

˜ 0 [1 − cos(t 2 )]   . ˜2 y (t 2 ) 1 + 

sin [θ(t 2 )] = √

(20)

0

˜ 0 cos(t 2 )ˆx + zˆ is an eigenvector One can show that v(t 2 ) = − of U0 (t 1 , t 2 ) with corresponding eigenvalue 1. Notice that v(t 2 ) is the vector appearing in the evolution equation of r0 (τ ) given in (10). Hence, from (16) we have, to good approximation, that the evolution of the Bloch vector is given by a rotation around the y-axis followed by a rotation around the vector v(t 2 ). 3.1. Accuracy of the multiple-scales solution We now determine the accuracy of the multiple-scales solution r00 (τ ) in (16). Equation (10) was solved numerically and the maximum relative error between the exact (numerical) solution r0 (τ ) and the multiple-scales solution r00 (τ ) was determined during the ˜ 1 for any initial condition inside or time interval 0 ≤ τ ≤ 5(2π )/ω on the Bloch sphere (the sphere centered at the origin with radius 1). In other words, we calculated the following quantity:

   |r (τ ) − r (τ )| 2π 0 00 , :0≤τ ≤5 |r0 (τ )| ω˜ 1  |r0 (0)| ≤ 1 .

˜ 0 , ω˜ 1 ) = max E (

(21)

˜ 0 , ω˜ 1 ) for 1 ≤  ˜ 0 ≤ 4 and 0.005 ≤ Fig. 1a shows contours of E ( ˜ 0 ≤ 1 and ω˜ 1 ≤ 0.02, while Fig. 1b does the same for 0.3 ≤  ˜ 1 ≤ 0.1. Observe that the multiple-scales solution is ac0.02 ≤ ω ˜ 0 ≥ 1 if curate for long times and for large qubit-field couplings  ω˜ 1 is sufficiently small. In particular, it can describe accurately the ˜1  1 dispersive deep strong-coupling regime [18] characterized by ω ˜ 0 /4 ≤ 1 or, using (9), ω1  ω0 and ω1 < 0 /2 ≤ ω0 ˜1 <  and ω (notice that the parameter g of [18] is 0 /2). Notice that the contours in Fig. 1 have approximately the form ˜ 0 = k with k > 0 a constant. This shape can be explained if ω˜ 1  one determines how close the multiple-scales solution is to satisfying the differential equation in (10). Substituting (16) in (10) one obtains d dτ





˜ 0 cos(ω˜ 1 τ )ˆx + zˆ × r00 (τ ) = F(t 1 , t 2 )r0 (0), (22) r00 (τ ) − −

where

˜ 0 sin(t 2 ) E(t 1 , t 2 ), F(t 1 , t 2 ) = −ω˜ 1  y (t 2 )3/2 ⎞ ⎛ ˜ 0 cos(t 2 )sin(t 1 ) E 13 E 11  E(t 1 , t 2 ) = ⎝ 0 0 0 ⎠, E 31 −sin(t 1 ) E 33 and

˜0 

E 11 =  [1 − cos(t 1 )cos(t 2 )] , ˜2 1+ 0

Fig. 1. The figures show the maximum relative error (21) that the multiple-scales ˜ 1 . Fig. 1a considers solution (16) can have during the time interval 0 ≤ τ ≤ 5(2π )/ω ˜ 0 ≤ 4 and 0.005 ≤ ω˜ 1 ≤ 0.02, while Fig. 1b illustrates the parameter the values 1 ≤  ˜ 0 ≤ 1 and 0.02 ≤ ω˜ 1 ≤ 0.1. range 0.3 ≤ 

E 31 =

˜ 2 cos(t 2 ) cos(t 1 ) +  0

E 13 = −



˜2 1+ 0

,

˜ 2 cos(t 1 )cos(t 2 ) 1+ 0

E 33 = 



˜0  ˜2 1+ 0

˜2 1+ 0

,

[cos(t 1 ) − cos(t 2 )] .

One can consider r00 (τ ) to be a good approximation to r0 (τ ) if the righthand side of (22) is very small. To determine when this happens one has to distinguish between two types of time intervals. If cos(t 2 )  ±1, then it is straightforward to show that the absolute value of each matrix element of F(t 1 , t 2 ) is bounded by ˜ 1 |sin(t 2 )|, which is going to be very small as long as ω˜ 1  1. 2ω The other time interval occurs when cos(t 2 )  0. In this case, it is straightforward to show that the absolute value of each matrix ˜ 0 , which is going to be ˜ 1 element of F(t 1 , t 2 ) is bounded by ω ˜ 0  1. Combining the two results one ˜ 1 very small as long as ω concludes that the multiple-scales solution r00 (τ ) is an accurate approximation if

˜ 0  1 ⇔ ω1  ω0 , 20 ω1  ω2 . ω˜ 1 , ω˜ 1  0 (23)

(24)

(25)

Notice that we used the definitions in (9) in the righthand side of (25). It is important to contrast the cases of red- and blue-detuning. For red-detuning (ω0 > ω1 ) it was established above that the multiple-scales solution is accurate if (25) holds. In particular, if ˜ 0 are small, then (25) is satisfied. For large qubit˜ 1 and  both ω ˜ 0 ≥ 1, one requires a very small field frequency field couplings 

2000

L.O. Castaños / Physics Letters A 383 (2019) 1997–2003

ω˜ 1 so (25) is satisfied. On the other hand, for blue-detuning (ω0 < ω1 ) it was found that, for fixed values of ω0 and ω1 , the multiplescales solution became more accurate for increasing qubit-field coupling 0 [17]. To illustrate why the multiple-scales method appears to have difficulties in describing the red-detuned regime, assume that ˜ 0 . Numerical calculations show that the components ω˜ 1  1   of the Bloch vector exhibit large, rapid changes in narrow regions ˜ 0 is not suf˜ 1 centered at the times where cos(t 2 ) = 0. When ω ficiently small, these regions are very narrow and the multiplescales solution (16) cannot describe the rapid changes. In this case, the numerical solution exhibits fast oscillations determined by t 1 , an envelop of these fast oscillations determined by t 2 (see the next section), and there appears to be an even slower time-scale that manifests itself as an envelop of the envelop determined by t 2 . In the blue-detuned case, it is better to first eliminate the free evolution of the qubit, that is, it is better to pass to the interaction picture (IP) defined by the unitary operator U I (t ) = exp(−i ω0 t σˆ z /2). The IP Bloch vector r I (t ) is determined by [17]

d dt

r I (t ) = w(t ) × r I (t ),

(26)

with





w(t ) = −20 cos(ω1 t ) cos(ω0 t )ˆx − sin(ω0 t )ˆy .

(27)

To get an idea if the multiple-scales method is going to work, one can take the slow time-scale to be constant in (27) and see if one is going to be able to solve (26). For blue-detuning, the slow timescale is ω0 t [17]. If one considers it to be constant in (27), then (26) is exactly solvable [17]. For red-detuning, the slow time-scale is ω1 t. If one considers it to be constant in (27), then (26) is still difficult to solve. Hence, the IP is useful only for blue-detuning.

˜ 0 = 1, ω˜ 1 = 0.03, and the initial condition ˜ 1 for  the time interval 0 ≤ τ ≤ 2π /ω r(0) = (0, 0, −1) T . The exact Bloch vector was obtained by solving numerically (10).

3.2. A special case In this section we assume that the qubit is initially in the ground state and that we are in a parameter regime where the multiple-scales solution is accurate (say, the maximum relative error is ≤ 15%). Then, the initial Bloch vector is r0 (0) = r(0) = (0, 0, −1)T and from (16) one finds that

⎛

r00 (τ ) =

Fig. 2. (Color online) The figures show the components of the exact Bloch vector r0 (τ ) = (α1 (τ ), α2 (τ ), α3 (τ )) T (red-solid lines) and, except for Fig. 2b, the envelops r−E (τ ) (blue-dashed lines) and r+E (τ ) (black-dotted line) in (29) during

⎞

cos ˜0 √(t 1 ) − cos(t 2 ) − ⎠. ⎝ y (t 2 )sin(t 1 )  √ ˜2  ˜ 0 cos(t 2 )cos(t 1 ) ˜ −1 +  y (t 2 ) 1 +  0 0

(28)

The envelops of the components of r00 (τ ) are obtained by eliminating the terms associated with the fast time scale t 1 :

⎞ ⎛ cos ˜0 √ (t 2 )  ⎠. ⎝ ± y (t 2 )  r±E (τ ) = √ ˜ 2 − ˜ 0 cos(t 2 ) ˜ −1 ±  y (t 2 ) 1 +  0 0

(29)

Fig. 2 shows the components of the exact solution r0 (τ ) = (α1 (τ ), α2 (τ ), α3 (τ ))T and the components of the envelops (29), ˜ 1 for except for Fig. 2b, during the time interval 0 ≤ τ ≤ 2π /ω ˜ 0 = 1, ω˜ 1 = 0.03, and the initial condition r(0) = (0, 0, −1) T . The  exact solution was obtained by solving (10) numerically and the relative error between the exact solution and the multiple-scales ˜ 1 and solution (28) is ≤ 4% during the time interval 0 ≤ τ ≤ 2π /ω < 13% during the extended time interval 0 ≤ τ ≤ 5(2π /ω˜ 1 ). Observe that the envelop of the first component α1 of r0 (τ ) is not actually an envelop but the midpoint of the oscillation of α1 , see Fig. 2a. Also, the envelops of the second component α2 were not included because they are simply horizontal lines, see (29), and there is a slight oscillation in α2 .

Using the definition of α3 in (5) and the fact that the density operator ρ (t ) has trace equal to 1, it follows from (28) that the probability P 2 (τ ) to find the qubit in the excited state is

P 2 (τ ) =

=

1 + α3 (τ ) 2 1 2

−

,

˜ 2 cos(t 2 )cos(t 1 ) 1+ 0 2

√



˜2 y (t 2 ) 1 +  0

.

(30)

This probability has fast oscillations determined by the fast time scale t 1 that are contained within envelops P ±2 (τ ) that depend on the slow time scale t 2 :

P ±2 (τ ) =

˜ 2 cos(t 2 ) 1± 0  − √ . 2 ˜2 2 y (t 2 ) 1 +  0 1

(31)

Since, P 2 (τ ) = 0 when α3 (τ ) = −1, see (30), from Fig. 2c one can observe that the envelops P ±2 (τ ) spend more time near their absolute minimum than anywhere else. In addition, P ±2 (τ ) have a form similar to the probability to find the qubit in the excited state in the large blue-detuned regime ω0  ω1 [16]. ˜ 1 that Notice that P ±2 (τ ) are periodic functions of period 2π /ω can be used to analyze P 2 (τ ). In particular, one has the following properties:

• P +2 (τ ) = 0 ⇔ ω˜ 1 τ = 2π n with n ≥ 0 an integer. Also, ˜ 1 τ = π + 2π n with n ≥ 0 an integer. In the P −2 (τ ) = 0 ⇔ ω neighborhood of these values the graph of the corresponding envelop is almost flat, see Fig. 2c with α3 = −1.

L.O. Castaños / Physics Letters A 383 (2019) 1997–2003

2001

the semiclassical Rabi model and the aforementioned conditions can be satisfied [1]. The approximate solution is also useful to describe the interaction of a single-mode, classical electromagnetic field with an ion [14], a molecule [15], an atom [17], and, in particular, a cold atom [18–20]. In the case of alkali atoms, transitions between Zeeman sublevels of the ground state configuration are especially relevant because spontaneous emission is negligible and the separation between the levels can be controlled using the Zeeman effect and a static magnetic field [17,19,20]. In this case, the dispersive deep strong-coupling regime [18] with ω1  ω0 and ω1 < 0 /2 ≤ ω0 can be considered [17,19,20]. 4. Conclusions Fig. 3. (Color online) The figure shows the trajectory (red-solid line) of the Bloch ˜ 1 with vector r0 (τ ) on the Bloch sphere during the time interval 0 ≤ τ ≤ 2π /ω ˜ 0 = 0.3, ω˜ 1 = 0.1, and the initial condition r0 (0) = (0, 0, −1) T . The south pole  appears as a black dot and r0 (τ ) was calculated numerically form (10).

• If cos(ω˜ 1 τ ) = −1, then P +2 (τ ) =

˜2  0 . ˜2 1+

(32)

0

˜ 0 so One can use (32) to determine the range of values of  that there are times τ for which P 2 (τ ) ≥ k with 0 < k < 1. One finds that



˜0≥ P 2 (τ ) ≥ k ⇔ 

k 1−k

(33)

.

In particular, there are times τ for which P 2 (τ ) ≥ 0.5 if and ˜ 0 ≥ 1 (that is, 0 ≥ ω0 /2) and there are times τ for only if  ˜ 0 ≥ 3 (that is, 0 ≥ 1.5ω0 ). which P 2 (τ ) ≥ 0.9 if and only if  These results are similar to the blue-detuned ω1 > ω0 case where one requires 0 ≥ πω1 /2 so that there are times τ such that P 2 (τ ) = 1 and one requires 0 ≥ πω1 /4 so that there are times τ such that P 2 (τ ) ≥ 0.5 [17]. The trajectory of the Bloch vector can be characterized easily ˜ 0  1. Makwhen the qubit-field coupling is small. Assume that  ˜ 0 in (28) one finds that ing a linear approximation in 

⎛

r00 (τ ) = ⎝

˜ 0 cos(t 2 )  0

−1

⎞

⎛

⎞

cos(t 1 ) ˜ 0 ⎝ sin(t 1 ) ⎠ . ⎠− 0

(34)

Observe that (34) is the parametrization of a circumference in the ˜ 0 whose center moves along the x-axis z = −1 plane of radius  ˜ 0 to  ˜ 0 . Moreover, the Bloch vector r00 (τ ) moves rapidly from − along the circumference (it depends on the fast time scale t 1 ), while the center of the circumference moves slowly (it depends on the slow time scale t 2 ). Fig. 3 shows that this is a very good approximation to the actual trajectory of the Bloch vector calculated numerically from (10). We note that, for the parameters used in Fig. 3, the relative error between the multiple-scales solution (28) and the exact solution is < 3.5%. 3.3. Applications to physical systems The results of this article are applicable whenever the physical system under consideration is described by the Hamiltonian given (1) and the conditions ω1  ω0 and 2ω1 0  ω02 hold. We note that the familiar rotating-wave-approximation (RWA) is not applicable because it requires a quasi-resonant situation ω1  ω0 [8]. In particular, the results are of use in the area of circuit quantum electrodynamics, since there are experiments described by

In this article we considered the semiclassical Rabi model in the large, red-detuned regime, that is, a qubit (i.e., a two level system) interacting with a single-mode, classical field in the case where the transition frequency ω0 of the qubit is larger than the frequency ω1 of the field. Using the method of multiple-scales we obtained a simple, analytic, and approximate solution that accurately describes the evolution of the system for long times if ω1  ω0 and 2ω1 0  ω02 with 0 the qubit-field coupling. In particular, it describes accurately the dispersive deep strong-coupling regime with ω1  ω0 and ω1 < 0 /2 ≤ ω0 . The analytic solution was used to determine the probability P 2 to find the qubit in the excited state and the trajectory of the associated Bloch vector. In particular, P 2 is described conveniently in terms of envelops and one requires strong qubit-field couplings 0 ≥ ω0 /2 so that one can have P 2 ≥ 0.5. On the other hand, the Bloch vector has as a trajectory a circumference whose center oscillates slowly along the x-axis when 20 < ω0 . In addition, it was discussed that the results of the article are applicable in the areas of circuit quantum electrodynamics, cold atoms, and atomic and molecular physics. It is interesting to contrast the multiple-scales solution in the red-detuned ω1 < ω0 regime with that in the blue-detuned ω1 > ω0 regime presented in [17]. In the blue-detuned regime it was found that, for a fixed detuning |ω1 − ω0 |, the approximate solution is increasingly accurate for larger couplings 0 . In the reddetuned regime, the opposite behavior was found, since the approximate solution is increasingly accurate for smaller couplings. It is important to mention that the fast and slow time-scales deduced in this article hold even when the multiple-scales solution is not accurate. Numerical calculations show that there appears an even slower time-scale involved in the evolution of the system. If this slower time-scale was found, then a simple, analytic solution could be obtained to describe the system in a larger parameter regime. 5. Appendix In this appendix we show how to derive the fast time-scale t 1 presented in (12) and the multiple-scales solution given in (16). Expanding the cross product in (10), one can write

d dτ

˜ 1 τ )r0 (τ ), r0 (τ ) = A(ω

where

⎛

(35)

⎞

0 −1 0 ˜ 0 cos(ω˜ 1 τ ) ⎠ . 0  A(ω˜ 1 τ ) = ⎝ 1 ˜ 0 cos(ω˜ 1 τ ) 0 − 0

(36)

˜ 1  1, so t 2 = ω˜ 1 τ given in (11) is the From (8) and (9) one has ω slow time-scale. Nevertheless, one does not know the fast timescale t 1 , so one takes t 1 = f (τ ),

(37)

2002

L.O. Castaños / Physics Letters A 383 (2019) 1997–2003

where f (τ ) is a function that satisfies the following properties: 1. f (τ ) ≥ 0, because we consider times t 1 , τ ≥ 0. 2. f  (τ ) > 0, because a time-scale is strictly increasing. ˜ 1 τ  f (τ ) if ω˜ 1 ↓ 0, because t 1 is the fast time-scale and t 2 3. ω is the slow time-scale. 4. f (0) = 0, because we demand that t 1 starts at zero just as t 2 and τ start at zero. 5. f (τ ) is continuously differentiable, because the calculations that follow require this property. To apply the multiple-scales method one has to define a new Bloch vector R(t 1 , t 2 ) by

R [t 1 (τ ), t 2 (τ )] = r0 (τ ).

(38)

Substituting (38) into (35) and using (37), one obtains

∂ Rk (t 1 , t 2 ) = −i Q(t 2 )DQ† (t 2 )Rk (t 1 , t 2 ) ∂ t1 ∂ Rk−1 1 −√ (t 1 , t 2 ), y (t 2 ) ∂ t 2 Rk (0, 0) = r0 (0)δk0 ,

(39)

Now, the fast time-scale in (37) is determined by the eigenvalues of A [t 2 (τ )]. It has to be chosen so that the time-dependence of the eigenvalues of A [t 2 (τ )] is eliminated when one divides by df /dτ . Diagonalizing A(t 2 ) one obtains †

A(t 2 ) = −i y (t 2 )Q(t 2 )DQ (t 2 ),

(40)

√1 1 ⎝ i y (t 2 ) Q(t 2 ) = √ 2 y (t 2 ) ˜ 0 cos(t 2 )  ⎞ ⎛

√ ˜ 0 cos(t 2 ) − 2

1 0 0 D = ⎝ 0 0 0 ⎠. 0 0 −1

√0

2

√1

where U0 (t 1 , t 2 ) = Notice that the initial condition R0 (0, 0) = r0 (0) cannot be applied yet. For k = 1 one obtains

R1 (t 1 , t 2 ) = U0 (t 1 , t 2 )R1 (0, t 2 ) 1 −√ Q(t 2 )e −i Dt1 E(t 1 , t 2 ), y (t 2 )

df dτ

(τ ) =



y [t 2 (τ )].

t1 E(t 1 , t 2 ) =

−i y (t 2 ) ⎠ , ˜ 0 cos(t 2 ) 

(42)

˜ 1τ ) = 1 − Integrating (42) from 0 to τ , using the identity cos2 (ω ˜ 1 τ ), and applying the condition f (0) = 0 in item 4 above, sin2 (ω one obtains the result in (12) and it is straightforward to show that it satisfies the properties in items 1-5 above. Omitting the τ dependence of t 1 and t 2 , it follows from (39)-(42) that

∂R ω˜ 1 ∂ R (t 1 , t 2 ) + √ (t 1 , t 2 ) ∂ t1 y (t 2 ) ∂ t 2



dt 1 e i Dt1 Q† (t 2 )

0

t1 =

Notice that Q(t 2 ) is a unitary matrix √ and that the eigenvalues of A(t 2 ) are the diagonal entries of −i y (t 2 )D . From (39)-(41), it follows that the time-dependence of the eigenvalues of A [t 2 (τ )] disappears if one chooses

(48)

where

⎞

(41)

(47)

Q(t 2 )e −i Dt1 Q† (t 2 ) is given in (17) and (18).

where y (t 2 ) is defined in (18) and

⎛

(46)

for k = 0, 1, 2, .... Here δk0 is the Kronecker delta and R−1 (t 1 , t 2 ) = 0. The differential equations in (46) can be solved using the exponential of a matrix [28]. For k = 0 one gets

R0 (t 1 , t 2 ) = U0 (t 1 , t 2 )R0 (0, t 2 ),

∂R df ∂R ˜1 (τ ) [t 1 (τ ), t 2 (τ )] + ω [t 1 (τ ), t 2 (τ )] dτ ∂ t1 ∂ t2 = A [t 2 (τ )] R[t 1 (τ ), t 2 (τ )].



˜ 1 = ω1 /ω0  1 is a perturbation parameter, see (8) Recall that ω and (9). We are interested in obtaining a one-term approximation to R, that is, we want to calculate R0 . Substituting (45) into (43) and (44) and equating equal powers ˜ 1 , one obtains the following initial value problems: of ω



dt 1 e i Dt1 Q† (t 2 )

0



+t 1

dQ dt 2

∂ R0  (t , t 2 ), ∂ t2 1 

(t 2 )e −i Dt1 Q† (t 2 )R0 (0, t 2 )

 ∂ R0 (t 2 )R0 (0, t 2 ) + Q (t 2 ) (0, t 2 ) . dt 2 ∂ t2

dQ†

†

(49)

Recall that the derivative of a matrix is obtained by calculating the derivative of all of its entries. From (48) and (49) it is straightforward to show that, if t 2 is fixed and t 1 is variable, then all the terms in the righthand side of (48) are bounded except those that are multiplied by t 1 in the second line of (49). These terms increase without bound as t 1 increases and destroy de order of the asymptotic expansion in (45). Hence, they are secular terms and must be eliminated. The dependence of R0 (0, t 2 ) on t 2 is used to eliminate these secular terms. One demands

dQ† dt 2

(t 2 )R0 (0, t 2 ) + Q† (t 2 )

∂ R0 (0, t 2 ) = 0. ∂ t2

(50)

(43)

Since Q(t 2 ) is a unitary matrix, one can solve for ∂ R0 /∂ t 2 to obtain

In addition, using (38) and that t 1 (0) = t 2 (0) = 0, it follows that R(t 1 , t 2 ) is subject to the initial condition

(51)

= −i Q(t 2 )DQ† (t 2 )R(t 1 , t 2 ).

R(0, 0) = r0 (0).

(44)

⎛ ⎞ ˜ 0 sin(t 2 ) 0 0 −1 ∂ R0  ⎝ 0 0 0 ⎠ R0 (0, t 2 ). (0, t 2 ) = − ∂ t2 y (t 2 ) 1 0 0

We now solve (43) and (44) to good approximation. Assume that R has an asymptotic expansion of the form

Given that the dependence on t 2 is factored out as a scalar multiplying a matrix in the righthand side of (51), this equation can be solved using the exponential of a matrix [28]:

˜ 1 R1 + ω˜ 12 R2 + .... R ∼ R0 + ω

R0 (0, t 2 ) = U1 [θ(t 2 )] r0 (0),

(45)

(52)

L.O. Castaños / Physics Letters A 383 (2019) 1997–2003

where U1 [θ(t 2 )] is defined in (19) and (20). Notice that we applied the initial condition R0 (0, 0) = r0 (0). From (47) and (52) one finally obtains that

R0 (t 1 , t 2 ) = U0 (t 1 , t 2 )U1 [θ(t 2 )] r0 (0).

(53)

Using the definition of R(t 1 , t 2 ) in (38), the asymptotic expansion in (45), and the expression for R0 (t 1 , t 2 ) in (53), one concludes that a one-term approximation to r0 (τ ) is given by (16), since

r0 (τ ) = R [t 1 (τ ), t 2 (τ )]  R0 [t 1 (τ ), t 2 (τ )] ,

= U0 [t 1 (τ ), t 2 (τ )] U1 {θ [t 2 (τ )]} r0 (0).

(54)

2003

[11] M. Devoret, S. Girvin, R. Schoelkopf, Circuit-qed: how strong can the coupling between a Josephson junction atom and a transmission line resonator be?, Ann. Phys. 16 (2007) 767, https://doi.org/10.1002/andp.200710261. [12] T. Niemczyk, et al., Circuit quantum electrodynamics in the ultrastrongcoupling regime, Nat. Phys. 6 (2010) 772, https://doi.org/10.1038/nphys1730. [13] P. Forn-Díaz, J. Lisenfeld, D. Marcos, J.J. García-Ripoll, E. Solano, C.J.P.M. Harmans, J.E. Mooij, Observation of the Bloch-Siegert shift in a qubit-oscillator system in the ultrastrong coupling regime, Phys. Rev. Lett. 105 (2010) 237001, https://doi.org/10.1103/PhysRevLett.105.237001. [14] J.S. Pedernales, I. Lizuain, S. Felicetti, G. Romero, L. Lamata, E. Solano, Quantum Rabi model with trapped ions, Sci. Rep. 5 (2015) 15472, https://doi.org/ 10.1038/srep15472. [15] S.H. Autler, C.H. Townes, Stark effect in rapidly varying fields, Phys. Rev. 100 (1955) 703, https://doi.org/10.1103/PhysRev.100.703. [16] K.K. Shakov, J.H. McGuire, Population control of 2s-2p transitions in hydrogen, Phys. Rev. A 67 (2003) 033405, https://doi.org/10.1103/PhysRevA.67.033405.

References [1] K. Dai, et al., Quantum simulation of the general semi-classical Rabi model in regimes of arbitrarily strong driving, Appl. Phys. Lett. 111 (2017) 242601, https://doi.org/10.1063/1.5006745. [2] R. Grimaudo, A.S.M. de Castro, H. Nakazato, A. Messina, Classes of exactly solvable generalized semi-classical Rabi systems, Ann. Phys. 530 (2018) 1800198, https://doi.org/10.1002/andp.201800198. [3] B. Cheng, Q.H. Wang, R. Joynt, Transfer matrix solution of a model of qubit decoherence due to telegraph noise, Phys. Rev. A 78 (2008) 022313, https:// doi.org/10.1103/PhysRevA.78.022313. [4] I. Goychuk, P. Hanggi, Quantum two-state dynamics driven by stationary nonMarkovian discrete noise: exact results, Chem. Phys. 324 (2006) 160, https:// doi.org/10.1016/j.chemphys.2005.11.026. [5] M.A.C. Rossi, M.G.A. Paris, Non-Markovian dynamics of single- and two-qubit systems interacting with Gaussian and non-Gaussian fluctuating transverse environments, J. Chem. Phys. 144 (2016) 024113, https://doi.org/10.1063/1. 4939733. [6] L. Cywinski, Dynamical-decoupling noise spectroscopy at an optimal working point of a qubit, Phys. Rev. A 90 (2014) 042307, https://doi.org/10.1103/ PhysRevA.90.042307. [7] D. Braak, et al., Semi-classical and quantum Rabi models: in celebration of 80 years, J. Phys. A, Math. Theor. 49 (2016) 300301, https://doi.org/10.1088/17518113/49/30/300301. [8] C. Cohen-Tannoudji, J. Dupont-Roc, G. Grynberg, Atom-Photon Iteractions: Basic Processes and Applications, Wiley-VCH, 1998. [9] H. Walther, B.T.H. Varcoe, B.G. Englert, T. Becker, Cavity quantum electrodynamics, Rep. Prog. Phys. 69 (2006) 1325, https://doi.org/10.1088/0034-4885/ 69/5/R02. [10] J.M. Raimond, M. Brune, S. Haroche, Manipulating quantum entanglement with atoms and photons in a cavity, Rev. Mod. Phys. 73 (2001) 565, https://doi.org/ 10.1103/RevModPhys.73.565.

ˆ [17] L.O. Castanos, Simple, analytic solutions of the semiclassical Rabi model, Opt. Commun. 430 (2019) 176–188, https://doi.org/10.1016/j.optcom.2018.08.046. [18] S. Felicetti, et al., Quantum Rabi model in the Brillouin zone with ultracold atoms, Phys. Rev. A 95 (2017) 013827, https://doi.org/10.1103/PhysRevA.95. 013827. [19] A. Dareau, Y. Meng, P. Schneeweiss, A. Rauschenbeutel, Observation of ultrastrong spin-motion coupling for cold atoms in optical microtraps, Phys. Rev. Lett. 121 (2018) 253603, https://doi.org/10.1103/PhysRevLett.121.253603. [20] P. Schneeweiss, A. Dareau, C. Sayrin, Cold-atom-based implementation of the quantum Rabi model, Phys. Rev. A 98 (2018) 021801(R), https://doi.org/10. 1103/PhysRevA.98.021801. [21] M.J. Hwang, R. Puebla, M.B. Plenio, Quantum phase transition and universal dynamics in the Rabi model, Phys. Rev. Lett. 115 (2015) 180404, https://doi. org/10.1103/PhysRevLett.115.180404. [22] P.K. Aravind, J.O. Hirschfelder, Two-state systems in semiclassical and quantized fields, J. Phys. Chem. 88 (1984) 4788, https://doi.org/10.1021/j150665a002. [23] Y. Wu, X. Yang, Strong-coupling theory of periodically driven two-level systems, Phys. Rev. Lett. 98 (2007) 013601, https://doi.org/10.1103/PhysRevLett.98. 013601. [24] J.C.A. Barata, W.F. Wreszinski, Strong-coupling theory of two-level atoms in periodic fields, Phys. Rev. Lett. 84 (2000) 2112, https://doi.org/10.1103/ PhysRevLett.84.2112. [25] D.F.V. James, J. Jerke, Effective Hamiltonian theory and its applications in quantum information, Can. J. Phys. 85 (2007) 625, https://doi.org/10.1139/p07-060. [26] J.B. Marion, S.T. Thornton, Classical Dynamics of Particles and Systems, 4th edition, Brooks/Cole, 1995. [27] M.H. Holmes, Introduction to Perturbation Methods, 2nd edition, Springer, 2013. [28] E.A. Coddington, R. Carlson, Linear Ordinary Differential Equations, Society for Industrial and Applied Mathematics, 1987.