text:
# Spatial evolution of short laser pulses under coherent population trapping

## I INTRODUCTION

The electromagnetically induced transparency (EIT) can be used to make optically thick media transparent to resonant laser radiation . The EIT is the result of various quantum interference effects such as nonlinear interference , coherent population trapping (CPT) , adiabatic population transfer (APT) . Optical characteristics of the matter undergo drastic changes under those effects to such an extent that they can now be manipulated. A lot of interesting applications based on that have been proposed and experimentally realized (see, e.g. ).

Interesting and unusual phenomena caused by the above indicated effects can be observed when laser pulses propagate in a resonant three-level medium. Propagation of pulses under EIT conditions was studied for example in . As a rule, the situations are considered when both pulses have identical forms and their duration is longer than the relaxation time of the intermediate resonant state (matched pulses ; dressed field pulses ) or when the duration of the coupling radiation considerably exceeds that of the probe radiation (adiabatons ). A theoretical study of certain features of spatial evolution under APT conditions is presented in . Propagation of soliton-like pulses in a three-level system is studied in . A three-level system with equal oscillator strengths is under consideration in all the above mentioned studies, whereas in actual fact the transition oscillator strengths most often are different.

In this paper, the spatial and temporal evolution is studied of two overlapping short laser pulses propagating in a resonant optically thick medium that consists of three-level $`\mathrm{\Lambda }`$-atoms. Pulses of such configuration are widely used to enhance the efficiency of nonlinear generation processes . The two pulses are assumed to have identical shapes but different durations $`(T_2>T_1)`$ as shown in Fig.1. The pulse durations are much shorter than any of the times of relaxation in the medium (short pulses). It is also assumed that the pulse envelopes satisfy the adiabaticity criterion :

$`(G_2\dot{G_1}G_1\dot{G_2})/G^31,`$ (1)
where $`G_{1,2}`$ are the Rabi frequencies of the respective fields, $`G=\sqrt{G_1^2+G_2^2}`$; the dot refers to time derivatives. Condition (1) is easy to satisfy by making one of the pulse amplitudes or both of them large even for short pulses. This will induce a strong coherence at the Raman transition resulting in the effect of CPT. The latter considerably decreases the absorption of the propagating resonant pulses. The dynamics of propagation of such pulses is studied here without restriction of the relationship between oscillator strengths of transitions.

Our theoretical model involves a set of coupled Schrödinger equations and a set of reduced wave equations, allowing thus a simultaneous description of temporal and spatial evolution of the atomic system and the radiation. The equations are analyzed in approximation (1) with the account of the first non-adiabatic correction. It will be shown that the dynamics of propagation strongly depends on the oscillator strengths ratio. Also the spatial and temporal behavior is analyzed of the atomic Raman coherence. A possibility of localizing the atomic coherence spatially has been established. The results obtained are compared to the results reported in .

The paper is organized as follows. In Section II, we describe the model and present the basic equations. Section III contains self-consistent solutions of those equations in the adiabatic limit and describes the temporal behaviour of the level populations and the atomic Raman coherence in the optically thin medium. The spatial evolution of pulses in the optically thick medium is described in Section IV for various oscillator strength ratios. In Section IV, we also discuss the adiabaticity criterion and demonstrate the effect of spatial localization of the atomic coherence. Finally we summarize the results obtained.

## II BASIC EQUATIONS

The three-level system under consideration is shown in Fig.1 together with the temporal configuration of the pulses as they enter the medium. The pulses travel along the same direction, $`z`$. States $`|0`$, $`|1`$ and $`|2`$ are connected by laser pulses $`E_1=1/2E_1(t)exp[i(\omega _1tk_1z)]+c.c.`$ and $`E_2=1/2E_2(t)exp[i(\omega _2tk_2z)]+c.c.`$, respectively. In our further consideration we shall refer to the first pulse $`E_1`$ as the probe and to the second pulse $`E_2`$ as the coupling pulse. The Rabi frequency of probe pulse is comparable with that of the coupling pulse. The pulses are sent simultaneously into an atom. The pulse durations $`T_{1,2}`$ $`(T_2>T_1)`$ are assumed to be much less than any of the relaxation times of atoms. The transition $`|0|2`$ is electric dipole forbidden. The intermediate state $`|1`$ is in a one-photon resonance with each field interacting only with the respective transition.

The following standard set of equations describes the spatial and temporal dynamics of the probability amplitudes of atomic states $`b_{0,1,2}`$ and slowly varying Rabi frequencies $`G_1=d_{10}E_1(t)/2\mathrm{}`$, $`G_2=d_{21}E_2(t)/2\mathrm{}`$ in the local-time coordinate system $`\tau =tz/c`$:

$`{\displaystyle \frac{b_0}{\tau }}`$ $`=`$ $`iG_1^{}b_1\mathrm{exp}(ik_1z),`$ (2)
$`{\displaystyle \frac{b_2}{\tau }}`$ $`=`$ $`iG_2^{}b_1\mathrm{exp}(ik_2z),`$ (3)
$`{\displaystyle \frac{b_1}{\tau }}`$ $`=`$ $`iG_1b_0\mathrm{exp}(ik_1z)+iG_2b_2\mathrm{exp}(ik_2z).`$ (4)
$`{\displaystyle \frac{G_1}{z}}`$ $`=`$ $`iK_1b_1b_0^{}\mathrm{exp}(ik_1z),`$ (5)
$`{\displaystyle \frac{G_2}{z}}`$ $`=`$ $`iK_2b_1b_2^{}\mathrm{exp}(ik_2z).`$ (6)
Here we assumed zero one-photon detunings. $`K_{1,2}=\pi \omega _1\left|d_{10,12}\right|^2N/\mathrm{}c`$ are the propagation coefficients, $`N`$ is the atomic concentration, $`d_{10,12}`$ are the dipole transition matrix elements, $`\omega _{1,2}`$, $`k_{1,2}`$ are the frequencies and the wave numbers of the interacting waves in vacuum, $`c`$ is the light velocity in vacuum. All atoms are assumed to be initially in the ground state $`|0`$: $`b_0(\mathrm{},z)=1`$, $`b_{1,2}(\mathrm{},z)=0`$. We use Gaussian pulses at the medium entrance $`z=0`$ for the purpose of numerical simulation: $`G_1(\tau )=G_1^0\mathrm{exp}(\tau ^2ln2/T_1^2)`$, $`G_2(\tau )=G_2^0\mathrm{exp}[\tau ^2ln2/T_2^2]`$.

In terms of $`a_0=b_0\mathrm{exp}(ik_1z)`$, $`a_2=b_2\mathrm{exp}(ik_2z)`$, $`a_1=ib_1`$, equations (2) and (5) can be written as:

$`{\displaystyle \frac{a_0}{\tau }}`$ $`=`$ $`G_1^{}a_1,`$ (7)
$`{\displaystyle \frac{a_2}{\tau }}`$ $`=`$ $`G_2^{}a_1,`$ (8)
$`{\displaystyle \frac{a_1}{\tau }}`$ $`=`$ $`G_1a_0G_2a_2.`$ (9)
$`{\displaystyle \frac{G_1}{z}}`$ $`=`$ $`K_1a_1a_0^{},`$ (10)
$`{\displaystyle \frac{G_2}{z}}`$ $`=`$ $`K_2a_1a_2^{}.`$ (11)
The coupled equations (7) and (10) give a complete semiclassical description of the resonant different-wavelengths propagation problem we are dealing with.

## III TEMPORAL DYNAMICS OF LEVEL POPULATIONS AND RAMAN COHERENCE IN THE ADIABATIC APPROXIMATION (OPTICALLY THIN MEDIUM)

In this section we study the temporal dynamics of populations and the atomic coherence in the given time-dependent field, assuming that the medium is optically thin. Based on that, $`G_{1,2}`$ will not depend on the coordinate $`z`$. One can show that condition (1) for Gaussian pulses reduces to $`G_2^0T_11`$ when $`T_2/T_1>\sqrt{2}`$. With the first non-adiabatic correction, the solution of Eq. (7) takes on the form:

$`a_0`$ $``$ $`{\displaystyle \frac{G_2(\tau )}{G(\tau )}},a_2{\displaystyle \frac{G_1(\tau )}{G(\tau )}},`$ (12)
$`a_1`$ $``$ $`{\displaystyle \frac{1}{G_1}}{\displaystyle \frac{(G_2/G)}{\tau }}{\displaystyle \frac{1}{G_2}}{\displaystyle \frac{(G_1/G)}{\tau }},`$ (13)
where $`G(\tau )=\sqrt{G_1^2(\tau )+G_2^2(\tau )}`$.

The solutions for probability amplitudes are convenient to be represented as:

