Thermal state of the general time-dependent harmonic oscillator

P

—journal of July 2003

physics pp. 7–20

Thermal state of the general time-dependent harmonic oscillator

JEONG-RYEOL CHOI

Department of Electronic Physics, College of Natural Science, Hankuk University of Foreign Studies, Yongin 449-791, Korea

Email: choiardor@hanmail.net

MS received 25 July 2002; revised 15 January 2003; accepted 21 February 2003

Abstract. Taking advantage of dynamical invariant operator, we derived quantum mechanical solu- tion of general time-dependent harmonic oscillator. The uncertainty relation of the system is always larger than^~=2 not only in number but also in the thermal state as expected. We used the diagonal elements of density operator satisfying Leouville–von Neumann equation to calculate various ex- pectation values in the thermal state. We applied our theory to a special case which is the forced Caldirola–Kanai oscillator.

Keywords. Time-dependent harmonic oscillator; thermal state; density operator.

PACS Nos 03.65.-w; 05.30.-d; 05.70.-a

1. Introduction

Vibration may be one of the most dominant physical aspect that we come upon in everyday life [1]. For small oscillation, it can be approximated to the motion of harmonic oscillator.

Harmonic oscillator that has time-dependent mass or frequency may be a good example of time-dependent Hamiltonian systems. Although a large number of dynamical systems have been investigated using approximation and perturbation method in the literature [2,3], we confine our concern to the exact quantum solution of the time-dependent system. There are about three kinds of methods to solve the quantum solution of time-dependent har- monic oscillator. These are propagator method [4–6], unitary transformation method [7–9]

and invariant operator method [10–16]. We will use invariant operator method and uni- tary transformation method together to evolve the quantum theory and investigate thermal state of the general time-dependent harmonic oscillator. The time-dependent harmonic oscillator has several applications such as electrical behavior of LC-circuit that has time- dependent parameter [17] and/or RLC-circuit [18], path-integral formulation of real-time finite-temperature field theory [19–21], dissipative quantum tunnelling effect in macro- scopic system [22–25] and quantum motion of an ion in a Paul trap [7,26,27].

When a system interacts with environment, its coupling parameters may explicitly de- pend on time. Even if the system is closed so that it is conserved as a whole, its subsystem

may implicitly depend on time through interaction with the remnant of the system. The main purpose of this paper is to evolve the thermal state of the general time-dependent harmonic oscillator. The Liouville–von Neumann equation [28,29] for non-equilibrium dynamics can be applicable to both time-dependent harmonic and unharmonic oscillator.

The density operator of the system can be obtained using the wave function satisfying Schr¨odinger equation and can be used to derive various expectation values of variables in the thermal state.

In^x2, we investigate quantum mechanical solution of the general time-dependent har- monic oscillator. The thermal state of the system is discussed in ^x3 on the basis of Liouville–von Neumann approach. In^x4, we will apply our theory for a special case which is the forced Caldirola–Kanai oscillator. Finally,^x5 summarizes this paper and concludes about some physical results of the system.

2. Hamiltonian, invariant operator and wave function

The Hamiltonian of general time-dependent harmonic oscillator can be written as

Hˆ⁽xˆ^;pˆ^;t⁾⁼A⁽t⁾pˆ²⁺B⁽t⁾⁽x ˆˆp⁺p ˆˆx⁾⁺C⁽t⁾pˆ⁺D₂⁽t⁾xˆ²⁺D₁⁽t⁾xˆ⁺D₀⁽t^); (1) where A⁽t⁾ C⁽t⁾and D_i⁽t⁾⁽i⁼0^;1^;2⁾are time-dependent coefficients. These coefficients are real and differentiable with respect to t and note that A⁽t⁾⁶⁼0. The corresponding equation of motion can be derived from Hamilton’s equation of motion as

¨ˆx A˙ A˙ˆx⁺

2 ˙AB

A 2 ˙B 4B²⁺4AD₂

ˆ x⁺AC˙

A C˙ 2BC⁺2AD₁⁼0^: (2) The introduction of invariant operator may save the labor of finding quantum mechanical solution of the system. We let the trial invariant operator as the form

Iˆ⁽t⁾⁼α₁^{(t)[pˆ p}p⁽t^)]2+α₂^{(t)f[xˆ x}p⁽t^)][pˆ p_p⁽t^)]