$`a_0=\mathrm{cos}\theta (\tau ),a_2=\mathrm{sin}\theta (\tau ),`$ (14)
where the mixing angle $`\theta (\tau )`$ is defined as $`\mathrm{tan}\theta (\tau )=G_1(\tau )/G_2(\tau )`$ (we shall discuss its meaning later on).

The expression for $`a_1`$ can be reduced to:

$`a_1=(G_2\dot{G_1}G_1\dot{G_2})/G^3=\dot{\theta }/G.`$ (15)
In adiabatic limit (1) $`|a_1|=|\dot{\theta }/G|1`$ ($`\dot{\theta }=\theta /\tau `$), i.e. the population of the intermediate state $`|1`$ is close to zero all the time during the interaction with pulses. This also implies that the resonant absorption of the light pulses is weak (electromagnetically induced transparency) and the population is mainly distributed between the initial $`|0`$ and the final $`|2`$ states:

$`|a_0|^2+|a_2|^21.`$ (16)
Equality (16) reflects the fact that atoms are trapped in the CPT state: $`a_{cpt}=(G_2/G)a_0(G_1/G)a_2=a_0\mathrm{cos}\theta a_2\mathrm{sin}\theta =1`$. This effect is responsible for the decrease in the resonant absorption of propagating pulses. Also the Raman coherence $`\rho _{20}=a_0a_2^{}`$ occurs:

$`\rho _{20}={\displaystyle \frac{1}{2}}\mathrm{sin}(2\theta )\text{or}\rho _{20}={\displaystyle \frac{G_1G_2}{G_1^2+G_2^2}}.`$ (17)
Obviously, the maximum coherence (in absolute value) $`|\rho _{20}|=1/2`$ is reached when $`\theta =\pi /4`$ ($`G_1^0=G_2^0`$). Figure 2 shows the temporal behavior of level populations $`|a_{0,2}(\tau )|^2`$, the atomic Raman coherence $`|\rho _{20}(\tau )|`$ and the mixing angle $`\theta (\tau )`$ for Gaussian pulses in the optically thin medium.

The above results can be interpreted in terms of the three-dimensional vector model where vector variables $`\stackrel{}{a}=(a_0,a_2,a_1)`$ and $`\stackrel{}{G}=(G_2,G_1,0)`$ (the torque vector) are introduced. Using these variables, we can rewrite (7) as

$`\dot{\stackrel{}{a}}=\stackrel{}{G}\times \stackrel{}{a},`$ (18)
where the sign ’$`\times `$’ means the vector product. The solution of equation (18) is vector $`\stackrel{}{a}=(G_2/G,G_1/G,\dot{\theta }/G)`$. Components of the vector $`\stackrel{}{a}`$ coincide with adiabatic solution (12). Figure 3 illustrates the dynamics of vectors $`\stackrel{}{G}`$ and $`\stackrel{}{a}`$ in the three-dimensional vector model. The torque vector $`\stackrel{}{G}`$ moves in the $`\stackrel{}{e}_1`$-$`\stackrel{}{e}_2`$ plane. And vector $`\stackrel{}{a}`$ having a small angle with respect to vector $`\stackrel{}{G}`$ ($`|\dot{\theta }/G|1`$) follows it. Thus one can see an absolute analogy with the adiabatic following in the case of a light pulse interacting with a two-level atom . Such a simple picture can be observed only in optically thin media. In optically dense media, $`G_{1,2}`$ and hence $`\theta `$ become dependent on the $`z`$ coordinate.

## IV SPATIAL EVOLUTION OF INTERACTING PULSES IN OPTICALLY DENSE MEDIA

### A General case: unequal oscillator strengths $`(K_1K_2)`$

The condition $`K_1K_2`$ means that the probability of the $`|0|1`$ transition is not equal to that of the $`|2|1`$ transition. We note that in the ideal adiabatic limit $`a_1=0`$, and the pulses would not change their shape as they propagate in a medium that is optically thick for each of the pulses (see (10)). However this is not the case. The non-adiabatic correction has to be introduced for real situations, which results in the induced dipole moments at the transitions $`|1|0`$ and $`|1|2`$ , and in the change of both pulses traveling in the medium. In order to attribute this effect to propagation of the interacting pulses in an optically dense medium, it is necessary to solve Eqs. (7) and (10) in a self-consistent way.

Use (14) and (15) to rewrite field equations (10) in the following form

$`{\displaystyle \frac{G_1}{z}}=K_1{\displaystyle \frac{\dot{\theta }}{G}}\mathrm{cos}\theta ,{\displaystyle \frac{G_2}{z}}=K_2{\displaystyle \frac{\dot{\theta }}{G}}\mathrm{sin}\theta .`$ (19)
From (19), one can show that

$`K_2G_1^2(\tau ,z)`$ $`+`$ $`K_1G_2^2(\tau ,z)`$ (20)
$`=`$ $`K_2G_1^2(\tau ,z=0)+K_1G_2^2(\tau ,z=0),`$ (21)
i.e. $`K_2G_1^2(\tau ,z)+K_1G_2^2(\tau ,z)`$ does not depend on the $`z`$ coordinate.

Equation (20) describes the Manley-Raw relation, i.e. the law of conservation of the total energy density during propagation under CPT conditions.

Using the definitions of $`\theta `$ and $`G`$, we obtain the following expressions for $`G_{1,2}(\tau ,z)`$

$`G_1(\tau ,z)=G(\tau ,z)\mathrm{sin}\theta (\tau ,z),`$ (22)
$`G_2(\tau ,z)=G(\tau ,z)\mathrm{cos}\theta (\tau ,z).`$ (23)
Substitution (22) into the Manley-Raw relation gives the following expression for $`G(\tau ,z)`$:

$`G^2(\tau ,z)=G_0^2(\tau ){\displaystyle \frac{K_2\mathrm{sin}^2(\theta _0(\tau ))+K_1\mathrm{cos}^2(\theta _0(\tau ))}{K_2\mathrm{sin}^2(\theta (\tau ,z))+K_1\mathrm{cos}^2(\theta (\tau ,z))}},`$ (24)
where $`G_0^2(\tau )=G_1^2(\tau ,0)+G_2^2(\tau ,0)`$ and $`\theta _0(\tau )=\theta _0(\tau ,0)`$ are the functions at the medium entrance $`z=0`$.

So the dynamics of level populations as well as the atomic Raman coherence and the evolution of the pulse shape are completely determined by the function $`\theta (\tau ,z)`$ that depends on both the time $`\tau `$ and the coordinate $`z`$. Differentiating $`\mathrm{tan}\theta =G_1/G_2`$ on the $`z`$ and using (19) we obtain the following equation for $`\theta (\tau ,z)`$:

$`{\displaystyle \frac{\theta }{\tau }}+{\displaystyle \frac{G^2(\tau ,z)}{K(\theta )}}{\displaystyle \frac{\theta }{z}}=0,`$ (25)
where $`K(\theta (\tau ,z))=K_1\mathrm{cos}^2(\theta (\tau ,z))+K_2\mathrm{sin}^2(\theta (\tau ,z))`$.

Eq.(25) is similar to the equations describing nonlinear waves with the sharpening of the wave front during propagation . The parameter $`u=G^2/K`$ can be treated as the ”nonlinear” velocity. The nonlinear velocity can be described as $`u(\tau ,z)=A(\tau )/K^2(\theta (\tau ,z))`$ where the first factor $`A(\tau )=G_0^2(\tau )[K_2\mathrm{sin}^2(\theta _0(\tau ))+K_1\mathrm{cos}^2(\theta _0(\tau ))]`$ is independent of the $`z`$ coordinate, and the second factor $`K(\theta (\tau ,z))`$ is maintained along a characteristic of Eq.(25) $`\theta (\tau ,z)=const`$. This allows to write down the characteristic curve equation in an evident form:

$`z={\displaystyle \frac{1}{K^2(\theta _0)}}{\displaystyle _{\tau _0}^\tau }A(\tau ^{})𝑑\tau ^{}`$ (26)
Here $`\tau _0`$ is the time at which the characteristic curve goes out of the medium boundary.

The solution for $`\theta (\tau ,z)`$ has the form:

$`\theta (\tau ,z)=\theta _0(\tau _0,0)`$ (27)
Here $`\tau _0`$ is to be determined from Eq.(26).

### B The case of equal oscillator strengths $`(K_1=K_2)`$

In the case of equal oscillator strengths $`K_1=K_2K`$ the $`G(z,\tau )`$ function is not subject to changes during propagation ($`G(\tau ,z)=G_0(\tau )`$ – the Manley-Raw relation for this case). Therefore Eq. (25) substantially simplifies:

$`{\displaystyle \frac{\theta }{\tau }}+{\displaystyle \frac{G_0^2(\tau )}{K}}{\displaystyle \frac{\theta }{z}}=0.`$ (28)
The solution of Eq. (28) can be written in the following form:

$`\theta (\tau ,z)=\theta _0(Z^1(Z(\tau )z),0),`$ (29)
where $`Z(\tau )=K^1_{\mathrm{}}^\tau 𝑑\tau ^{^{}}G^2(0,\tau ^{^{}})`$, $`Z^1(\tau )`$ is the inverse function of $`Z(\tau )`$.

It is not difficult to show that in this case $`G=\sqrt{G_1^2(\tau ,z)+G_2^2(\tau ,z)}`$ coincides with the definition of dressed field pulses : $`G_{}=a_0G_2a_2G_1`$. Thus the pulses in our case can be identified as dressed field pulses (only at $`K_1=K_2`$). It would be interesting to note that other combination $`G_+=a_0G_1+a_2G_20`$ (see also ). The concept of dressed field pulses cannot be applied to the case of $`K_1K_2`$, but in both cases the Manley-Raw relation remains valid.

### C The adiabaticity criterion

The above results were obtained on the assumption that the adiabaticity criterion (1) (or $`|\dot{\theta }/G|1`$) remains valid during propagation of pulses. However it is not necessary the case. Therefore we investigate the adiabaticity condition with the account of propagation. Differentiating (27) on $`\tau `$, we can write the following expression for the adiabaticity criterion:

$`{\displaystyle \frac{\dot{\theta }(\tau ,z)}{G(\tau ,z)}}`$ $`=`$ $`{\displaystyle \frac{\theta _0}{\tau _0}}{\displaystyle \frac{G(\tau ,z)}{G_0^2(\tau _0)}}{\displaystyle \frac{K(\theta (\tau ,z))}{K(\theta _0(\tau _0))}}`$ (30)
$`\times `$ $`[1+\mathrm{sin}(2\theta _0(\tau _0)){\displaystyle \frac{2(K_2K_1)z}{G_0^2(\tau _0)}}{\displaystyle \frac{\theta _0}{\tau _0}}]^11`$ (31)
As follows from (30), the adiabaticity condition is destroyed ($`\dot{\theta }/G\mathrm{}`$) when the factor in the square brackets tends to zero. Since $`\mathrm{sin}(2\theta _0(\tau _0))>0`$ in the entire range of change of $`\theta `$ ($`0\theta \pi /4`$), relation (30) is not fulfilled under the following conditions

$`{\displaystyle \frac{2(K_1K_2)z}{G_0^2(\tau _0)}}{\displaystyle \frac{\theta _0}{\tau _0}}\mathrm{sin}(2\theta _0(\tau _0))=1,`$ (32)
$`(K_1K_2){\displaystyle \frac{\theta _0}{\tau _0}}>0,`$ (33)
Evidently, condition (32) is not fulfilled at $`K_1=K_2`$, and the adiabaticity criterion holds throughout the propagation process in this case. It should be mentioned that numerical analysis of Eq. (32) reveals that, generally speaking, there is a range of change of the $`q=K_1/K_2`$ parameter when the first condition in (32) is not fulfilled either: $`q_{min}<q<q_{max}`$, $`q_{min}<1`$, $`q_{max}>1`$. The values $`q_{max}`$ and $`q_{min}`$ depend on the ratio $`a=T_2/T_1`$ and the shape of pulses. For example, the $`q`$ parameter for Gaussian pulses with $`T_2/T_1=3`$ may vary within the limit $`0.40<q<1.35`$, as shown by the numerical simulation.

In this case for the thick medium the adiabaticity criterion has the form:

$`{\displaystyle \frac{\theta _0}{\tau _0}}{\displaystyle \frac{G(\tau ,z)}{G_0^2(\tau _0)}}1.`$ (34)
It can be readily shown that condition (34) for Gaussian pulses is fulfilled when $`T_2/T_1>\sqrt{3}`$ and $`G_2^0T_11`$ (compare with the case of the thin medium).

For Gaussian pulses ($`T_2>T_1`$) and the initial conditions $`a_0(\mathrm{})=1`$, $`a_2(\mathrm{})=0`$, we have $`\dot{\theta _0}>0`$ for $`\mathrm{}<\tau _0<0`$. With $`K_1>K_2`$ ($`q>q_{max}`$) the adiabaticity criterion begins to break down at the leading edge: the mixing angle front becomes steeper. Condition (30) gets destroyed at the trailing edge (where $`\dot{\theta _0}<0`$) when $`K_1<K_2`$ ($`q<q_{min}`$). Let s introduce a critical length $`z_c`$ at which the adiabaticity condition (30) cannot be satisfied. Using (32), one can obtain the following simple estimation for $`z_c`$ for the case $`K_1K_2`$:

$`z_c{\displaystyle \frac{G_0^2T_1}{2|K_1K_2|}}.`$ (35)
Let us now consider the experimental parameters: $`N=10^{15}cm^3`$; $`1/\lambda _{1,2}=10000,20000(cm^1)`$; $`1/\gamma _{1,2}=10,50ns`$ are the relaxation times of levels $`|1`$ and $`|2`$, respectively; $`T_1=0.1ns`$; $`G_{1,2}^0T_1=20`$; $`a=T_2/T_1=3`$; $`fg_{10}=0.1`$ is the oscillator strength of the probe transition. From (35) we obtain $`z_c=10^5z_02cm`$, where $`z_0`$ is the linear absorption length. This estimation agrees with the results presented in Figs. 5-7.

Using the characteristic equation (26) it is easy to find the area where the adiabaticity criterion breaks down. Characteristic curves for various values of $`q=K_1/K_2`$ are shown in Fig 4. The thickening of the characteristic curves means the sharpening of the mixing angle front at a certain medium depth. At the point of intersection (see Fig. 4 b,c), $`\dot{\theta }\mathrm{}`$, condition (30) collapses. However there are no such points in the case when $`q_{min}<q<q_{max}`$ (see Fig.4a), i.e. here the adiabaticity criterion (30) is maintained during propagation of pulses.

The above discussion is illustrated by Fig.5 where the temporal behavior of the mixing angle $`\theta `$ is presented for $`q=K_1/K_2=1,0.25,4`$ and for different normalized propagation lengths. On the Fig.5a-b the analytical results are presented (formulae (27) and (29)). Figure 5c shows the numerical solution for $`\theta `$ at $`K_1>K_2`$.

One can see that the evolution of the $`\theta `$ parameter at $`K_1>K_2`$ is different from that at $`K_1<K_2`$. At $`K_1>K_2`$ the adiabaticity condition fails for all values of $`z`$, beginning from the critical lengths $`z_c`$ defined by Eq. (35). Here, the analytical theory does not apply at the very late stage of nonlinear wave propagation. In the case of $`K_1<K_2`$ , nonadiabaticity develops at the trailing edge (the front becomes steeper at a certain propagation length), but it does not go deeper into the medium. A good agreement between the analytical and numerical solutions for $`\theta `$ at $`K_1K_2`$ is observed over the entire propagation length. This leads us to conclude that the interaction adiabaticity is fairly sensitive to the ratio of the oscillator strengths of the transitions interacting with the pulses.

### D Discussion of results

The solutions obtained have been used to analyze the temporal dynamics and spatial behavior of the propagating EIT pulses and the atomic coherence for various oscillator strength ratios. Figure 6 illustrates the temporal and spatial evolution of normalized Rabi frequencies of the both pulses $`g_{1,2}(\tau )=G_{1,2}/\sqrt{G_{1}^{0}{}_{}{}^{2}+G_{2}^{0}{}_{}{}^{2}}`$ as they propagate inside the optically thick medium. The temporal evolution of pulses can be seen to depend on the ratio between the transition oscillator strengths. In the case of $`K_1K_2`$ both pulses undergo reshaping as they propagate in the medium (Fig.6a and b). The probe pulse is gradually depleted and the coupling gets stronger. Note that the pulse shape at the initial stage of propagation shows very little change along the length of the medium, which may exceed the linear absorption length. Complete reemitting of the probe pulse into the coupling one during propagation is possible. Using the Eq. (26) one can obtain the following expression for the maximal distance $`z_m`$ on which the probe pulse propagates into the medium:

$`z_m={\displaystyle \frac{1}{K^2(\theta _0(\mathrm{}))}}{\displaystyle _{\mathrm{}}^{\mathrm{}}}A(\tau )𝑑\tau .`$ (36)
For the Gaussian pulses we have

$`z_m=\sqrt{{\displaystyle \frac{\pi }{2ln2}}}(G_2^0)^2T_1{\displaystyle \frac{K_2(G_1^0/G_2^0)^2+K_1a}{K_1^2}}.`$ (37)
Here $`a=T_2/T_1`$.