+[pˆ p_p⁽t^)][xˆ x_p⁽t^)]g+α₃^{(t)[xˆ x}p⁽t^)]2; (3) whereα₁^(t) α₃^(t)are time-variable functions which should be determined afterwards and x_p⁽t⁾is a particular solution of the equation of motion in ˆx space, (eq. (2)) and p_p⁽t⁾is the corresponding particular solution of the equation of motion in ˆp space. We can choose the dimension of ˆI⁽t⁾the same as that of the Hamiltonian. By virtue of its definition, the invariant operator must satisfy the following relation:

d ˆI⁽t⁾ dt ⁼

∂^Iˆ(t)

∂^{t +} 1

i^~[Iˆ⁽t^);Hˆ^]=0^: (4)

Substituting eqs (1) and (3) in the above equation gives

α₁^(t)=c₁ρ₁^2(t)+c₂ρ₁^(t)ρ₂^(t)+c₃ρ₂^{2(t); (5)} α₂^{(t)= 1}

4A^f4^[c₁ρ₁^2(t)+c₂ρ₁^(t)ρ₂^(t)+c₃ρ₂^2(t)]B

[2c₁ρ₁^(t)ρ^˙₁^(t)+c₂ρ^˙₁^(t)ρ₂^(t)+c₂ρ^˙₂^(t)ρ₁^(t)+2c₃ρ₂^(t)ρ^˙₂^{(t)]g; (6)}

α₃^{(t)= 1} 2A²

f

1

2^[c₁ρ˙₁^2(t)+c₂ρ^˙₁^(t)ρ^˙₂^(t)+c₃ρ^˙₂^2(t)]

B^[2c₁ρ₁^(t)ρ^˙₁^(t)+c₂ρ^˙₁^(t)ρ₂^(t)+c₂ρ^˙₂^(t)ρ₁^(t)+2c₃ρ₂^(t)ρ^˙₂^(t)]

+2B^2[c₁ρ1²⁽t⁾⁺c₂ρ₁^(t)ρ₂^(t)+c₃ρ2²⁽t^)]g; (7) where c₁–c₃are constants and ¨ρ₁;2⁽t⁾are two independent solutions of the following dif- ferential equation

ρ¨₁;2⁽t⁾ A˙

Aρ˙₁;2⁽t⁾⁺

2 ˙AB

A 2 ˙B 4B²⁺4AD₂

ρ₁;2⁽t⁾⁼0^: (8) To simplify the invariant operator, we introduce the unitary operator defined as

Uˆ_t⁼Uˆ⁰⁰Uˆ⁰Uˆ^; (9)

Uˆ ⁼exp

i

~

x_ppˆ

exp

i

~

p_pxˆ

; (10)

Uˆ⁰⁼exp

i α₂ 2α₁^~xˆ2

; (11)

Uˆ⁰⁰⁼exp

i 4^~

(x ˆˆp⁺p ˆˆx⁾ln⁽2α₁⁾

: (12)

We can transform the invariant operator using the above operator as

Iˆ⁰⁼Uˆ_tI ˆˆU_t^†: (13)

Then, ˆI⁰reduces to the following simple form Iˆ⁰⁼1

2pˆ²⁺1

2ω^{2xˆ2; (14)}

where

ω²⁼⁴⁽α₁α₃ α2²⁾⁼

1

4A²⁽ρ₁ρ^˙₂ ρ^˙₁ρ₂^)2(4c₁^c₃ ^c22⁾⁼constant^: (15) For convenience, we only discuss the system forω^2>0. Since eq. (14) is the same as that of the Hamiltonian of the simple harmonic oscillator with unit mass, we can introduce the ladder operators defined as

ˆb⁼

rω

2^~xˆ⁺ i

p

2ω^{~pˆ; (16)}

ˆb^†=

rω

2^~xˆ i

p

2ω^{~pˆ: (17)}

These satisfy the boson commutation relation^[ˆb^;ˆb^†]=1. In terms of these operators, eq.

(14) can be simplified to

Iˆ^0=~ω

ˆb^†ˆb⁺1 2

: (18)

The eigenvalue equation for ˆI⁰can be written as

Iˆ^0jn⁰⁽t⁾ⁱ⁼λn^jn⁰⁽t^)i: (19)