For the same parameters as in the previous sub-section, we have $`z_m=4.810^5z_09cm`$ for $`K_1/K_2=1`$ and $`z_m=8.410^5z_016cm`$ for $`K_1/K_2=1/4`$. These values agree with the results shown in Figs. 5-7.

An interesting feature of the spatial distribution of the probe pulse is illustrated in Fig. 6a,b: in some areas of the medium, the field in the tail of the pulse is different from zero. This is believed to be the result of the spatial compression of the probe caused by the slowing down of the group velocity of the probe pulse .

In the case of $`K_1>K_2`$, the adiabaticity condition is maintained over the $`z_c`$ range, which can also be much longer than the length of the linear absorption. The leading edge of the probe pulse undergoes gradual depletion and the pulse amplitudes display only small changes in that range. Outside that range, the pulse splits into several peaks. This occurs due to the nonadiabaticity of interaction (see Fig. 6c where numerical simulations are presented for $`g_{1,2}(\tau )`$. Here the first three curves correspond to the adiabatic interaction).

The pulse evolution, described above, is determined by spatial and temporal behaviour of the atomic Raman coherence $`\rho _{20}`$ shown in Fig.7 (see also ). Figures 7a,b reveal an unusual spatial and temporal behavior of the atomic coherence, which we interpret as the slowing down, stopping and localization of the atomic coherence in the medium. The probe pulse is transferred into and stored in the collective atomic excitation under the control of the coupling. One can say that a phase grating is created in the atomic medium. The phase grating is preserved throughout the entire period of relaxation of the atomic coherence. During that period of time, the information stored in the atomic excitations can be transferred back to the radiation using another coupling pulse of the same or of a different frequency.

We find these effects to be similar to the ones predicted and demonstrated in , but for some differences as indicated below. In our case, both the coupling pulse and the probe pulse are strong. Unlike , we used the boundary conditions for the pulse envelopes yielding the time distribution of pulses at the medium boundary, $`z=0`$. We believe, this condition is more natural than the one used in  where the authors use the probe pulse distribution in the medium at fixed time as the initial condition. Also they do not take into consideration evolution of the coupling pulse. The effect takes place in the case $`K_1K_2`$ and is not observed in the case of $`K_1>K_2`$.

## V CONCLUSION

The propagation of two short overlapping pulses with durations $`T_2>T_1`$ in optically thick three-level media under CPT conditions has been studied for the general case of unequal transition oscillator strengths. An analytical solution has been obtained for the set of reduced wave equations under the adiabatic following condition. Also it has been shown how spatial evolution of pulses depends on the oscillator strength ratio.

The condition of adiabaticity provided at the medium entrance preserves for any value of propagation lengths if $`K_1K_2`$ ($`q<q_{max}`$) and breaks down at $`K_1>K_2`$ ($`q>q_{max}`$). In the range $`q<q_{max}`$, the probe pulse is completely depleted and reemitted into the coupling pulse during propagation. This is not possible in the case of $`q>q_{max}`$. It has been established that to provide for the adiabaticity condition in an optically thin medium, the restriction $`T_2/T_1>\sqrt{2}`$ has to be ensured, whereas in a thick medium $`T_2/T_1>\sqrt{3}`$.

We also have studied the spatial behavior of the atomic coherence $`\rho _{20}`$, which plays a significant role, for example, in nonlinear mixing processes. It has been found that a strong coherence can be maintained over a length equal to several hundreds of thousand of one-photon absorption lengths during propagation. The effect of localization of the atomic coherence is demonstrated.

FIGURES

FIG. 1. (a) The three-level Lambda-type system coupled by two resonant pulses with Rabi frequencies $`G_1`$ and $`G_2`$. (b) The shapes of probe $`G_1`$ and coupling $`G_2`$ pulses at the medium entrance.

FIG. 2. The time evolution of (a) level populations $`|a_{0,2}(\tau )|^2`$ and the atomic Raman coherence $`|\rho _{20}(\tau )|`$, (b) the mixing angle $`\theta (\tau )`$ in an optically thin medium for the Gaussian pulses. $`G_1^0T_1=G_2^0T_1=20`$, (a) $`T_2/T_1=3`$; (b) $`T_2/T_1=\sqrt{2}`$ (A), $`\sqrt{3}`$ (B), $`3`$ (C), $`10`$ (D).

FIG. 3. The vector model of adiabatic interaction of two short pulses with a three-level Lambda-type system.

FIG. 4. The characteristic curves for the Eq.(25): (a) $`K_1/K_2=1`$, (b) $`K_1/K_2=0.25`$, (c) $`K_1/K_2=4`$.

FIG. 5. The time evolution of the parameter $`\theta `$ for different relationships between $`K_1`$ and $`K_2`$ at different propagation lengths. (a) $`K_1=K_2`$, (b) $`K_1/K_2=0.25`$, (c) $`K_1/K_2=4`$. $`G_{1,2}^0T_1=20`$, $`T_2/T_1=3`$. Here and in all the other figures the time $`\tau `$ is measured in the units of the pulse duration $`T_1`$, and the propagation length $`z`$ of pulses in the medium is measured in the units of the length of linear absorption of the probe radiation determined in accordance with the Beer’s low. In (c) the numerical solution for the case $`K_1/K_2=4`$ is presented.

FIG. 6. The time evolution of the normalized Rabi frequencies $`g_{1,2}=G_{1,2}/\sqrt{G_{1}^{0}{}_{}{}^{2}+G_{2}^{0}{}_{}{}^{2}}`$ of the probe and coupling pulses for different relationships between $`K_1`$ and $`K_2`$ at different propagation lengths. (a) $`K_1=K_2`$, (b) $`K_1/K_2=0.25`$, (c) $`K_1/K_2=4`$. $`g_{1,2}^0=20`$, $`T_2/T_1=3`$.

FIG. 7. The time evolution of the atomic Raman coherence for different relationships between $`K_1`$ and $`K_2`$ at different propagation lengths. (a) $`K_1=K_2`$, (b) $`K_1/K_2=0.25`$, (c) $`K_1/K_2=4`$. $`G_{1,2}^0T_1=20`$, $`T_2/T_1=3`$.


format: text

text:
# Stability analysis of cosmological models through Liapunov’s method

## 1 Introduction

A major difficulty in analyzing general relativistic problems, namely in cosmology and black hole physics, stems from the fact that the relevant field equations are nonlinear. This seriously limits the possibility of obtaining exact solutions, and makes it difficult to assess the degree of generality of the behaviour and special features of those exact solutions which are known. On the other hand, there is an increasing realization of the importance of the asymptotic behaviour of models, which provides the relevant features to be compared with the physical data available for each era of the universe.

As a consequence, in recent years the emphasis in the study of cosmological models has drifted from the search for exact solutions to the analysis of the qualitative properties of the equations and of the long term behaviour of the solutions . The main tools used in this context have been perturbation methods and qualitative theory, especially linear stability analysis. Liapunov’s method has been used in relatively few instances  as a means to prove stability.

In this paper, we argue in favour of the use of Liapunov’s method for the stability analysis of the possible attractors in a large class of cosmological models, as an alternative to linear approximation. One main advantage of the method is that it provides information not only about stability but also about the basins of attraction. Moreover, Liapunov’s method is also applicable to non-autonomous systems .

In its stronger version, Liapunov’s method depends on the construction of a phase function $`F`$ with the properties of being strictly decreasing along the orbits and having a minimum at the equilibrium point. Although stated in most introductory books in dynamical systems, this method is seldom used in the applications because it is difficult in general to find the form of the function $`F`$. However, there are some simple models in physics which indicate an obvious candidate for $`F`$ such as in the case of the classical textbook application, the damped oscillator: the energy of the system is a positive definite form, and decreases along the orbits, thus verifying the condition of applicability of the method in its stronger version.

A wide class of cosmological models share with the latter example the same form and qualitative dynamics. The equations of motion for all these models can be generically written as

$`3{\displaystyle \frac{\dot{a}^2}{a^2}}`$ $`=`$ $`{\displaystyle \underset{i=0}{\overset{N}{}}}{\displaystyle \frac{1}{2}}\dot{b}_i^2+M(a,b_1,\mathrm{},b_N)`$ (1)
$`\ddot{b}_i`$ $`=`$ $`3{\displaystyle \frac{\dot{a}}{a}}\dot{b}_i{\displaystyle \frac{M}{b_i}}(a,b_1,\mathrm{},b_N),i=1,2,\mathrm{},N`$ (2)
where $`a`$ is non-negative.

The degrees of freedom can be divided into two sets. The first corresponds to the global scale factor $`a`$ which has, in the simpler models, a monotone behaviour in time, and the second set is formed by the quantities $`b_i`$, $`i=1,2\mathrm{},N`$, which in many cases behave as damped oscillators and are thus especially well suited to the study and analysis through the construction of Liapunov functions.

The paper is structured as follows. In Section 2, we briefly recall the basic results underlying Liapunov’s method. In Section 3 we shall consider various examples for which a Liapunov function may be constructed that coincides essentially with the energy of the damped oscillator modes. For each case, the relevant dynamical information provided by the method is presented. First, in Section 3.1, we study the basin of attraction of the de Sitter behaviour in Friedman-Robertson-Walker models with an inflaton field having a non-vanishing vacuum. We also assess the conditions for the existence of parametric resonance when the scalar field equation exhibits a periodic forcing term. In Section 3.2 we consider the specific model of single-field reheating  and show that the results derived on the parametric resonance define the admissible class of potentials. Next, in Section 3.3, we consider the convergence of scalar-tensor gravity theories to general relativity in the case of Friedman-Robertson-Walker models . Subsequently, in Section 3.4, we restrict ourselves to the simple case of diagonal Bianchi type-IX cosmologies with no matter and a positive cosmological constant and study the validity of the cosmic no-hair theorem for this model. We shall recover in a simple manner, for this model, the result proved in . Finally, the new information provided in each case by Liapunov’s method is summed up in Section 4.

## 2 Reminder of the basic theorems

Let us start by stating some general definitions.

Given a smooth dynamical system $`\dot{x}=f(x)`$, $`x\mathrm{𝗅𝖱}^n`$, and an equilibrium point $`x_0`$, we say that a continuous function $`F:\mathrm{𝗅𝖱}^n\mathrm{𝗅𝖱}`$, in a neighborhood $`U`$ of $`x_0`$, is a Liapunov function for the point $`x_0`$ if

1. $`F`$ is differentiable in $`U\backslash \{x_0\}`$,
2. $`F(x)>F(x_0)`$ and
3. $`\dot{F}(x)0`$ for every $`xU\backslash \{x_0\}`$.

The existence of a Liapunov function $`F`$ guarantees the stability of $`x_0`$; furthermore, if the strict inequality $`\dot{F}(x)<0`$ holds in $`U\backslash \{x_0\}`$ then $`x_0`$ is asymptotically stable (see for instance )

The neighborhood $`U`$ is essentially constrained by condition (iii). So, when (iii) holds strictly, the method provides information not only about the asymptotically stability of the equilibrium point but also about its basin of attraction, which must contain the set U. This kind of information cannot be obtained neither by linear stability analysis nor by first order perturbation theory.

Again in the case when (iii) holds strictly, the existence of a Liapunov function for the equilibrium point $`x_0`$ has important consequences for the behaviour of time dependent perturbations of the above autonomous system. Consider the perturbed system

$$\dot{x}=f(x)+g(x,t),$$
(3)
where $`g(x,t)`$ is smooth and bounded. Malkin’s theorem  states that the solutions $`x(t;x(t_0),t_0)`$ of (3) will remain for all time inside a prescribed neighborhood $`B_ϵ`$ of $`x_0`$, provided that we pick initial conditions $`x(t_0)`$ sufficiently close to $`x_0`$ and that the magnitude of the perturbation $`g(x,t)`$ is sufficiently small for all $`t>t_0`$ and all $`xB_ϵ`$. This theorem guarantees the stability of $`x_0`$ for any sufficiently small time dependent perturbation, but it does not imply that the solution $`x(t;x(t_0),t_0)`$ tends to $`x_0`$ when $`t\mathrm{}`$, since it does not even require that $`g(x,t)`$ should vanish as $`t\mathrm{}`$. If this is the case, then a related result also due to Malkin (see ) shows that an asymptotic property persists, in the sense that, as time increases, the neighborhood $`B_ϵ`$ may be taken progressively smaller.

Liapunov’s method can also be used in the more general setting of time dependent systems which are not perturbations of autonomous systems (see ). Given a smooth non-autonomous system $`\dot{x}=f(x,t)`$, $`x\mathrm{𝗅𝖱}^n`$, $`t\mathrm{𝗅𝖱}`$, and an equilibrium point $`x_0`$ such that $`f(x_0,t)=0`$, if there exists a smooth function $`F:U\mathrm{𝗅𝖱}^n\times \mathrm{𝗅𝖱}\mathrm{𝗅𝖱}`$, $`U`$ a neighborhood of $`x_0`$, such that $`F(x_0,t)=F_0=const`$, $`F(,t)F_0`$ is positive definite on $`U`$ for every $`t\mathrm{𝗅𝖱}`$ and $`\dot{F}(x,t)0`$ on $`U\times \mathrm{𝗅𝖱}`$ then $`x_0`$ is stable.

This result is analogous to the weak version of the theorem for the autonomous case. In the non-autonomous case, however, the conditions for asymptotic stability are more restrictive. The point $`x_0`$ is asymptotically stable if, moreover, $`F(x,t)`$ is bounded from above, for all $`t`$, by a continuous increasing function of the distance $`xx_0`$, and $`\dot{F}(x,t)`$ is negative definite in $`U`$ for all $`t\mathrm{𝗅𝖱}`$.

## 3 Stability analysis

### 3.1 The basin of attraction of de Sitter inflation

During inflation the leading contribution to the energy momentum tensor is given by the inflaton scalar field $`\varphi `$. The evolution of the Friedman-Robertson-Walker universes is described by

$`3{\displaystyle \frac{\dot{a}^2}{a^2}}+3{\displaystyle \frac{k}{a^2}}`$ $`=`$ $`{\displaystyle \frac{1}{2}}\dot{\varphi }^2+V(\varphi )`$ (4)
$`\ddot{\varphi }`$ $`=`$ $`3{\displaystyle \frac{\dot{a}}{a}}\dot{\varphi }{\displaystyle \frac{\mathrm{d}V}{\mathrm{d}\varphi }}(\varphi )`$ (5)
where $`a`$ is the (positive) scale factor of the universe, and $`k\{1,0,1\}`$ distinguishes the various spatial curvature cases (notice also that we adopt units in which $`8\pi G=1=c`$ throughout the paper). Under some assumptions on the scalar field potential $`V`$ and on the initial conditions $`\varphi _0`$, $`\dot{\varphi }_0`$ and $`a_0`$, there exist solutions where the friction term $`3H\dot{\varphi }`$ dominates in (4) over $`\ddot{\varphi }`$ and the potential term in (5) dominates over the kinetic term . This is the inflationary stage, where the universe evolves with accelerated expansion. For definiteness, we shall consider the simplest model with $`V(\varphi )=\lambda +m^2\varphi ^2/2`$, $`\lambda 0`$.

With the change of variable $`H=\dot{a}/a`$ and $`b=1/a`$ equations (4) and (5) become

$`3H^2`$ $`=`$ $`3kb^2+{\displaystyle \frac{1}{2}}\dot{\varphi }^2+V(\varphi )`$ (6)
$`\dot{b}`$ $`=`$ $`Hb`$ (7)
$`\ddot{\varphi }`$ $`=`$ $`3H\dot{\varphi }{\displaystyle \frac{\mathrm{d}V}{\mathrm{d}\varphi }}(\varphi ).`$ (8)
For future use we shall also need

$`\dot{H}=kb^2{\displaystyle \frac{1}{2}}\dot{\varphi }^2.`$ (9)
Our first goal in this section will be to apply Liapunov’s method for autonomous systems to the characterization of the late time behaviour of the model associated to equations (6), (7) and (8). This system has an equilibrium point at $`(\varphi =0,\dot{\varphi }=0,b=0)`$ corresponding to $`a=\mathrm{}`$ where $`H=\sqrt{\lambda /3}`$. Since $`H=\dot{a}/a`$, this equilibrium point corresponds to an infinite universe with exponential asymptotic expansion rate whenever $`\lambda 0`$. We shall therefore call the point $`(\varphi =0,\dot{\varphi }=0,b=0)`$ the de Sitter spacetime.

Let us prove that for $`k=1,0`$ any initially expanding model will approach the de Sitter spacetime asymptotically for any $`\varphi _0`$, $`\dot{\varphi }_0`$ and $`b_0`$, and that for $`k=1`$ the basin of attraction of the de Sitter attractor is given by $`b_0<\sqrt{\lambda /3}`$.

From (6), for an initially expanding model, $`H`$ can be written as a function of $`b`$, $`\varphi `$ and $`\dot{\varphi }`$ as

$`H(\varphi ,\dot{\varphi },b)={\displaystyle \frac{1}{\sqrt{3}}}\left(3kb^2+{\displaystyle \frac{1}{2}}\dot{\varphi }^2+\lambda +{\displaystyle \frac{1}{2}}m^2\varphi ^2\right)^{1/2}.`$ (10)
Notice that the argument of the square root is always positive for $`k`$ non-positive. For the case $`k=1`$ this is also true if the initial value of the scale factor satisfies $`b_0<\sqrt{\lambda /3}`$, because then $`b`$ is a monotonically decreasing function of time.