We can easily identify the eigenstate in ˆx space as

hxˆ^jn⁰⁽t⁾ⁱ⁼

ω

~π

1⁼4 1

p

2ⁿn!Hn

rω

~

ˆ x

exp

ω

2^~xˆ²

: (20)

The eigenstate of the untransformed invariant operator can be obtained from

hxˆ^jn⁽t⁾ⁱ⁼Uˆ_t^†hxˆ^jn⁰⁽t^)i: (21) Substituting eqs (9) and (20) into the above equation gives

hxˆ^jn⁽t⁾ⁱ⁼

ω

2α₁^~π

1⁼4

1

p

2ⁿn!H_n

r ω

2α₁^{~(xˆ xp)}

exp

i

~

p_pxˆ

exp

1 2α₁^~

ω 2 ⁺iα₂

(xˆ x_p⁾²

: (22)

The ˆx space Schr¨odinger solution^hxˆ^jψnⁱof the Hamiltonian, eq. (1), is the same as the eigenstate of ˆI, except for some time-dependent phase factor,εn⁽t⁾[30]

hxˆ^jψn⁽t⁾ⁱ⁼exp^[iεn⁽t^)]hxˆ^jn⁽t^)i: (23) Inserting the above equation into Schr¨odinger equation, we derive the relation

~ε˙n⁽t^)=hn⁽t^)j

i^~∂

∂^{t Hˆ}

jn⁽t^)i: (24)

Using eqs (1), (22) and (24),εn⁽t⁾can be obtained as εn⁽t⁾⁼ ω

n⁺1 2

Zt 0

A⁽t⁰⁾ α₁ ^dt

0

1

~ Zt

0

H_p⁽x_p⁽t^0);p_p⁽t^0);t⁰⁾dt^0; (25) where Hp⁽xp⁽t^);pp⁽t^);t⁾is defined as

H_p⁽x_p⁽t^);p_p⁽t^);t⁾⁼A⁽t⁾p²_p⁽t⁾⁺C⁽t⁾p_p⁽t⁾ D₂⁽t⁾x²_p⁽t⁾⁺D₀⁽t^): (26) Substituting eq. (25) into (23), we can obtain the exact wave function as

hxˆ^jψn⁽t^)i=hxˆ^jn⁽t⁾ⁱexp

iω

n⁺1 2

Zt 0

A⁽t⁰⁾ α₁ ^dt

0

i

~ Zt

0

H_p⁽x_p⁽t^0);p_p⁽t^0);t⁰⁾dt⁰

: (27)

The ˆp space wave function is related to the ˆx space by the Fourier transformation

hpˆ^jψn⁽t⁾ⁱ⁼ 1

p

2π^~

Z ∞

∞

hxˆ^jψn⁽t⁾ⁱexp

ip ˆˆx

~

d ˆx^: (28)

Using (28), the above equation can be calculated as

hpˆ^jψn⁽t⁾ⁱ⁼⁽ i⁾ⁿ

2ωα₁

~π

1⁼4 1

p

2ⁿn!