The asymptotic stability of the de Sitter attractor is obtained, as in the case of the classic damped oscillator, by using a Liapunov function given by

$`F(b,\varphi ,\dot{\varphi })={\displaystyle \frac{3}{2}}b^2+{\displaystyle \frac{1}{2}}\dot{\varphi }^2+V(\varphi )>0.`$ (11)
From equations (7), (8) and (9) we have

$`\dot{F}(b,\varphi ,\dot{\varphi })=3H\left(\dot{\varphi }^2+b^2\right)<0,`$ (12)
since $`H(\varphi ,\dot{\varphi },b)>0`$. Then, $`F`$ given by (11) satisfies (i), (ii) and (iii) for equations (7), (8), and (iii) holds strictly. This means that for $`k=1,0`$ and $`\lambda 0`$ the de Sitter spacetime is a global attractor, while for $`k=1`$ and $`\lambda >0`$ its basin of attraction contains all initial conditions $`(\varphi _0,\dot{\varphi }_0,b_0)`$ such that $`b_0<\sqrt{\lambda /3}`$.

Consider now the following modification of the equations of motion

$`3H^2`$ $`=`$ $`3kb^2+{\displaystyle \frac{1}{2}}\dot{\varphi }^2+V(\varphi )`$ (13)
$`\dot{b}`$ $`=`$ $`Hb`$ (14)
$`\ddot{\varphi }`$ $`=`$ $`3H\dot{\varphi }m^2\left(1+\epsilon f(t)\right)\varphi ,`$ (15)
where $`f`$ is a bounded, continuous, periodic function with period close to the free oscillating period of the unperturbated equation (15).

We shall take equations (13), (14), (15), as a toy model for the study of the problem of reheating, or parametric resonance, and defer to the following subsection its application to a specific model.

Because (iii) holds strictly Malkin’s theorem can be applied to this case. It states that the solutions of (13), (14) and (15), will remain for all time inside a prescribed neighborhood of the infinity attractor provided that we pick initial conditions sufficiently close to it and that the magnitude of the perturbation is sufficiently small for all $`t`$. This theorem guarantees the stability of the infinity attractor for any sufficiently small time-dependent perturbation, but it does not imply that the solution tends to it as time goes to infinity. Nevertheless, we can look for some hypothesis concerning $`H`$ which will entail the asymptotic stability of the equilibrium point $`\varphi =\dot{\varphi }=b=0`$. Let

$`F(b,\varphi ,\dot{\varphi })={\displaystyle \frac{1}{2}}\dot{\varphi }^2+{\displaystyle \frac{1}{2}}\left(3\alpha H(\varphi ,\dot{\varphi },b)+m^2\right)\varphi ^2+\alpha \varphi \dot{\varphi }+{\displaystyle \frac{3}{2}}b^2,`$ (16)
which is positive definite for $`0<\alpha <m`$. We obtain, using (13), (14), (15),

$`\dot{F}=\left(3H\alpha \right)\dot{\varphi }^2`$ $``$ $`\alpha \left(m^2{\displaystyle \frac{3}{2}}\dot{H}+m^2\epsilon f(t)\right)\varphi ^2m^2\epsilon f(t)\varphi \dot{\varphi }`$ (17)
$``$ $`3Hb^2.`$
$`\dot{F}`$ will be negative definite if

$`\left(3H\alpha \right)\left(m^2{\displaystyle \frac{3}{2}}\dot{H}+m^2\epsilon f(t)\right){\displaystyle \frac{1}{4}}m^4\epsilon ^2f^2(t)>0.`$ (18)
Since $`\dot{H}(\varphi =0,\dot{\varphi }=0,b=0)=0`$, this condition will hold for $`t`$ large enough and $`ϵ`$ small enough whenever $`\underset{t+\mathrm{}}{lim}H(t)>\alpha /3`$.

There are two basic conclusions that can be obtained from these facts. First, parametric resonance can only occur when (18) does not hold, that is, parametric resonance can only occur if $`H`$ tends to zero as time goes to infinity, because, otherwise, the Liapunov function (16), ensures that the equilibrium point $`\varphi =\dot{\varphi }=b_0=0`$ is asymptotically stable. This implies, for this toy model, that reheating can only be obtained from a potential that has a null minimum $`V(0)=\lambda =0`$. Second, even if $`H`$ tends to zero as time goes to infinity, parametric resonance can only be observed if $`\epsilon `$ is greater than some critical value because otherwise Malkin’s theorem implies the stability of the infinity attractor.

### 3.2 Single field reheating

We now examine a specific model of reheating after inflation. According to the inflationary theory, almost all elementary particles populating the universe were created during the process of reheating after inflation . Reheating occurs due to particle production by the oscillating scalar field $`\varphi `$. In the simplest models, this field $`\varphi `$ is the same inflaton field that drives inflation at the early stages of the evolution of the universe. After inflation, the inflaton oscillates near the minimum of its potential and this induces the quantum process of creation of particles . At the root of this process is a phenomenon of parametric resonance which gives rise to the excitation of certain modes of the quantum fluctuations of fields .

We shall consider a simple model with the scalar field potential of the form $`V(\varphi )=V_0\varphi ^n`$, with $`V_0>0`$, characteristic of the chaotic inflation scenario, and of the original studies of reheating. Our purpose is to ilustrate the possibilities underlying the aplication of the results of the previous subsection to this model. We aim at determining under what conditions does the mechanism of resonance occur, avoiding the drastic assumption often made in the analytic treatment found in the literature which consists in neglecting the Hubble rate of expansion . Furthermore, we also wish to show that the theorem on parametric resonance requiring the damping term to be vanishing in time and the amplitude of the periodic perturbation in the scalar field equation to be bounded from below permits to constrain the form of the potential in the case under consideration.

We consider a FRW metric with non-positive curvature and shall focus on perturbations of the single scalar field ignoring the accompanying metric fluctuations. This is well justified on scales smaller than the Hubble distance, that is, satisfying $`paH`$, where $`p`$ is the wave number of the perturbation, since relativistic effects play a lesser role on these scales . On scales larger than the Hubble radius this approximation is more restrictive and, in particular, it has been recently shown to miss the resonant behaviour of the metric perturbations . With this caveat we shall thus consider the set of equations of motion

$`3H^2`$ $`=`$ $`3kb^2+{\displaystyle \frac{1}{2}}\dot{\varphi }^2+V(\varphi )`$ (19)
$`\dot{b}`$ $`=`$ $`Hb`$ (20)
$`\ddot{\varphi }`$ $`=`$ $`3H\dot{\varphi }nV_0\varphi ^{n1}+b^2_\stackrel{}{x}^2\varphi ,`$ (21)
with $`n>1`$ and $`\stackrel{}{x}`$ the spatial coordinates.

Following , we decompose the inflaton field into the sum of a classical inflaton field $`\phi `$ and a quantum fluctuation $`\delta \varphi `$ in the form $`\varphi (\stackrel{}{x},t)=\phi (t)+\delta \varphi (\stackrel{}{x},t)`$. We next concentrate in a specific mode of the fluctuation $`\delta \varphi `$. This yields, keeping the first order terms in $`\delta \varphi _p`$,

$`\ddot{\phi }+3H\dot{\phi }+nV_0\phi ^{n1}=0`$ (22)
$`\delta \ddot{\varphi }_p+3H\delta \dot{\varphi }_p+\left[n(n1)V_0\phi ^{n2}+p^2b^2\right]\delta \varphi _p=0.`$ (23)
Notice that $`n`$ should be strictly greater than two for the existence of resonant behaviour.

Let $`d\tau /dt=b`$ and denote the derivative $`d/d\tau `$ by a prime. The equations of motion become

$`3H^2`$ $`=`$ $`3kb^2+{\displaystyle \frac{1}{2}}b^2\phi ^2+V(\phi )`$ (24)
$`b^{}`$ $`=`$ $`H`$ (25)
$`\phi ^{\prime \prime }`$ $`=`$ $`2{\displaystyle \frac{b^{}}{b}}\phi ^{}nb^2V_0\phi ^{n1}`$ (26)
$`\delta \varphi _p^{\prime \prime }`$ $`=`$ $`2{\displaystyle \frac{b^{}}{b}}\delta \varphi _p^{}\left[n(n1)b^2V_0\phi ^{n2}+p^2\right]\delta \varphi _p`$ (27)
with $`n>2`$.

Let $`T=\dot{\phi }^2/2=b^2\phi ^2/2`$ and $`\gamma =2T/(T+V)`$. Since $`2\overline{T}=n\overline{V}`$ holds asymptotically , we have $`\overline{\gamma }=2n/(n+2)`$ and a simple argument shows that $`Hb^{3\overline{\gamma }/2}`$ if $`k=0`$, $`Hb`$ if $`k=1`$. Consequently $`b^{}/b`$ and $`b^2\phi ^{n2}`$ approach the infinity attractor as

$`b^{}/b`$ $``$ $`\{\begin{array}{cc}b^{(2n2)/(n+2)}\hfill & \text{if }k=0\hfill \\ b\hfill & \text{if }k=1\hfill \end{array}`$ (28)
$`\phi ^{n2}b^2`$ $``$ $`b^{(4n16)/(n+2)}.`$ (29)
From the previous discussion of our toy model we know that, in order to have parametric resonance, the damping term in (27) must tend to zero and the perturbing term in (27) must be bounded away from zero. The first condition is always satisfied, in view of (28) and of the fact that $`n`$ must be greater than $`2`$. From (29), the second condition is satisfied provided that $`n4`$. Notice that for $`n=4`$ the term $`\phi ^{n2}b^2`$ in (27) oscillates with constant amplitude for all time. This is the case usually studied in the literature (see for example ). However, (29) shows that reheating should be more effective for potentials with exponents closer to 2 in the allowed interval $`(2,4]`$.

### 3.3 General Relativity as a Cosmological attractor <br>of Scalar-Tensor Gravity Theories

General scalar-tensor gravity theories are based on the Lagrangian

$`L_\mathrm{\Phi }=\mathrm{\Phi }R{\displaystyle \frac{\omega (\mathrm{\Phi })}{\mathrm{\Phi }}}g^{\mu \nu }\mathrm{\Phi }_{,\mu }\mathrm{\Phi }_{,\nu }+2U(\mathrm{\Phi })+16\pi L_m`$ (30)
where $`R`$ is the Ricci curvature scalar of a space time endowed with a metric $`g_{\mu \nu }`$, $`\mathrm{\Phi }`$ is a scalar field, $`\omega (\mathrm{\Phi })`$ is a dimensionless coupling function, $`U(\mathrm{\Phi })`$ is a function of $`\mathrm{\Phi }`$, and $`L_m`$ is the Lagrangian for the matter fields. They provide a most natural generalization of Einstein’s general relativity, and their investigation enables a generic model-independent approach to the main features and cosmological implications of the unifications schemes. The archetypal of these theories is Brans-Dicke theory  for which $`\omega (\mathrm{\Phi })`$ is a constant and the most distinctive feature of the theories defined by (30) is the variation of the gravitational constant, since $`G=\mathrm{\Phi }^1`$.

Consider the homogeneous and isotropic Friedman-Robertson-Walker (FRW) universes and assume that the matter content of the universe is a perfect fluid with the usual equation of state $`p=(\gamma 1)\rho `$, where $`0\gamma 2`$. Here $`p`$ and $`\rho `$ are, respectively, the pressure and the energy-density measured by a comoving observer. Performing the conformal transformation of the metric  $`\stackrel{~}{g}_{ab}=(\mathrm{\Phi }/\mathrm{\Phi }_{})g_{ab}`$, where $`\mathrm{\Phi }_{}`$ is a constant which normalizes the gravitational constant, to the so-called Einstein frame which corresponds to replacing the dynamical variables $`(a,t,\mathrm{\Phi })`$ with $`(\sqrt{\frac{\mathrm{\Phi }}{\mathrm{\Phi }_{}}}a,\sqrt{\frac{\mathrm{\Phi }}{\mathrm{\Phi }_{}}}t,\phi )`$ where

$`{\displaystyle \frac{d\mathrm{ln}\mathrm{\Phi }}{d\phi }}=\sqrt{{\displaystyle \frac{16\pi }{\mathrm{\Phi }_{}(2\omega (\phi )+3)}}},`$ (31)
the field equations reduce to the simple form

$`3{\displaystyle \frac{\dot{a}^2}{a^2}}+3{\displaystyle \frac{k}{a^2}}`$ $`=`$ $`{\displaystyle \frac{1}{2}}\dot{\phi }^2+a^{3\gamma }M(\phi )+V(\phi )`$ (32)
$`\ddot{\phi }+3{\displaystyle \frac{\dot{a}}{a}}\dot{\phi }`$ $`=`$ $`a^{3\gamma }M^{}(\phi ){\displaystyle \frac{\mathrm{d}V}{\mathrm{d}\phi }}(\phi ),`$ (33)
where $`V(\phi )=\mathrm{\Phi }_{}^2U\left(\mathrm{\Phi }(\phi )\right)/\left(8\pi \mathrm{\Phi }^2(\phi )\right)`$, $`M(\phi )=\mu \left(\mathrm{\Phi }(\phi )/\mathrm{\Phi }_{}\right)^{2+3\gamma /2}`$ and $`\mu `$ is a constant which fixes the initial conditions for the energy-density of the perfect fluid. We shall assume that $`M(\phi )0`$ (which corresponds to having $`G0`$) and that $`V(\phi )=\lambda +\frac{1}{2}m^2\phi ^2`$, i.e., that $`V(\phi )`$ has a non-degenerate positive global minimum at the origin.

Let $`H=\dot{a}/a`$ and $`b=1/a`$. It is easy to see that the former system is equivalent to

$`3H^2`$ $`=`$ $`3kb^2+{\displaystyle \frac{1}{2}}\dot{\phi }^2+b^{3\gamma }M(\phi )+V(\phi )`$ (34)
$`\dot{b}`$ $`=`$ $`Hb`$ (35)
$`\ddot{\phi }`$ $`=`$ $`3H\dot{\phi }b^{3\gamma }M^{}(\phi ){\displaystyle \frac{\mathrm{d}V}{\mathrm{d}\phi }}(\phi ).`$ (36)
These equations define a dynamical system of dimension $`3`$. We are interested in studying the stability of the equilibrium point at $`a=+\mathrm{}`$ corresponding to $`b=0`$, $`H=\sqrt{\lambda /3}`$ and $`\phi =\dot{\phi }=0`$, which is the de Sitter solution .

Let us prove that for $`k=1,0`$ any initially expanding model will approach the de Sitter spacetime asymptotically for any $`\phi _0`$, $`\dot{\phi }_0`$ and $`b_0`$, and that for $`k=1`$ the basin of attraction of the de Sitter attractor is given by $`b_0<\sqrt{\lambda /3}`$. Consider the initial conditions $`H_0>0`$ associated with an expanding model, and $`b_0`$, $`\phi _0`$ and $`\dot{\phi }_0`$ arbitrary. Using the hypothesis $`M(\phi )0`$, the expanding branch $`H=H(\phi ,\dot{\phi },b)`$ of (35)is given by

$`H={\displaystyle \frac{1}{\sqrt{3}}}\left\{3kb^2+{\displaystyle \frac{1}{2}}\dot{\phi }^2+b^{3\gamma }M(\phi )+V(\phi )\right\}^{1/2}>0.`$ (37)
As for the case $`k=1`$, we take again $`H_0>0`$ and $`b_0\sqrt{\lambda /3}`$ and, since, from (35), $`b`$ is a strictly decreasing function of time, it is possible for these initial conditions to write $`H=H(\phi ,\dot{\phi },b)`$ as in (37).

The function defined by

$`F(b,\phi ,\dot{\phi })={\displaystyle \frac{1}{2}}\dot{\phi }^2+\lambda +{\displaystyle \frac{1}{2}}m^2\phi ^2+b^{3\gamma }M(\phi )`$ (38)
is a Liapunov function for the equilibrium point $`(b,\phi ,\dot{\phi })=(0,0,0)`$. Indeed using (34), (35), (36) and (37), it is easy to see that

$`\dot{F}(b,\phi ,\dot{\phi })=3H\left(\dot{\phi }^2+\gamma b^{3\gamma }M(\phi )\right)<0,`$ (39)
and it is also straightforward to check that (i) and (ii) hold for $`F`$. Thus this proves that the equilibrium point $`(b,\phi ,\dot{\phi })=(0,0,0)`$ is asymptotically stable and shows that an attractor mechanism indeed exists when the matter content of the universe is a perfect fluid and the potential has a non-degenerate non-negative minimum. This means that for $`k=1,0`$ and $`\lambda 0`$ the de Sitter spacetime is a global attractor, while for $`k=1`$ and $`\lambda >0`$ its basin of attraction contains all initial conditions $`(\phi _0,\dot{\phi }_0,b_0)`$ such that $`b_0<\sqrt{\lambda /3}`$. A similar attracting mechanism should also exist in the non-isotropic case.