(ω ²ⁱα₂⁾ⁿ

(ω⁺²ⁱα₂⁾ⁿ⁺¹

1⁼2

H_n

"s

2α₁ω

~(ω²⁺⁴α₂^{2)(pˆ pp)}

#

exp

i

~

x_p⁽pˆ p_p⁾ α₁⁽pˆ pp⁾2

~(ω⁺²ⁱα₂⁾

exp

iω

n⁺1 2

Zt 0

A⁽t⁰⁾ α₁ ^dt

0

i

~ Zt

0

H_p⁽x_p⁽t^0);p_p⁽t^0);t⁰⁾dt⁰

: (29)

To express ˆI⁽t⁾in a simple form, we introduce another ladder operator as ˆ

a⁽t⁾⁼ 1

pω^~α₁

hω 2 ⁺iα₂

(xˆ x_p⁾⁺iα₁^{(pˆ p}p⁾

i

; (30)

ˆ

a^†(t⁾⁼ 1

pω^~α₁

h ω 2 iα₂

(xˆ x_p⁾ iα₁^{(pˆ p}p⁾

i

: (31)

These operators also satisfy^[aˆ^;aˆ^†]=1. In terms of eqs (30) and (31), (3) can be expressed as

Iˆ⁽t^)=~ω

ˆ a^†aˆ⁺1

2

: (32)

From eqs (30) and (31), we can confirm that the coordinate and the momentum can be expressed as

ˆ x⁼

r

~α₁

ω ^{(aˆ+aˆ†)+xp; (33)} ˆ

p⁼i

s

~

α₁ω

hω 2 ⁺iα₂

ˆ a^†

ω 2 iα₂

ˆ a

i

+p_p^: (34)

Using eqs (27), (33) and (34), we can calculate the following expectation values

hψn^jxˆ^jψnⁱ⁼xp^; (35)

hψn^jpˆ^jψnⁱ⁼p_p^; (36)

hψn^jxˆ^2jψnⁱ⁼

~α₁

ω ^{(2n+1)+x2p; (37)}

hψn^jpˆ^2jψnⁱ⁼

~

α₁ω

ω² 4

+α2²

(2n⁺1⁾⁺p²_p^; (38)

hψn^j(x ˆˆp⁺p ˆˆx^)jψnⁱ⁼

2α₂^~

ω ^{(2n+1)+2xppp: (39)} Then, we can easily identify the uncertainty relation as

∆x∆ˆ pˆ^=[hψn^jxˆ^2jψn^{i (h}ψn^jxˆ^jψnⁱ⁾2

]

1⁼2

[hψn^jpˆ^2jψn^{i (h}ψn^jpˆ^jψnⁱ⁾2

]

1⁼2

=

~

ω

q

ω²⁺⁴α₂²

n⁺1 2

: (40)

This is always larger than^~=2 as expected. The uncertainty relation of q-deformed har- monic oscillator also differs from the uncertainty relation for the simple harmonic oscillator [31].

By performing a similar procedure, we obtain the expectation value of Hamiltonian as

hψn^jHˆ^jψnⁱ⁼

~

ω^[α₁^D₂ ²α₂^B+α₃^A](2n+1)+Hp⁽x_p⁽t^);p_p⁽t^);t^): (41)

3. Thermal state

We consider an ensemble of particles that satisfies the given general time-dependent har- monic oscillator motion. Let us assume that these particles conform to the Bose–Einstein distribution function.

Density operator of the system may satisfy Liouville–von Neumann equation as

∂ρ^ˆ(t)

∂^{t +} 1

i^~[ρˆ^{(t);Hˆ]=0: (42)}

Then, we can express the density operator in ˆx space as ρˆ^{(xˆ;xˆ0;t)= 1}

Z⁽t⁾

∑

hxˆ^jψn⁽t⁾ⁱexp

~ω kT

n⁺1 2

hψn⁽t^)jxˆ^0i; (43) where k is the Boltzmann constant and T the temperature of the system at initial time.

The partition function of the system can be given by Z⁽t⁾⁼

∑

hψn⁽t^)je ^Iˆ(t)=(kT)jψn⁽t^)i: (44) Using eq. (27), the partition function, eq. (44) and the density operator, eq. (43) can be calculated as

Z⁽t⁾⁼ 1

2 sinh^[~ω^{=(2kT)]; (45)}

ρˆ^{(xˆ;xˆ0;t)=}

ω

2π^~α₁^tanh

~ω 2kT

1⁼2

exp

(

iα₂

2^~α₁^{[(xˆ xp)2 (xˆ}

0 x_p^)2]

ω 8^~α₁

(xˆ⁺xˆ⁰ 2x_p⁾²tanh

~ω 2kT

+(xˆ xˆ⁰⁾²coth

~ω 2kT

)

exp

i

~

(xˆ xˆ⁰⁾pp

: (46)

At high temperature, eq. (46) reduces to ρˆ^{(xˆ;xˆ0;t)'}ω

2 1

(α₁π^kT)1=2exp

i

~

(xˆ xˆ⁰⁾p_p

exp

iα₂ 2^~α₁^{[xˆ2 xˆ}

02 2⁽xˆ xˆ⁰⁾x_p^] kT 4^~2α₁^{(xˆ xˆ}

0

)

2

: (47) If the difference between ˆx and ˆx⁰ are sufficiently small compared to⁽2^~α₁⁼α₂^{)1=2 and}

~=p_p, and 1⁼⁽kT⁾approaches zero, the density operator may be simply represented as ρˆ^{(xˆ;xˆ0;t)'~}ω

kT δ^{(xˆ xˆ0): (48)}

On the other hand, at low temperature, it becomes ρˆ^{(xˆ;xˆ0;t)'}

ω

2π^~α₁

1⁼2

exp

i

~

p_p⁽xˆ xˆ⁰⁾

exp

(

iα₂ 2^~α₁^{[xˆ2 xˆ}

02 2⁽xˆ xˆ⁰⁾x_p^] ω

4^~α₁^(xˆ2+xˆ

02

+2x²_p 2 ˆxxp 2 ˆx⁰xp⁾

)

: (49)

The diagonal element of the density operator, eq. (46), can be written as f⁽xˆ⁾⁼

ω

2π^~α₁^tanh

~ω 2kT

1⁼2

exp

ω

2^~α₁^tanh

~ω 2kT

(xˆ xp⁾2

:

(50) The above equation represents the probability that the mass of the oscillator reside at ˆx. As temperature increases, it becomes

f⁽xˆ^)'ω 2

1

(πα₁^kT)1=2exp

ω²

4kTα₁^{(xˆ xp)2}

: (51)

Equation (50) can be used to calculate the expectation value in coordinate space as

hxˆ^li_T⁼

Z∞

∞xˆ^lf⁽xˆ⁾d ˆx^: (52)

For example, we can obtain for l⁼1^;2 as

hxˆⁱ_T⁼x_p^; (53)

hxˆ²ⁱ_T^=~α₁ ω ^coth

~ω 2kT

+x²_p^: (54)

When considering eq. (44), the expectation value of ˆI in the thermal state can be derived from

hIˆⁱ_T⁼kT² ∂

∂^{Tln Z(t): (55)}

Making use of eq. (45), the above equation becomes

hIˆⁱ_T⁼1 2^~ω^coth

~ω 2kT

: (56)

Using the same procedure in the ˆx space, the ˆp space representation of the density operator can be obtained as

ρ⁽pˆ^;pˆ^0;t⁾⁼

2α₁ω

~π⁽ω²⁺⁴α₂^2)tanh

~ω 2kT

1⁼2

exp

i

~

xp⁽pˆ pˆ⁰⁾

exp

(

α₁ω

~(ω²⁺⁴α₂²⁾

(

2iα₂ ω ^{[pˆ2 pˆ}

02 2p_p⁽pˆ pˆ^0)]

1 2

"

(pˆ⁺pˆ⁰ 2p_p⁾²tanh

~ω 2kT

+(pˆ pˆ⁰⁾²coth

~ω 2kT

#))

:

(57) At high temperature, the above equation becomes

ρ^{(pˆ;pˆ0;t)'}

α₁ω² π⁽ω²⁺⁴α₂^2)kT

1⁼2

exp

i

~

x_p⁽pˆ pˆ⁰⁾

exp

(

α₁

~(ω²⁺⁴α₂²⁾

(

2iα₂^{[pˆ2 pˆ02 2p}p⁽pˆ pˆ^0)]

kT

~

(pˆ pˆ⁰⁾²

))

: (58)

On the other hand, at low temperature it can be expressed as ρ^{(pˆ;pˆ0;t)'}

2α₁ω

~π⁽ω²⁺⁴α₂²⁾

1⁼2

exp

i

~

x_p⁽pˆ pˆ⁰⁾

exp

(

α₁ω

~(ω²⁺⁴α₂²⁾

(

2iα₂ ω ^{[pˆ2 pˆ}

02 2p_p⁽pˆ pˆ^0)]

(pˆ²⁺pˆ⁰² 2 ˆpp_p 2 ˆp⁰p_p⁺2p²_p⁾

))

: (59)

The probability that the mass of the oscillator resides at ˆp is obtained taking the diagonal elements of eq. (57) as

f⁽pˆ⁾⁼

2α₁ω

~π⁽ω²⁺⁴α₂^2)tanh

~ω 2kT

1⁼2

exp

2α₁ω

~(ω²⁺⁴α₂^2)tanh

~ω 2kT

(pˆ p_p⁾²

: (60)

Using eq. (60), we can calculate the expectation value in ˆp space as

hpˆⁱ_T⁼p_p^; (61)

hpˆ²ⁱ_T^=~(ω²⁺⁴α₂²⁾ 4α₁ω ^coth

~ω 2kT

+p²_p^: (62)

We can also write the expectation value of ˆI in the thermal state as

hIˆⁱ_T⁼α₁^(hpˆ2i_T ^p2p⁾⁺2α₂^{(hx ˆˆpi}_T ^xpp_p⁾⁺α₃^(hxˆ2i_T ^x2p^): (63) Substituting eqs (54), (56) and (62) into the above equation gives

hx ˆˆpⁱ_T^{= ~}ω 2α₂

1 4

1

ω²⁽α2²⁺α₁α₃⁾

coth

~ω 2kT

+x_pp_p^: (64) Then, the expectation value of the Hamiltonian, eq. (1), can be calculated as

hHˆⁱ_T^{= ~}

ω⁽α₃^{A 2}α₂^B+α₁^D₂^)coth

~ω 2kT

+H_p⁽x_p⁽t^);p_p⁽t^);t⁾ (65) and the uncertainty relation in thermal state can be calculated as

(∆xˆ∆pˆ⁾_T^=[hxˆ²ⁱ_T ^hxˆⁱ²_T^)(hpˆ²ⁱ_T ^hpˆⁱ²_T^)]1=2

=

~

2ω

q

ω²⁺⁴α₂^2coth

~ω 2kT

: (66)

By comparing the above equation with eq. (40), we can confirm that the uncertainty rela- tion in thermal state varies as time goes by, with the same fashion in number state.

4. Forced Caldirola–Kanai oscillator

We can apply our theory to various kinds of time-dependent Hamiltonian systems. As an example, let us see for the forced Caldirola–Kanai oscillator [32,33]. For this system, the time-dependent coefficients in eq. (1) are given by

A⁽t⁾⁼ 1

2me ^βt; (67)

D₂⁽t⁾⁼1

2mω0²e^βt; (68)

D₁⁽t⁾⁼ F⁽t⁾e^βt; (69)

B⁽t⁾⁼C⁽t⁾⁼D₀⁽t⁾⁼0^; (70)

so that we can rewrite the Hamiltonian as Hˆ⁼e ^βt pˆ²

2m⁺e^βt1

2mω₀²xˆ² e^βtF⁽t⁾xˆ^; (71) where m is the mass,β the damping constant and F⁽t⁾the arbitrary time-dependent driving force. Equation (8) becomes

ρ¨₁;2⁺βρ^˙₁;2⁺ω₀²ρ₁;2⁼0^: (72) The two classical solutions of the above equation can be written as

ρ₁^(t)=ρ₁^{(0)e βt=2eiωt; (73)} ρ₂^(t)=ρ₂^{(0)e βt=2e iωt; (74)} whereωis given by

ω⁼

r

ω₀² β²

4 ^: (75)

We choose c₁–c₃in eqs (5)–(7) as

c₂⁼ 1

2mρ₁⁽⁰⁾ρ₂^{(0); c1=c3=0: (76)}

The particular solutions x_pand p_psatisfy the following relations:

¨

x_p⁺β^x˙p⁺ω₀^2xp⁼

F⁽t⁾ m

; (77)

¨

p_p β^p˙p⁺ω0²p_p⁼e^βtF˙⁽t^): (78)

The solutions of the above equations depend on F⁽t⁾. If, we choose F⁽t⁾as

F⁽t⁾⁼F₀t^; (79)

the solutions of eqs (77) and (78) will be x_p⁽t⁾⁼ F₀

mω₀^2t β^F₀

mω₀^{4; (80)}

p_p⁽t⁾⁼ F₀

ω₀^{2eβt: (81)}

We will also investigate the system driven by the exponentially decaying force:

F⁽t⁾⁼F₀e ^γt; (82)

whereγis an arbitrary real constant. In this case, the particular solutions are given by x_p⁽t⁾⁼ F₀⁼m

γ² βγ⁺ω₀^{2e γt; (83)}

p_p⁽t⁾⁼ γ^F₀ γ² βγ⁺ω₀^2e

(β γ⁾t

: (84)