In general, linear stability analysis only provides information about the stability and asymptotic behaviour in a undetermined neighborhood of the equilibrium point. This can be easily seen in the following example. Consider the equations of motion (35), (36) for $`\gamma =4/3`$ (for $`\gamma 4/3`$ the linear stability method can not be applied without additional hypothesis on $`M(\phi )`$). The eigenvalues associated with the de Sitter attractor are $`\lambda _b=\sqrt{\lambda /3}`$, $`\lambda _\phi =\sqrt{3\lambda }\sqrt{3\lambda 4m^2}`$ and $`\lambda _{\dot{\phi }}=\sqrt{3\lambda }+\sqrt{3\lambda 4m^2}`$. The linear approximation confirms the result obtained by the Liapunov method, the real part of the eigenvalues is always negative, and shows that there exist two basic behaviours for the field associated to the monotone approach of the scale factor to the de Sitter attractor: an oscillatory mode (damped oscillations) if $`0<\lambda <4m^2/3`$ and a monotone behaviour (over-damped) for the field if $`\lambda 4m^2/3`$ which corresponds, respectively, to a stable focus and a stable node in the plane $`(\phi ,\dot{\phi },b=0)`$. No information is obtained concerning the size of the basin of attraction for the de Sitter attractor with this method.

Also notice that when $`\gamma <1/3`$ the second member of equation (36) is not differentiable at $`b=0`$, which precludes the application of linear stability analysis in its simplest form. This fact emphasizes the advantage of the use of the Liapunov’s method.

### 3.4 Asymptotic behaviour for diagonal Bianchi type-IX cosmologies in the presence of a positive cosmological constant without matter.

By definition, the general Bianchi type-IX spacetime has topology $`\mathrm{𝗅𝖱}\times 𝐒^3`$, with a simply transitive action of the isometry group $`SU(2)`$ on $`𝐒^3`$ spatial slices. The metric of a general Bianchi type-IX model can be put in the form

$`ds^2=dt^2+a^2(t){\displaystyle \underset{i,j=1}{\overset{3}{}}}\left[e^\beta \right]_{ij}d\omega ^id\omega ^j,`$ (40)
where $`\omega ^i`$, $`i=1,2,3`$, are isometry invariant one-forms on the three-sphere, $`a`$ is a scalar, and $`\beta `$ is a traceless $`3\times 3`$ matrix. Both $`a`$ and $`\beta `$ are functions of the proper time $`t`$ only.

For a diagonal spacetime, let $`\beta _i`$, $`i=1,2,3`$, denote the diagonal elements of the matrix $`\beta `$. Only two of these quantities are independent since $`\beta _1+\beta _2+\beta _3=0`$ on account of the tracelessness of $`\beta `$. We choose the independent variables to be

$`b_1`$ $`=`$ $`{\displaystyle \frac{1}{2\sqrt{6}}}\beta _3`$ (41)
$`b_2`$ $`=`$ $`{\displaystyle \frac{1}{6\sqrt{2}}}\left(\beta _1\beta _2\right).`$ (42)
The vacuum Einstein equations with the cosmological constant take the form

$`3{\displaystyle \frac{\dot{a}^2}{a^2}}={\displaystyle \frac{1}{2}}\left(\dot{b}_1^2+\dot{b}_2^2\right)+{\displaystyle \frac{1}{a^2}}V(b_1,b_2)+\mathrm{\Lambda }`$ (43)
$`\ddot{b}_i+3{\displaystyle \frac{\dot{a}}{a}}\dot{b}_i+{\displaystyle \frac{1}{a^2}}{\displaystyle \frac{V}{b_i}}=0,i=1,2`$ (44)
with

$`V(b_1,b_2)`$ $`=`$ $`e^{\sqrt{\frac{2}{3}}b_1}\mathrm{cosh}(\sqrt{2}b_2)+{\displaystyle \frac{1}{4}}e^{4\sqrt{\frac{2}{3}}b_1}`$
$`+{\displaystyle \frac{1}{2}}e^{2\sqrt{\frac{2}{3}}b_1}\left[\mathrm{cosh}(2\sqrt{2}b_2)1\right]`$
Let $`H=\dot{a}/a`$ and $`b=1/a`$. It is easy to see that the former system is equivalent to

$`3H^2`$ $`=`$ $`{\displaystyle \frac{1}{2}}\left(\dot{b}_1^2+\dot{b}_2^2\right)+\mathrm{\Lambda }+b^2V(b_1,b_2)=0`$ (46)
$`\dot{b}`$ $`=`$ $`Hb`$ (47)
$`\ddot{b}_i+3H\dot{b}_i+b^2{\displaystyle \frac{V}{b_i}}=0,i=1,2.`$ (48)
These equations define a dynamical system of dimension $`5`$. This system has two finite equilibrium points $`b=0`$, $`H=\sqrt{\mathrm{\Lambda }/3}`$ and $`b=2\sqrt{\mathrm{\Lambda }/3}`$, $`H=0`$ both with $`b_1=b_2=\dot{b}_1=\dot{b}_2=0`$. We are interested in the stability of the equilibrium point $`b=0`$ associated with the de Sitter attractor.

Notice that near the minimum $`(b_1,b_2)=(0,0)`$ the potential (3.4) takes the form

$`V(b_1,b_2)`$ $`=`$ $`{\displaystyle \frac{3}{4}}+\left(b_1^2+b_2^2\right)+O_3(b_1,b_2)`$ (49)
$`=`$ $`{\displaystyle \frac{3}{4}}+v(b_1,b_2),`$
where $`v(0,0)=0`$.

We take $`H_0>0`$ and $`b_02\sqrt{\mathrm{\Lambda }/3}`$. Since, from (47), $`b`$ is a strictly monotone decreasing function with time the expanding branch of $`H=H(b,b_1,b_2,\dot{b}_1,\dot{b}_2)`$ can be written for these initial conditions as

$`H={\displaystyle \frac{1}{\sqrt{3}}}\left\{{\displaystyle \frac{1}{2}}\left(\dot{b}_1^2+\dot{b}_2^2\right)+b^2v(b_1,b_2)+\mathrm{\Lambda }{\displaystyle \frac{3}{4}}b^2\right\}^{1/2}>0.`$ (50)
The Liapunov function for the equilibrium point $`(b,b_1,b_2,\dot{b}_1,\dot{b}_2)=(0,0,0,0,0)`$ will now be

$`F(b,b_1,b_2,\dot{b}_1,\dot{b}_2)={\displaystyle \frac{1}{2}}\left(\dot{b}_1^2+\dot{b}_2^2\right)+b^2v(b_1,b_2)+{\displaystyle \frac{3}{2}}b^20.`$ (51)
Indeed using the equations of motion we have

$`\dot{F}=3H\left[\dot{b}_1^2+\dot{b}_2^2+{\displaystyle \frac{2}{3}}b^2v(b_1,b_2)+b^2\right]<0`$ (52)
and thus it is a simple matter to check that all the conditions (i), (ii) and (iii) are satisfied. Therefore the equilibrium point $`(b,b_1,b_2,\dot{b}_1,\dot{b}_2)=(0,0,0,0,0)`$ is asymptotically stable and the basin of attraction for de Sitter spacetime contains all initial conditions $`(b_0,b_1(0),b_2(0),\dot{b}_1(0),\dot{b}_2(0))`$ such that $`b_0<2\sqrt{\mathrm{\Lambda }/3}`$. In this case the linear stability analysis reveals only that the eigenvalues of the de Sitter attractor are all negative.

## 4 Summary and conclusions

In this paper we have applied the Liapunov’s method to the study of the general asymptotic behaviour of cosmological models, namely: FRW cosmologies with inflaton field; scalar-tensor gravity theories in FRW; diagonal Bianchi-IX vacuum cosmologies in the presence of a positive cosmological constant. For these models we recover through this method the well known stability result for the de Sitter solution and we show that the method also provides information about its basin of attraction.

Furthermore, we have applied Liapunov’s method to another problem of a completely different nature, namely the analysis of the conditions for the existence of the parametric resonance mechanism underlying the reheating process. We have concluded that in order to have parametric resonance, the damping term in the quantum perturbation equation must tend to zero and the forcing term must be bounded away from zero. This implies that for the basic model with potentials of the form $`V(\phi )=V_0\phi ^n`$, reheating may occur only when $`n(2,4]`$. As far as we know, albeit the simplicity of the reheating model considered, this conclusion is new and fully illustrates the extension to which Liapunov’s method may be utilized to draw rigorous results in cosmology.

## Acknowledgements

The authors are grateful to Alfredo B. Henriques for useful discussions and wish to acknowledge the finantial support from Fundação de Ciência e Tecnologia under the grant PBIC/C/FIS/2215/95.