text:
# Glass-like dynamical behavior in hierarchical models submitted to continuous cooling and heating processes

## I Introduction

In recent years, there has been quite a considerable amount of work in models in which glassy behavior is generated not by quenched disorder, but by kinetic constraints. The kinetic restrictions are responsible for the slow relaxation, since the state of a particle or a group of particles can only change if some condition of its environment is fulfilled. In particular, “facilitated” models have been considered, both for structural glasses  and for granular systems . The characteristic feature of facilitated models is that a particle (spin) can only change its state if a certain number of its neighbors is in an excited state. Also, hierarchically constrained models have been used to study stretched exponential relaxation in glasses . In these models, the system is structured in levels and a particle in a given level can only make a transition if a given cluster of particles in the lower level is in a certain subset of configurations. Then, the dynamics of the several levels are coupled, and the characteristic relaxation times increase with the level index. Hierarchically constrained dynamics may be relevant for those complex systems in which the time evolution of the slowest modes is controlled by the relaxation of the fastest ones. This qualitative picture is adequate to describe, among other problems, proteins relaxation  and the densification of powders and structural glasses at high pressure . Very recently, a kind of hierarchically constrained dynamics has been shown to exhibit, in quite a natural way, logarithmic relaxation . This kind of “anomalous”, highly nonexponential, decay is observed in a wide variety of complex systems, including spin-glasses , granular materials , structural glasses  and protein models .

The aim of this work is to study the dynamical behavior of a general class of hierarchically constrained models when submitted to more complicated processes. In particular, we are interested in the behavior of a system with hierarchical constraints when the temperature changes in time, what makes the coupling between the levels be time-dependent. The consideration of time-dependent temperature processes requires an extension of the original model as formulated by Palmer et al. . This will be done in a very simple, but natural, way: the coupling between the levels vary in time because the probability of a cluster configuration allowing to relax a particle in the next level depends on the temperature. On the other hand, neither the number of particles in a given level nor the length of the clusters “facilitating” the relaxation depend on the temperature. They are considered as quantities defined in the coarse-grained description of the system introduced to model the physical problem at hand.

Let us note that hierarchical models can also be applied to the analysis of non-thermal systems, as granular materials in the dense regime. For those materials, thermal energy is not enough to make the system explore the phase space of configurations. Then, the system must be externally excited, for instance vibrating it, in order to be able to evolve. In these situations, the role of the temperature is played by the intensity of the external driving. If the stationary state reached by the system in the long time limit can be described by Edward’s theory , the compactivity $`X`$, which is the analogous of the temperature in thermal systems, will be a function of the intensity of the external force. Then, by exploiting the analogies of Edward’s theory, i. e., substituting volume by energy and compactivity by temperature, it is possible to incorporate non-thermal systems in our formulation.

Processes in which the temperature is time-dependent are physically relevant because they can be used to study some characteristic dynamical aspects of glasses. For instance, when a supercooled liquid is cooled down to very low temperatures, a laboratory glass transition is observed. A dramatic change in the behavior of the system takes place, and it departs from the equilibrium curve, getting “frozen” in a far from equilibrium state. This transition appears as a consequence of the fast increase of the relaxation time with decreasing temperature. In order to characterize the cooling process, experimental physicists often use the residual value of the relevant physical properties, i. e., the difference between their actual values over the cooling curve and the value obtained by extrapolation of the equilibrium curve to the very low temperature region . If the system is reheated from the nonequilibrium state, it returns to equilibrium for high temperatures, but it follows a different curve from the cooling one, giving rise to hysteresis effects. This phenomenon is related to the “nonlinearity” of glassy relaxation: the approach towards the equilibrium curve depends on the configuration of the system, measured by the so-called fictive temperature . Narayanaswami’s theory provides a phenomenological explanation of this behavior . Interestingly, a similar behavior has been found in vibrated granular materials when the tapping intensity is varied in a cyclic way , although the hysteresis effects are more evident when the heating process begins in a loosely packed state, referred as to the “irreversible” branch in the experiments.

We will start from a very general hierarchical spin model, in which pseudospins are organized into levels, labelled by an index $`n`$. The pseudospins are assumed to correspond to some coarse-grained description of the system. They can take only two values, representing, for instance, two possible densities of a certain small subvolume of the system. One pseudospin in level $`n+1`$ can only flip between the two possible values if a cluster of $`\mu _n`$ spins in level $`n`$ is in a given subset of configurations. This is the basic characteristic of hierarchically constrained models as introduced in Ref. , and it slows down the relaxation in level $`n+1`$, as compared with that of level $`n`$. Here we will consider the simple choice that all the spins in the cluster must be in the up (excited) state. The exact dynamical equation for the evolution of the pseudospins involves very complicated moments of the probability distribution. To get an exactly solvable model, a “mean-field” approximation will be introduced. Then, the characteristic relaxation time $`\tau _n`$ of level $`n`$ is seen to increase both with the index label $`n`$, due to the hierarchical constraint, and also with decreasing temperature, since the configurations allowing the system to relax become less probable when the temperature is lowered.

Some exact dynamical results for systems described with master equations with time-dependent transition rates are known. In particular, the existence of a “normal” solution, i. e., a solution of the master equation that is approached by all the others, has been proved on a very general basis . The main required conditions are the irreducibility of the Markov process for long enough times and that the transition rates be externally controlled, so that they do not depend on the probability distribution of the system. Moreover, it has been established that the normal solution tends to the equilibrium curve for very high temperatures, in continuous heating processes . The linear correction of the normal solution with respect to the equilibrium curve has been computed using Hilbert’s method . It is not evident whether these results still hold when approximations are introduced in the dynamics of the system. In particular, in “mean-field” type approximations, the demonstrations in Refs.  are not valid, since the transition rates become functionals of the probability distribution function. Nevertheless, we will show that all the above mentioned properties apply to our simplified model. This is a good test of the plausibility of the approximations carried out, and perhaps an indication that the results obtained here are more general than the derivations in Refs. .

The organization of the paper is as follows. In Sec. II the hierarchical model is introduced, and the exact evolution equation for the average spin is obtained. By introducing a mean-field approximation, this equation can be closed. Afterwards, a specific, but quite general, choice for the functions defining the model is made. This allows to introduce a continuous limit in which the relaxation of the system at constant temperature is solved in Sec. III. Time-dependent temperature processes are considered in Sec. IV, where the general solution for the evolution of the probability distribution and the average energy are obtained. The general solution is similar to the expression proposed by Narayanaswami on a phenomenological basis . Section IV A is devoted to the analysis of Hilbert’s expansion, which is valid in the very high temperatures regime. Cooling processes are addressed in Sec. V where, for the sake of simplicity, a concrete cooling law is studied, for which the residual properties are analytically calculated. A qualitative analysis of the glass-like transition is presented in Sec. V A. It allows to give very good estimates of the residual properties, and leads to the introduction of the concept of fictive temperature in a very natural way. The behavior of the system when it is reheated from low temperatures is considered in Sec. VI. The main role played by the normal solution for the understanding of the hysteresis effects shows up. Moreover, the analysis clearly indicates that the relevance of the normal solution is not restricted to near equilibrium situations. Finally, a discussion of the main points in this work is given in Sec. VII.

## II Dynamics of hierarchically constrained models

In this section a general kind of spin models with hierarchically constrained dynamics will be introduced. We will focus on the evolution of the average value of the spin, which is supposed to be the relevant variable. For instance, in a thermal system it will be directly related to the mean energy. Then, let us consider a system whose degrees of freedom can be classified into levels, labelled by an index $`n=0,1,2,\mathrm{},n_{\text{m}ax}`$. The degrees of freedom in level $`n`$ will be represented by $`N_n`$ pseudospins, $`\sigma _i^{(n)}=\pm 1`$, $`i=1,2,\mathrm{},N_n`$. The hamiltonian of the system is assumed to have the form

$$=h\underset{n=0}{\overset{n_{\text{m}ax}}{}}\underset{i=1}{\overset{N_n}{}}m_i^{(n)},m_i^{(n)}=\frac{1+\sigma _i^{(n)}}{2}.$$
(1)
Note that $`m_i^{(n)}`$ is the occupation number of the “up” (+1) state of the corresponding site. In order to write Eq. (1) we have supposed that there is no interaction between the pseudospins, but there is an “external field” $`h`$. For a thermal system, $``$ gives the energy of a given microstate of the system, while for a non-thermal system, like a powder, it could be interpreted as the volume of a given, mechanically stable, configuration of “grains”. Using the terminology for thermal systems, the average value of the dimensionless energy per spin over the ensemble of systems considered is

$$\epsilon =\frac{}{Nh}=\frac{1}{N}\underset{n=0}{\overset{n_{\mathrm{max}}}{}}\underset{i=1}{\overset{N_n}{}}p_i^{(n)},$$
(2)
where

$$p_i^{(n)}=m_i^{(n)}=\frac{1+\sigma _i^{(n)}}{2}$$
(3)
is the probability that the $`i`$-th spin of level be in the up state, $`N=_{n=0}^{n_{\mathrm{max}}}N_n`$ is the total number of pseudospins, and the angular brackets denote statistical ensemble average.

The system is considered to be in contact with a heat bath at temperature $`T`$, so that the equilibrium average value of the pseudospins does not depend either on $`i`$ nor on $`n`$, and it is given by

$$\sigma _e\sigma _i^{(n)}_e=\mathrm{tanh}\left(\frac{1}{T^{}}\right).$$
(4)
Here $`T^{}`$ is a dimensionless temperature, $`T^{}=2k_BT/h`$, $`k_B`$ being Boltzmann’s constant. For the sake of concision, we will drop the asterisk in the following. The above average value of the spin follows from the equilibrium probability for the “up” state of any spin,

$$p_ep_{i,e}^{(n)}=\frac{e^{\frac{1}{T}}}{e^{\frac{1}{T}}+e^{\frac{1}{T}}}.$$
(5)
In the limit of infinite temperature, or zero external field, both states of the pseudospins are equiprobable, $`p_e=1/2`$, and $`\sigma _e=0`$. From Eq. (2), it follows that the equilibrium value of $`\epsilon `$ is

$$\epsilon _e=p_e.$$
(6)
For granular materials, the role of the temperature $`T`$ is played by the compactivity , which is linked to the intensity of the perturbation allowing the system to explore the configuration space.

The dynamics of the model is formulated by means of a master equation with single-spin flip Glauber transition rates . Let us consider the flip of a given spin $`\sigma _i^{(n)}`$. This transition connects a given configuration $`𝝈`$ of the whole system with the configuration $`R_i^{(n)}𝝈`$, where $`R_i^{(n)}`$ is the operator which rotates the spin $`\sigma _i^{(n)}`$, keeping all the other spins the same. The transition rate for the flip of the spin $`\sigma _i^{(n)}`$ in configuration $`𝝈`$ is

$$W_i^{(n)}(𝝈)=\frac{1}{2}\alpha _i^{(n)}(𝝈)\left[1+\sigma _i^{(n)}\mathrm{tanh}\left(\frac{1}{T}\right)\right].$$
(7)
The characteristic relaxation rate $`\alpha _i^{(n)}`$ of the spin $`\sigma _i^{(n)}`$ depends on the configuration $`𝝈`$ of the system through the hierarchical constraint

$$\alpha _i^{(n)}(𝝈)=\left[\alpha _{k_i}^{(n1)}(𝝈)\alpha _{k_i+1}^{(n1)}(𝝈)\mathrm{}\alpha _{k_i+\mu _{n1}1}^{(n1)}(𝝈)\right]^{1/\mu _{n1}}\underset{j=k_i}{\overset{k_i+\mu _{n1}1}{}}\delta _{\sigma _j^{(n1)},+1},$$
(8)
where $`\delta _{ij}`$ is the Kronecker delta. This expression implies that the spin $`\sigma _i^{(n)}`$ needs, in order to flip, that all the spins belonging to a cluster of length $`\mu _{n1}`$ starting at a given spin $`k_i`$ of level $`n1`$ must be in the up (+1) state. Besides, the characteristic flip rate of the spin is the average of the characteristic flip rates of the spins belonging to the cluster determining its possibility of change. This restricting condition is schematically depicted in Fig. 1. Note that the possibility of a given spin in level $`n`$ to flip is restricted by the state of a set of clusters in all levels $`n^{}<n`$, the number of clusters involved in each level increasing as $`n^{}`$ decreases. The hierarchical constraint implies that the configuration with all the spins in the down (-1) state is completely absorbent, i. e., the system does not evolve in time from that configuration.

We are interested in the time evolution of the average spin $`\sigma _i^{(n)}`$, which is given in Glauber dynamics by

$$\frac{d}{dt}\sigma _i^{(n)}=2\sigma _i^{(n)}W_i^{(n)}(𝝈),$$
(9)
and substitution of Eqs. (7) and (8) into this expression yields

$$\frac{d}{dt}\sigma _i^{(n)}=\left[\alpha _{k_i}^{(n1)}(𝝈)\alpha _{k_i+1}^{(n1)}(𝝈)\mathrm{}\alpha _{k_i+\mu _{n1}1}^{(n1)}(𝝈)\right]^{1/\mu _{n1}}\left(\sigma _i^{(n)}\sigma _e\right)\underset{j=k_i}{\overset{k_i+\mu _{n1}1}{}}\delta _{\sigma _j^{(n1)},+1},$$
(10)
where we have taken into account that $`(\sigma _i^{(n)})^2=+1`$ for all $`i`$, $`n`$. This equation is rather involved, since it couples the evolution of $`\sigma _i^{(n)}`$ to moments of the probability distribution containing an increasing number of spins of all the levels $`n^{}`$ such that $`0n^{}<n`$. The levels $`0n^{}n2`$ enter into the equation through the rates $`\alpha _j^{(n1)}(𝝈)`$. Then, we introduce at this stage a sort of “mean field” approximation for the transition rates, which is closely related, in spirit, to the seminal work of Palmer et al. . Within this coarsening dynamics, it is reasonable to expect that spins in level $`n`$ evolve over a time scale quite larger than that characteristic of level $`n1`$. This means that spins in level $`n1`$ change many times their state before a transition in level $`n`$ takes place. Thus, we replace the product of the Kronecker deltas in Eq. (8) by its average value, i. e., by the probability $`P^{(n1)}(\sigma _{k_i}^{(n1)}=+1,\mathrm{},\sigma _{k_i+\mu _{n1}1}^{(n1)}=+1)`$ that all the $`\mu _{n1}`$ spins of the given cluster be in the up state,

$`\alpha _i^{(n)}(𝝈)`$ $`=`$ $`\left(\alpha _{k_i}^{(n1)}(𝝈)\alpha _{k_i+1}^{(n1)}(𝝈)\mathrm{}\alpha _{k_i+\mu _{n1}1}^{(n1)}(𝝈)\right)^{1/\mu _{n1}}`$ (12)
$`\times P^{(n1)}(\sigma _{k_i}^{(n1)}=+1,\mathrm{},\sigma _{k_i+\mu _{n1}1}^{(n1)}=+1).`$
Moreover, we will restrict ourselves to situations where there is spatial homogeneity within each of the levels, so that the dependence on the specific site considered in a given level can be dropped, obtaining

$$\alpha ^{(n)}(𝝈)=\alpha ^{(n1)}(𝝈)P^{(n1)}(\sigma _1^{(n1)}=+1,\mathrm{},\sigma _{\mu _{n1}}^{(n1)}=+1).$$
(13)
Iteration of the above relation gives

$$\alpha ^{(n)}(𝝈)=\alpha ^{(0)}\underset{j=0}{\overset{n1}{}}P^{(j)}(\sigma _1^{(j)}=+1,\mathrm{},\sigma _{\mu _j}^{(j)}=+1),$$
(14)
with $`\alpha ^{(0)}`$ being a constant, that characterizes the relaxation rate of the spins belonging to level $`n=0`$, whose dynamics is not constrained. As $`\alpha ^{(0)}`$ determines the basic time scale, which is arbitrary, we will take $`\alpha ^{(0)}=1`$ in the following. In the mean field approximation just introduced, the time evolution of the average value of the spin, Eq. (10), takes the form

$$\frac{d}{dt}\sigma ^{(n)}=\underset{j=0}{\overset{n1}{}}P(\sigma _1^{(j)}=+1,\mathrm{},\sigma _{\mu _j}^{(j)}=+1)\left(\sigma ^{(n)}\sigma _e\right).$$
(15)
Now, a new approximation will be made. The probability $`P^{(j)}`$ in Eq. (15) will be substituted by its equilibrium value. This would be exact in linear response around equilibrium, but it will be taken here as an approximation leading to the basic equation of our hierarchically constrained model, namely

$$\frac{d}{dt}\sigma ^{(n)}=\underset{j=0}{\overset{n1}{}}p_e^{\mu _j}\left(\sigma ^{(n)}\sigma _e\right),$$
(16)
where $`p_e`$ is the equilibrium probability of any spin being in the up state, given by Eq. (5). Again, on physical grounds, this is a sensible approximation due to the separation of the characteristic time scales of the different levels. Because of the hierarchically constrained dynamics, spins in level $`n1`$ reach equilibrium over a time scale in which spins in level $`n`$ have not begun to evolve. However, as a consequence of the last approximation, the configuration with all the spins in the down state is no longer absorbent.

From Eq. (16), an equivalent equation can be written for the evolution of the probability $`p^{(n)}`$ of the up state in level $`n`$, defined in Eq. (3), i. e.,

$$\frac{d}{dt}p^{(n)}=\alpha _n\left(p^{(n)}p_e\right),$$
(17)
being $`\alpha _n`$ the characteristic relaxation rate of level $`n`$,

$$\alpha _n=p_e^{g_n},g_n=\underset{j=0}{\overset{n1}{}}\mu _j.$$
(18)
Equation (17) implies that, due to the hierarchical constraints, the spin relaxation slows down with increasing level $`n`$, since $`\alpha _n`$ is a decreasing function of $`n`$, because $`p_e<1`$. Eq. (17) is the main result in this section. In the following, we will explore its implications, considering first processes at constant temperature in the next section, and cooling and heating processes in the remainder of the paper.

## III Relaxation at constant temperature

For the case of constant temperature $`T`$, and therefore constant $`\alpha _n`$, Eq. (17) is easily solved,

$$p^{(n)}(t)p_e=\left[p^{(n)}(0)p_e\right]e^{\alpha _nt}.$$
(19)
Then, each spin relaxes exponentially to equilibrium with the rate characteristic of its level.

For the homogenous situations within each level we are considering, the dimensionless mean energy per spin  defined in Eq. (2) simplifies to

$$\epsilon (t)=\underset{n=0}{\overset{n_{\mathrm{max}}}{}}w_np^{(n)}(t),$$
(20)
where $`w_n=N_n/N`$ is the fraction of spins in level $`n`$, verifying $`_{n=0}^{n_{\mathrm{max}}}w_n=1`$.

Use of Eq. (19) into Eq. (20) yields

$$\epsilon (t)=\epsilon _e(T)+\underset{n=0}{\overset{n_{\mathrm{max}}}{}}w_n\left[p^{(n)}(0)p_e\right]e^{\alpha _nt},$$
(21)
for the relaxation of the energy at constant temperature. In order to proceed, we will consider the simple case in which the initial probability distribution $`p^{(n)}(0)`$ does not depend on the index level $`n`$. This will be the situation, for instance, when the initial state corresponds to equilibrium at a different temperature $`T+\mathrm{\Delta }T`$. Thus, the relaxation function of the physical property described by the hamiltonian of the system is given by

$$\varphi (t)\frac{\epsilon (t)\epsilon _e}{\epsilon (0)\epsilon _e}=\underset{n=0}{\overset{n_{\mathrm{max}}}{}}w_ne^{\alpha _nt}.$$
(22)
This equation is a generalization of the result derived by Palmer et al. in their pioneering work in hierarchically constrained dynamics , which corresponds to the choice $`p_e=1/2`$. This is formally equivalent to the particularization of Eq. (22) for $`T\mathrm{}`$. If other positive values of the temperature are considered, the effect is an increase of the relaxation times

$$\tau _n=\alpha _n^1,$$
(23)
since $`p_e`$ is a decreasing function of the temperature. A main advantage of the formulation of the hierarchical models as presented here, aside from its larger generality, is that it allows the analysis of processes in which the temperature of a thermal system (or the vibration intensity in a granular system) changes in time. This kind of processes will be addressed in the next section. A mean relaxation time $`\tau `$ can be defined as

$$\tau =_0^{\mathrm{}}𝑑t\varphi (t)=\underset{n=0}{\overset{n_{\mathrm{max}}}{}}w_n\tau _n,$$
(24)
providing a quantitative measure of the time it takes the system to relax to equilibrium at temperature $`T`$.

Let us consider that the fraction of spins in level $`n`$, $`w_n`$, and the number of “facilitating” spins in level $`n`$, $`\mu _n`$, depend very smoothly on $`n`$, i. e., they can be expressed as functions of the form

$$w_n=w(n\eta )\mu _n=\mu (n\eta ),$$
(25)
where $`\eta 1`$. These seem to be sensible conditions when modelling a real system, in which the introduction of the levels and the pseudospins is associated to some coarse-grained description. By defining

$$x_n=n\eta ,$$
(26)
that is a continuous variable in the limit $`\eta 0`$, the sums over $`n`$ can be replaced by integrals. The relaxation rate of level $`n`$, given by Eq. (18), becomes a function of the continuous variable $`x`$,

$$\alpha (x)=p_e^{g(x)},$$
(27)
with

$$g(x)=\frac{_{k=0}^{n1}\mu _k\eta }{_{k=0}^{n_{\mathrm{max}}}w_k\eta }=\frac{_0^x𝑑x^{}\mu (x^{})}{_0^{x_{\mathrm{max}}}𝑑x^{}w(x^{})},$$
(28)
where $`x_{\mathrm{max}}=n_{\mathrm{max}}\eta `$, and the normalization of the weights $`w_n`$ has been used. Therefore, the relaxation function $`\varphi (t)`$, given by Eq. (22), becomes

$$\varphi (t)=\frac{_{n=0}^{n_{\mathrm{max}}}\eta w_ne^{\alpha _nt}}{_{n=0}^{n_{\mathrm{max}}}\eta w_n}=\frac{_0^{x_{\mathrm{max}}}𝑑xw(x)e^{\alpha (x)t}}{_0^{x_{\mathrm{max}}}𝑑xw(x)},$$
(29)
in the continuous limit.

In general, Eq. (29) is mathematically rather involved, since it depends both on the functions $`\mu (x)`$ and $`w(x)`$. The simplest possibility appears to be $`\mu (x)`$ proportional to $`w(x)`$, i. e., $`\mu _n`$ proportional to $`w_n`$ or, equivalently, to $`N_n`$. In other words, the simplest kind of hierarchically constrained models shows up when the number of “facilitating” spins at a level is an extensive function of the number of spins at the same level . This condition is expressed as

$$\mu (x)=aw(x),$$
(30)
with $`a`$ being a constant, independent of $`x`$. In this case, it is useful to define the new variable

$$u=\frac{_0^x𝑑x^{}w(x^{})}{_0^{x_{\mathrm{max}}}𝑑x^{}w(x^{})},$$
(31)
measuring the fraction of spins belonging to levels up to $`n=x/\eta `$. In terms of $`u`$, the relaxation rates of Eq. (27) are given by

$$\alpha (u)=p_e^{g(u)},g(u)=au,$$
(32)
the relaxation function is expressed as

$$\varphi (t)=_0^1𝑑ue^{\alpha (u)t},$$
(33)
and the mean relaxation time reads

$$\tau =_0^{\mathrm{}}𝑑t\varphi (t)=_0^1𝑑u\tau (u)=\frac{p_e^a1}{a|\mathrm{ln}p_e|},$$
(34)
with $`\tau (u)=\alpha ^1(u)`$. It is interesting to consider situations for which $`p_e^a1`$, so that the minimum relaxation rate $`\alpha (1)`$ is much smaller than the maximum one $`\alpha (0)=1`$ in our dimensionless time scale. In this case, the relaxation function $`\varphi (t)`$ is linear in $`\mathrm{ln}t`$ over an intermediate time window, $`1tp_e^a`$, namely

$$\varphi (t)1\frac{1}{a|\mathrm{ln}p_e|}\left(\gamma +\mathrm{ln}t\right),$$
(35)
where $`\gamma `$ stands for Euler’s constant, $`\gamma 0.577`$. This kind of linear logarithmic behavior is characteristic of a great variety of complex systems, including spin-glasses , granular materials , structural glasses , and protein models . In the present context, the condition $`p_e^a1`$ corresponds to a “low” temperature limit, in which the mean relaxation time, given by Eq. (34), is very large, i. e., the relaxation of the system becomes very slow.

It is worth noting that, in terms of the $`u`$ variable, the continuous limit is formally obtained by changing the functions of the index level $`n`$ by the corresponding functions of the variable $`u`$, and by making the replacement

$$\underset{n=0}{\overset{n_{\mathrm{max}}}{}}w_n_0^1𝑑u.$$
(36)
It follows that, in the continuous limit, the dynamical behavior of the system does not depend explicitly on the level populations $`N_n`$, but only on the relaxation rates $`\alpha `$ expressed as functions of $`u`$, as given by Eq. (32).

## IV Time-dependent temperature processes

In this section processes in which the temperature changes in time will be studied. The evolution equation (17) is now

$$\frac{d}{dt}p^{(n)}(t)=\alpha _n(T)\left[p^{(n)}(t)p_e(T)\right],$$
(37)
where $`T=T(t)`$. The general solution of Eq. (37) is

$$p^{(n)}(t)=\left[p^{(n)}(t_0)p_e(T_0)\right]\chi _n(t,t_0)+p_e(T)_{t_0}^t𝑑t^{}\frac{dp_e(T^{})}{dT^{}}\frac{dT^{}}{dt^{}}\chi _n(t,t^{}).$$
(38)
Here $`T_0=T(t_0)`$ is the initial value of the temperature, $`T=T(t)`$, $`T^{}=T(t^{})`$, $`T^{\prime \prime }=T(t^{\prime \prime })`$, and we have introduced the function

$$\chi _n(t_1,t_2)=e^{_{t_2}^{t_1}𝑑t\alpha _n(t)}.$$
(39)
The above equation is valid for any law of variation for the temperature. Taking into account Eq. (20), the average energy per spin $`\epsilon `$, is given by

$`\epsilon (t)`$ $`=`$ $`{\displaystyle \underset{n=0}{\overset{n_{\mathrm{max}}}{}}}w_n\left[p^{(n)}(t_0)p_e(T_0)\right]\chi _n(t,t_0)`$ (41)
$`+\epsilon _e(T){\displaystyle _{t_0}^t}𝑑t^{}{\displaystyle \frac{d\epsilon _e(T^{})}{dT^{}}}{\displaystyle \frac{dT^{}}{dt^{}}}{\displaystyle \underset{n=0}{\overset{n_{\mathrm{max}}}{}}}w_n\chi _n(t,t^{}),`$
where $`\epsilon _e`$ is the average equilibrium energy, defined in Eq. (6).

Let us assume that the system is initially at equilibrium with $`T=T_0`$. Then the first term on the rhs of Eq. (41) vanishes, and

$$\epsilon (t)=\epsilon _e(T)_{t_0}^t𝑑t^{}\frac{d\epsilon _e(T^{})}{dT^{}}\frac{dT^{}}{dt^{}}M(t,t^{}),$$
(42)
where

$$M(t,t^{})=\underset{n=0}{\overset{n_{\mathrm{max}}}{}}w_n\chi _n(t,t^{}).$$
(43)
is a memory function. In the case of constant temperature, $`M(t,t^{})`$ is equal to the relaxation function $`\varphi (tt^{})`$, as seen by comparing Eq. (43) with Eq. (22). The structure of Eq. (42) is the same as that of Narayanaswami’s phenomenological theory of glasses . A similar result was obtained some years ago for the one dimensional Ising model with Glauber dynamics .

### A High temperature limit: Hilbert’s method

We are going to look for a solution of Eq. (37) by means of Hilbert’s method. A special solution

$$p_H^{(n)}(t)=\underset{k=0}{\overset{\mathrm{}}{}}p_H^{(n),k}(t)$$
(44)
is constructed in an iterative way as follows. We take

$$p_H^{(n),0}(t)=p_e(T),,$$
(45)
while for $`k1`$

$$p_H^{(n),k}(t)=\alpha _n^1(T)\frac{dp^{(n),k1}(t)}{dt}.$$
(46)
Equation (45) shows that Hilbert’s expansion agrees with the equilibrium distribution to the lowest order. Besides, for $`k=1`$ we get from Eq. (46)

$$p_H^{(n),1}(t)=\tau _n(T)\frac{dp_e(T)}{dT}\frac{dT}{dt}.$$
(47)
This equation indicates the main limitation of Hilbert’s method. Due to the divergence of the relaxation times $`\tau _n=\alpha _n^1`$ in the low temperature limit, see Eq. (18), also $`p^{(n),1}`$ diverges in that limit. As a consequence, Hilbert’s solution is only accurate in the high temperature regime, in which an expansion around equilibrium provides a good approximation. Restricting ourselves to high temperatures, we approximate

$$p_H^{(n)}(T)p_e(T)\tau _n(T)\frac{dp_e(T)}{dT}\frac{dT}{dt},$$
(48)
and, from Eq. (20),

$$\epsilon _H(T)\epsilon _e(T)\frac{d\epsilon _e(T)}{dT}\frac{dT}{dt}\underset{n=0}{\overset{n_{\mathrm{max}}}{}}w_n\tau _n(T).$$
(49)
Taking into account the definition of the average relaxation time $`\tau `$, Eq. (24), the above expression is seen to be equivalent to

$$\epsilon _H(T)\epsilon _e(T)\frac{d\epsilon _e(T)}{dT}\frac{dT}{dt}\tau (T),$$
(50)
which agrees with the high temperature behavior of Eq. (42). In Eqs. (48) and (50), $`p_H^{(n)}`$ and $`\epsilon _H`$ depend on time only through the temperature $`T(t)`$. Thus, Hilbert’s method provides a “normal” solution, in the sense often used in kinetic theory. The validity of a expression identical to Eq. (50), also in the high temperature limit, has been established for a quite general class of systems whose dynamics is described by a master equation . Although we have made here several drastic approximations in order to get a closed equation for the average spin, the high temperature limit of the solution, given by Hilbert’s method, remains formally the same as that of the exact solution of the original model. Certainly, this is a good property of those approximations.

What is the physical meaning of the failure of Hilbert’s expansion for low temperatures? Due to the divergence of the characteristic relaxation times, the system does not have enough time to relax to the equilibrium curve at very low temperatures, and it gets “frozen” in a far from equilibrium state. Since Hilbert’s method is an expansion around equilibrium, it fails in the low temperature region. In fact, the rhs of Eq. (48) becomes negative for low enough temperatures. On the other hand, Hilbert’s expansion is useful to estimate the values of the physical properties in the “frozen” state . Also, Hilbert’s method provides a qualitative understanding of the hysteresis effects appearing in thermal cycles (cooling and reheating). In cooling processes ($`dT/dt<0`$), it is $`\epsilon _H\epsilon _e`$, while in heating processes ($`dT/dt>0`$) it is $`\epsilon _H\epsilon _e`$. Then, $`\epsilon _H`$ lies to opposite sides of the equilibrium curve for cooling and heating processes, and hysteresis effects show up in thermal cycling experiments, as it will be discussed in more detail in Sec. VI.

## V Cooling processes

Next, we are going to study the continuous cooling of the system down to very low temperatures. The origin of time is taken at the beginning of the cooling process. The initial condition will be the equilibrium configuration at a “high” temperature $`T_0`$, i. e.,

$$p^{(n)}(0)=p_e(T_0).$$
(51)
Then, the first correction in Hilbert’s expansion, $`p^{(n),1}(t)`$, is very small as compared with $`p_e(T_0)`$ for $`TT_0`$. Particularization of Eq. (38) for the above initial condition gives

$$p^{(n)}(t)=p_e(T)_0^t𝑑t^{}\frac{dp_e(T^{})}{dT^{}}\frac{dT^{}}{dt^{}}\chi _n(t,t^{}).$$
(52)
Since $`p_e(T)`$ is an increasing function of the temperature, for continuous cooling processes it is

$$p^{(n)}(t)p_e(T)\text{for all }t\text{ and }n\text{.}$$
(53)
The possible deviations from the equilibrium distribution always lead to an increase of the probability of the spin being in the excited state. Moreover, Eq. (52) directly implies that

$$\epsilon (t)=\epsilon _e(T)_0^t𝑑t^{}\frac{d\epsilon _e(T^{})}{dT^{}}\frac{dT^{}}{dt^{}}\underset{n=0}{\overset{n_{\mathrm{max}}}{}}w_n\chi _n(t,t^{})\epsilon _e(T),$$
(54)
where we have used Eq. (20). This inequality has been experimentally observed in glass-forming liquids . Since the reported experiments were made at constant pressure, the quantity $`\epsilon `$ considered here must be interpreted as the enthalpy in that context.

In order to proceed further in our analysis, the continuous limit introduced in the study of the relaxation at constant temperature in Sec. III will be considered. As already mentioned, this continuous limit is expected to be closer to the description of real systems than the discrete level picture. Besides, for the sake of concreteness, we will restrict ourselves to those models verifying Eq. (30). The index level $`n`$ is substituted by the continuous variable $`u`$, defined in Eq. (31), representing the fraction of the total number of spins up to level $`n`$. With an obvious change of notation, Eq. (52) becomes

$$p(t;u)=p_e(T)_0^t𝑑t^{}\frac{dp_e(T^{})}{dT^{}}\frac{dT^{}}{dt^{}}\chi (t,t^{};u),$$
(55)
where

$$\chi (t,t^{};u)=e^{_t^{}^t𝑑t^{\prime \prime }\alpha (T^{\prime \prime };u)},$$
(56)
with $`\alpha (T;u)`$ given by Eq. (32), i. e.,

$$\alpha (T;u)=p_e(T)^{au}.$$
(57)
Also, using Eq. (36), it is found

$$\epsilon (t)=\epsilon _e(T)_0^t𝑑t^{}\frac{d\epsilon _e(T^{})}{dT^{}}\frac{dT^{}}{dt^{}}_0^1𝑑u\chi (t,t^{};u).$$
(58)
The time evolution of the probability $`p(t;u)`$ depends on the explicit form of the cooling law. In experiments, linear cooling is usually employed,

$$\frac{dT}{dt}=r_c,$$
(59)
where $`r_c>0`$ is the cooling rate determining the time scale, $`r_c^1`$, over which the temperature changes. Linear cooling implies that $`p_e(T)`$ depends on time in a rather involved way. From Eqs. (5) and (59) one gets

$$\frac{dp_e}{dt}=\frac{1}{2}r_cp_e(1p_e)\left(\mathrm{ln}\frac{p_e}{1p_e}\right)^2.$$
(60)
We are interested in cooling processes for which the temperature changes slowly in time, $`r_c1`$, so that the system departs from the equilibrium curve for very low temperatures, where $`p_e1`$. Then, a law equivalent to linear cooling, aside from logarithmic corrections, is

$$\frac{dp_e}{dt}=r_cp_e,$$
(61)
having the advantage that analytical calculations are much more simple with Eq. (62) . In the following, we will use the notation

$$p_e(T)=p_e,p_e(T^{})=p_e^{},p_e(T_0)=p_{e0}.$$
(62)
The time integrals in Eq. (55) can be transformed into integrals over $`p_e`$ by means of the cooling law (61), with the result

$$p(t;u)=p_e+_{p_e}^{p_{e0}}𝑑p_e^{}\mathrm{exp}\left(\frac{p_e^{au}p_e^{au}}{r_cau}\right),$$
(63)
The second term in this expression is dominant in the low temperature region, where $`p_e`$ is very small. Therefore, it follows that spins in any level $`u`$ fall in a nonequilibrium state for low enough temperatures. The details of this glass-like transition will be analyzed below, in a separate subsection.

One of the main quantities characterizing a given cooling process is the residual value $`f_{\text{res}}`$ of a relevant property $`f`$. The residual value measures the excess with respect to the equilibrium curve, extrapolated to very low temperatures,

$$f_{\text{res}}=\underset{T0}{lim}(ff_e).$$
(64)
In particular, for the probability $`p(t;u)`$ we have from Eq. (63)

$$p_{\text{res}}(u)=_0^{p_{e0}}𝑑p_e^{}\mathrm{exp}\left(\frac{p_e^{au}}{r_cau}\right),$$
(65)
which is easily transformed into

$$p_{\text{res}}(u)=\frac{1}{au}\left(r_cau\right)^{\frac{1}{au}}_0^{x_0}𝑑xx^{\frac{1}{au}1}e^x,$$
(66)
where $`x_0=p_{e0}^{au}/(r_cau)`$. The slow cooling limit is defined by the residual properties being independent of the initial conditions, and determined univocally by the cooling rate . In our case, slow cooling means that the upper integration limit in Eq. (66) can be substituted by infinity for all $`u`$, i. e.,

$$\frac{p_{e0}^{au}}{r_cau}1.$$
(67)
As the $`u`$ variable varies in the interval $`0u1`$, the slow cooling condition is

$$r_cap_{e0}^a<1.$$
(68)
Then, with an exponentially small error,

$$p_{\text{res}}(u)\frac{1}{au}\left(r_cau\right)^{\frac{1}{au}}\mathrm{\Gamma }\left(\frac{1}{au}\right)=\left(r_cau\right)^{\frac{1}{au}}\mathrm{\Gamma }\left(1+\frac{1}{au}\right).$$
(69)
This expression gives the probability that the spins in level $`u`$ be in the up state at very low temperatures. In Fig. 2, the residual probability $`p_{\text{res}}(u=1)`$ is plotted as a function of the cooling rate $`r_c`$, for $`a=1`$. The asymptotic result, given by Eq. (69), is compared with the numerical integration of Eq. (66) with $`p_{e0}=1/2`$, i. e., the system is taken initially at infinite temperature. The agreement is quite good up to $`r_c0.1`$, which is not very small.

The evolution of the average energy, for the cooling law (61), is given by

$$\epsilon (t)\epsilon _e(T)=_0^1𝑑u\left[p(t;u)p_e\right]=_0^1𝑑u_{p_e}^{p_{e0}}𝑑p_e^{}\mathrm{exp}\left(\frac{p_e^{au}p_e^{au}}{r_cau}\right).$$
(70)
The residual energy can be easily computed by particularizing this expression for $`T0`$,

$$\epsilon _{\text{res}}=_0^1𝑑up_{\text{res}}(u).$$
(71)
The integrand $`p_{\text{res}}(u)`$, given by Eq. (69), vanishes exponentially in the limit $`u0`$. A standard Laplace analysis can be made, with the result

$$\epsilon _{\text{res}}\mathrm{\Gamma }\left(1+\frac{1}{a}\right)\frac{\left(r_ca\right)^{\frac{1}{a}}}{\left|\mathrm{ln}\left(r_ca\right)^{\frac{1}{a}}\right|}.$$
(72)
The leading behavior is potential with $`r_c`$, since

$$\mathrm{ln}\epsilon _{\text{res}}\frac{1}{a}\mathrm{ln}\left(r_ca\right),$$
(73)
which comes from the upper limit of integration, $`u=1`$, corresponding to the largest relaxation time. The integral over the whole distribution function $`p_{\text{res}}(u)`$ gives a logarithmic correction $`\left|\mathrm{ln}(r_ca)\right|^1`$, which makes the residual energy smaller than the dominant term $`(r_ca)^{1/a}`$. This is due to the increasing behavior of $`p_{\text{res}}(u)`$ with $`u`$, $`p_{\text{res}}(u)p_{\text{res}}(u=1)`$. In Fig. 3 the residual value of the energy is plotted. The asymptotic expression, Eq. (72), is compared with the numerical results from Eqs. (71) and (66). A good agreement is found up to $`r_c=0.01`$. It is worth noting that, for the cooling law considered, the logarithmic correction to the potential in $`r_c`$ behavior is not present for other simple models of structural glasses previously studied .

### A Demarcation mode, fictive temperature and glass transition

The existence of non-vanishing residual properties is an indication of the departure of the system from equilibrium at low temperatures. Due to the divergence of the characteristic relaxation times $`\tau _n`$, for low enough temperatures the system does not have enough time to relax towards equilibrium, and a kinetic phenomenon resembling the laboratory glass transition  shows up. Next, we will try to understand the physical origin of this kinetic transition.

Let us consider again the time evolution of the probability distribution $`p(t;u)`$ in a cooling process, as given by Eqs. (55) and (56). The integral in $`\chi (t,t^{};u)`$,

$$I(t,t^{};u)=_t^{}^t𝑑t^{\prime \prime }\alpha (T^{\prime \prime };u)$$
(74)
is a measure of the average number of transitions occurring in level $`u`$ in the time interval between $`t^{}`$ and $`t`$. Consequently, a mathematical definition for the limit of slow cooling is that the condition

$$I(t^{},0;u)>>1$$
(75)
holds for all $`u`$, where $`t^{}`$ is the time for which the temperature vanishes if extrapolated accordingly to the prescribed cooling law, i. e., $`T(t^{})=0`$. Eq. (75) guarantees that the system experiments a large number of transitions before getting eventually frozen, so that it has enough time to forget the details of the initial condition. For the cooling law defined in Eq. (61) it is easily verified that Eq. (75) is equivalent to Eq. (67).

If we are dealing with a slow cooling process, Eq. (75) implies that there is a time window over which

$$I(t,0;u)>>1.$$
(76)
This is the time regime we are interested in. Let us analyze the behavior of $`\chi (t,t^{};u)`$ as a function of $`t^{}`$, $`t_0t^{}t`$, for a given time $`t`$ such that Eq. (76) holds. The function $`I(t,t^{};u)`$ changes from a very large value to zero when $`t^{}`$ goes from $`0`$ to $`t`$. Consequently, $`\chi (t,t^{};u)`$ increases from practically zero to unity when $`t^{}`$ moves in the above time interval. Let us define a time $`t_f(t;u)`$, prior to $`t`$, by

$$I(t,t_f;u)=1,$$
(77)
so that the average number of transitions taking place in level $`u`$ in the time interval between $`t_f`$ and $`t`$ equals unity. Then, $`\chi (t,t_f;u)=e^1`$, and the function $`\chi (t,t^{};u)`$ changes from zero to unity in a certain time interval around $`t_f`$. In order to proceed, we will assume that this change takes place in the vicinity of $`t_f`$ very rapidly, as compared with the variation of the rest of the integrand of Eq. (55). Decomposing Eq. (55) in the form

$$p(t;u)=p_e(T)_0^{t_f}𝑑t^{}\frac{dp_e(T^{})}{dT^{}}\frac{dT^{}}{dt^{}}\chi (t,t^{};u)_{t_f}^t𝑑t^{}\frac{dp_e(T^{})}{dT^{}}\frac{dT^{}}{dt^{}}\chi (t,t^{};u),$$
(78)
the first integral is subdominant with respect to the second one and, moreover,

$$_{t_f}^t𝑑t^{}\frac{dp_e(T^{})}{dT^{}}\frac{dT^{}}{dt^{}}\chi (t,t^{};u)_{t_f}^t𝑑t^{}\frac{dp_e(T^{})}{dT^{}}\frac{dT^{}}{dt^{}}=p_e(T)p_e[T_f(t;u)],$$
(79)
where $`T_f(t;u)`$ is the temperature of the system at time $`t_f(t;u)`$, i. e.,

$$T_f(t;u)=T[t_f(t;u)].$$
(80)
Note that $`T_f(t;u)>T`$, since the time instant $`t_f(t;u)<t`$. Substitution of Eq. (79) into Eq. (78), and use of the above approximations yields

$$p(t;u)=p_e[T_f(t;u)].$$
(81)
The arguments leading from Eq. (55) to Eq. (81) are formally equivalent to assume that

$$\chi (t,t^{};u)=\mathrm{\Theta }[t^{}t_f(t;u)],$$
(82)
where $`\mathrm{\Theta }(x)`$ is Heaviside’s step function. The result in Eq. (81) has a neat physical interpretation: the probability distribution of spins in level $`u`$ at time $`t`$ is given by the equilibrium distribution corresponding to the temperature of the system at a prior time $`t_f(t;u)`$, defined in Eq. (77). The temperature $`T_f(t;u)`$ is, therefore, the “fictive” temperature of level $`u`$ for the cooling program under consideration, in the sense used in the phenomenological theories of glass-forming liquids . In the present context, each level has its own fictive temperature, so that there is not a unique equilibrium distribution describing the whole state of the system at a given temperature. As a consequence, the behavior of different macroscopic properties of the system may be quite different, depending on the modes which dominate for the evaluation of each property.

The definition of the fictive temperature, Eq. (80), can be written as

$$I(t^{},t_f;u)=I(t^{},t;u)+I(t,t_f;u)=I(t^{},t;u)+1.$$
(83)
In the “high temperature” regime, where the number of transitions between temperature $`T`$ and $`T=0`$ in level $`u`$ is very large, we can neglect the unity on the rhs of Eq. (83) and we get $`T_f(t;u)T`$; the modes in level $`u`$ remain in equilibrium. On the other hand, in the “low temperature” limit, $`I(t^{},t;u)`$ becomes very small and $`T_f(t;u)`$ tends to a limiting value,

$$T_f^{}(u)=\underset{tt^{}}{lim}T_f(t;u)=T[t_f^{}(u)],$$
(84)
where $`t_f^{}(u)`$ is the time for which

$$I(t^{},t_f^{};u)=_{t_f^{}}^t^{}𝑑t\alpha (T;u)=1.$$
(85)
The low temperature region for level $`u`$ corresponds to temperatures such that the average number of transitions remaining to spins in level $`u`$ before formally reaching $`T=0`$ is smaller than unity. Since the rate $`\alpha (T;u)`$ decreases as a function of $`u`$, see Eq. (32), $`T_f^{}(u)`$ is an increasing function of $`u`$, i. e., the fictive temperature $`T_f^{}(u)`$ is larger the slower the modes. Then, the following picture of the kinetic glass-like transition appears. The relaxation modes associated to level $`u`$ are in equilibrium for $`TT_f^{}(u)`$, getting frozen in their equilibrium distribution corresponding to the fictive temperature $`T_f^{}(u)`$ for $`TT_f^{}(u)`$. A similar idea was first used in the context of the structural glass transition by Dyre , who characterized the transition by means of the concept of demarcation mode. For a given value of the temperature $`T`$, the modes in levels with $`T_f^{}(u)>T`$ will be frozen, while those with $`T_f^{}(u)<T`$ will be still able to relax. The demarcation mode at temperature $`T`$, $`u_d(T)`$, is defined by $`T_f^{}(u_d)=T`$. The modes with $`u>u_d(T)`$ are frozen at this temperature, and modes with $`u<u_d(T)`$ can still evolve. The glass transition can be understood as the movement of the demarcation mode from $`u=1`$, corresponding to the largest relaxation time, to $`u=0`$, which is the fastest mode. The temperature $`T_g`$ at which the glass transition begins for a given cooling program is that for which the slowest modes freeze, i. e., $`T_g=T_f^{}(u=1)`$.

The physical image developed in the above paragraph allows us to estimate the residual values in quite a simple way. For a given cooling law, first we calculate the limit value of the fictive temperature $`T_f^{}(u)`$. Then, the residual probability distribution is given by,

$$p_{\text{res}}(u)p_e[T_f^{}(u)],$$
(86)
and, consequently, the residual value of the energy is

$$\epsilon _{\text{res}}=_0^1𝑑up_{\text{res}}(u)_0^1𝑑up_e[T_f^{}(u)].$$
(87)
As a test, we have considered the cooling process of the hierarchical model studied in this paper with the law given by Eq. (61). In the slow cooling limit, it is easy to show that

$$T_f^{}(u)\frac{2au}{|\mathrm{ln}(r_ca)|},$$
(88)
i. e., the glass-like transition begins at

$$T_g\frac{2a}{|\mathrm{ln}(r_ca)|},$$
(89)
for which the slowest modes become frozen. In Fig. 4 we plot the evolution of the probability, as given by Eq. (63), of the slowest level $`u=1`$, for $`a=1`$ and $`r_c=10^4`$. For these parameters, it is $`T_g=T_f^{}(u=1)0.217`$, which is seen to be a good value for the fictive temperature. It is worth noting that the simple arguments used in this section lead to an estimate of the residual population

$$p_{\text{res}}(u)p_e\left[T_f^{}(u)\right]=\left(r_cau\right)^{\frac{1}{au}},$$
(90)
that gives the same dependence of $`p_{\text{res}}`$ on $`r_c`$ as Eq. (69), except for a factor of the order of unity.

In Fig. 5 the dimensionless average energy $`\epsilon `$, Eq. (70), is shown, for the same parameters as in Fig. 4. The predicted glass transition temperature $`T_g0.217`$ is also indicated. It provides a good estimate for the beginning of the departure from the equilibrium curve. The residual value of the energy calculated from an asymptotic analysis of Eqs. (87) and (90) reads

$$\epsilon _{\text{res}}=_0^1𝑑up_{\text{res}}(u)\frac{\left(r_ca\right)^{\frac{1}{a}}}{\left|\mathrm{ln}\left(r_ca\right)^{\frac{1}{a}}\right|}.$$
(91)
Comparison of Eqs. (91) and Eq. (72) also indicates a good agreement, the difference being again a factor of the order of unity. We can introduce a global fictive temperature $`\overline{T}_f^{}`$ for the energy, as the temperature for which the equilibrium energy is the same as the energy of the system in the cooling process, when extrapolated to $`T=0`$, i. e.,

$$\epsilon _{\text{res}}=\epsilon _e(\overline{T}_f^{}),$$
(92)
This global fictive temperature $`\overline{T}_f^{}`$ is related to the fictive temperatures of the levels $`T_f^{}(u)`$ by

$$\epsilon _e(\overline{T}_f^{})=_0^1𝑑up_e[T_f^{}(u)].$$
(93)
The global fictive temperature is then a kind of average of the fictive temperatures $`T_f^{}(u)`$, but the probability distribution of the frozen state is not the equilibrium distribution at the temperature $`\overline{T}_f^{}`$. Each level freezes at its corresponding fictive temperature $`T_f^{}(u)`$, and the probability of the frozen state $`p_{\text{res}}(u)`$ is approximately given by Eq. (90). For the values of the parameters considered in Fig. 5, Eq. (91) leads to $`\overline{T}_f^{}0.175`$, which is also a good approximation for the value of the temperature at which Eq. (92) holds. Finally, let us stress that similar results for the fictive temperatures can be derived by means of the qualitative reasonings developed in Ref. , and based on Hilbert’s expansion.

## VI Heating processes

This section will be devoted to the study of heating processes after the system has been previously cooled down to very low temperatures. Then, the system is not initially at equilibrium, but the initial probability distribution $`p(t_0;u)`$ corresponds to the final state reached in the cooling process. It is convenient to introduce a time $`t_h<t_0`$ such that $`T(t_h)=0`$ if the heating law is extrapolated to times shorter than $`t_0`$. From Eq. (38) it follows that

$`p(t;u)`$ $`=`$ $`\left[p(t_0;u)p_e(T_0)\right]\chi (t,t_0;u)+p_e(T)`$ (95)
$`{\displaystyle _{t_0}^t}𝑑t^{}{\displaystyle \frac{dp_e(T^{})}{dT^{}}}{\displaystyle \frac{dT^{}}{dt^{}}}\chi (t,t^{};u),`$
where $`\chi (t_1,t_2;u)`$ was defined in Eq. (56), and $`T_0=T(t_0)`$ is the initial temperature of the heating process. The first term on the rhs of Eq. (95) represents the decay of the initial nonequilibrium condition. For $`T(t)`$ in the high temperature region it is $`\chi (t,t_0;u)1`$, and this contribution can be neglected. This implies that $`p(t;u)`$ reaches a behavior which is independent of the initial condition or, equivalently, independent of the previous cooling program. Moreover, the lower limit $`t_0`$ in the integral can be replaced by $`t_h`$, with a relative error that decreases as $`t`$ increases. Therefore, $`p(t;u)`$ tends to the “normal” solution

$$p_N(T;u)=p_e(T)_{t_h}^t𝑑t^{}\frac{dp_e(T^{})}{dT^{}}\frac{dT^{}}{dt^{}}\chi (t,t^{};u).$$
(96)
Any arbitrary $`p(t;u)`$ approaches $`p_N(T;u)`$ for long enough times, corresponding to high temperatures. The existence of the normal solution for heating processes is not restricted to the hierarchical models considered in this paper, but it has been established for a quite general class of systems, whose dynamics is described in terms of a master equation . In this context, the existence of the normal solution in our simplified, mean-field description of the hierarchical models, provides a consistency test of the approximations we have introduced starting from the master equation formulation. However, it cannot be assured that the special choice of the initial conditions leading to the normal solution in our mean field type approximation remains the same for a more exact analysis.

In the limit of high temperatures, $`p_N(T;u)`$ must be closely related to Hilbert’s expansion solution $`p_H(T;u)`$, analyzed in Sec. IV A. In fact, since $`p_H(T;u)`$ does not refer to any particular initial condition, it is to be expected that the normal solution coincides with Hilbert’s expansion in their common range of validity. However, the normal solution is not restricted to near equilibrium situations, and the system can approach the normal curve in a temperature range for which Hilbert’s solution is not accurate. This point will be clearly illustrated below.

For the sake of simplicity, we will consider that the system is submitted to a cooling process given by Eq. (61), followed by a heating process of the form

$$\frac{dp_e}{dt}=r_hp_e,$$
(97)
where $`r_h>0`$ is the heating rate. Using this heating law, we can express Eq. (95) in the simpler way

$$p(t;u)=\left[p_0(u)p_{e0}\right]\mathrm{exp}\left(\frac{p_e^{au}p_{e0}^{au}}{r_hau}\right)+p_e_{p_{e0}}^{p_e}𝑑p_e^{}\mathrm{exp}\left(\frac{p_e^{au}p_e^{au}}{r_hau}\right).$$
(98)
Here we employ the same notation as in Eq. (62). Similarly,

$$p_N(T;u)=p_e_0^{p_e}𝑑p_e^{}\mathrm{exp}\left(\frac{p_e^{au}p_e^{au}}{r_hau}\right).$$
(99)
The high temperature regime for level $`u`$ is given by

$$\frac{p_e^{au}}{r_hau}1,$$
(100)
and, in this limit, Eq. (99) reduces to

$$p_N(T;u)p_er_hp_e^{1au},$$
(101)
that agrees with Hilbert’s solution. Eq. (48), particularized for the heating law (97).

Figure 6 shows the evolution of $`p(t;u=1)`$, as given by Eq. (98), in a heating process with $`a=1`$ and $`r_h=10^4`$. Two different initial conditions have been considered, corresponding to slow previous coolings of the system with $`r_c=10^4`$ and $`r_c=10^6`$, respectively. In the figure, it is seen that both heating curves approach the normal solution, reaching it in a region where the normal curve represents a clear nonequilibrium state. Hilbert’s expansion approximation, Eq. (101), is also plotted. It provides a good description for the linear correction around equilibrium of the normal curve, but fails for low enough temperatures. In fact, for the values of the parameters considered, Eq. (101) becomes negative for $`T<0.217`$. Therefore, in general, Hilbert’s solution cannot be used to estimate the normal curve over the whole range in which the system is well described by the normal solution, which depends on the details of the previous cooling process. In the figure, for the smallest cooling rate $`r_c=10^6`$, the system reaches the normal curve for a temperature at which Hilbert’s expansion is not valid. This illustrates how the normal solution is relevant for far from equilibrium states, in which a linear theory in the deviations from equilibrium is not accurate.

In the experiments, either with supercooled liquids or with granular materials, the distribution function cannot be directly measured, but instead the average values of the relevant physical properties. In our model, we can consider the dimensionless energy $`\epsilon (t)`$ , which is calculated by integrating $`p(t;u)`$ over $`u`$. For the specific heating law we are considering, it follows from Eq. (98) that

$`\epsilon (t)`$ $`=`$ $`{\displaystyle _0^1}𝑑u\left[p_0(u)p_{e0}\right]\mathrm{exp}\left({\displaystyle \frac{p_e^{au}p_{e0}^{au}}{r_hau}}\right)+\epsilon _e(T)`$ (103)
$`{\displaystyle _0^1}𝑑u{\displaystyle _{p_{e0}}^{p_e}}𝑑p_e^{}\mathrm{exp}\left({\displaystyle \frac{p_e^{au}p_e^{au}}{r_hau}}\right).`$
There is also a normal curve for the energy $`\epsilon (T)`$, that can be obtained by integration of the normal probability distribution, Eq. (99), i. e.,

$$\epsilon _N(t)=\epsilon _e(T)_0^1𝑑u_0^{p_e}𝑑p_e^{}\mathrm{exp}\left(\frac{p_e^{au}p_e^{au}}{r_hau}\right).$$
(104)
For long enough times, independently of the initial probability distribution $`p_0(u)`$, $`\epsilon (t)`$ will approach $`\epsilon _N(T)`$, as a consequence of the tendency of $`p(t;u)`$ towards $`p_N(T;u)`$. The time regime in which $`\epsilon (t)`$ practically agrees with $`\epsilon _N(T)`$ corresponds to the condition given by Eq. (100) being verified for all $`u`$, i. e., when it is fulfilled by the slowest modes, $`u=1`$. Therefore, the approach to the normal curve is controlled by the slowest level $`u=1`$. For very high temperatures, Hilbert’s result (50) holds and, particularizing for the heating law of Eq. (97), we get

$$\epsilon _N(T)\epsilon _H(T)=p_er_hp_e\frac{p_e^a1}{a|\mathrm{ln}p_e|}.$$
(105)
An important feature of the normal curve is that

$$\epsilon _N(T)\epsilon _e(T),$$
(106)
for all temperatures. This property, together with the inequality Eq. (54) for cooling processes, explains the hysteresis effects exhibited by the model when submitted to a thermal cycle. These effects are similar to the experimentally observed behavior in glass-forming liquids  and other complex systems, as granular materials . In Fig. 7 a hysteresis cycle of the energy is shown for $`r_c=r_h=10^4`$. In the cooling process the energy is greater than the corresponding equilibrium value, while in the heating process the tendency of the system to approach the normal curve, which verifies Eq. (106), makes the system overtake the equilibrium curve. Only for very high temperatures, where Hilbert’s expansion is accurate, the normal solution tends to the equilibrium one from below. Comparison of Figs. 6 and 7 shows that the separation of the normal solution for the energy from the equilibrium curve is smaller than that of the slowest level probability $`p(u=1)`$. This is due to the contribution of the other modes, which reach their equilibrium values for lower temperatures.

In our opinion, the explanation of the hysteretic behavior observed in real systems must be similar to the above discussion. This idea is supported by the fact that in other models based on a master equation formalism, analogous results have been found . Also, Hilbert’s expansion, although of limited validity, leads to the conclusion that cooling and heating curves are at opposite sides of the equilibrium curve.

## VII Discussion

In this paper, we have analyzed a simple spin model with hierarchically constrained dynamics. The spins are classified into levels, and a spin in level $`n+1`$ is able to relax only if a certain number of spins in level $`n`$ are in the up (excited) state. Starting from a master equation formulation of the dynamics, a mean-field approximation was introduced. The result is a generalization of Palmer’s et al. model , with a closed equation for the evolution of the probability of finding a given spin in the up state. Each level has a characteristic relaxation rate, which is a function of the temperature of the system. Two sets of parameters characterize the model: the number of spins in level $`n`$, $`N_n`$, and the number $`\mu _n`$ of spins in level $`n`$ involved in the facilitation of the relaxation of a spin in level $`n+1`$. We have chosen the simplest possibility, namely that $`\mu _n`$ and $`N_n`$ are proportional to each other, which leads to a relaxation behavior that is independent of the level populations $`N_n`$. In relaxation processes at constant temperature, the system displays linear logarithmic decay . This is a characteristic feature of the behavior of a wide variety of complex systems, including structural glasses, spin-glasses, protein models and powders.

A key point in our approach to the problem of hierarchically constrained dynamics is the introduction of the “mean-field” approximation, which allows us to reduce the initial, rather involved, problem to a solvable one. This approximation is based on the physical idea that hierarchical constraints render the characteristic time scales of the different levels clearly separated. This is indeed the case in the low temperature region, since the ratio $`\tau _n/\tau _{n1}`$ of the characteristic times of levels $`n`$ and $`n1`$, as given by Eq. (23), diverges for $`T0`$.

The model has also been used to study processes in which the temperature changes in time in an arbitrary way. The general solution for the time dependent distribution function can be explicitly written, and the average energy has a form resembling that of Narayanaswami’s phenomenological theory of glass-forming liquids . By means of Hilbert’s expansion we have constructed an approximate solution for the probability distribution. This expression is valid in the limit of very high temperatures, and for situations where the system is in the linear around equilibrium region. The behavior of Hilbert’s solution is formally identical to the one previously found for a very general class of systems described by master equations , despite the mean-field character of the simplified model considered here. Hilbert’s expansion predicts the existence of hysteresis effects when the system is first cooled down to low temperatures and, afterwards, reheated to high temperatures. This is because the average energy is at opposite sides of the equilibrium curve for cooling and heating processes.

Another point we have addressed is the evolution of the system in continuous cooling processes. A phenomenon similar to the laboratory glass transition shows up, and the system departs from equilibrium at low temperatures. The magnitude of this separation can be measured by the values of the residual properties, which have been analytically computed. A simple but physically appealing argument is presented, in order to understand the origin of this glassy behavior. Each level in the system becomes frozen at the equilibrium value corresponding to a temperature, called the fictive temperature of the level, such that the average number of transitions per spin in that level until reaching $`T=0`$, following the prescribed cooling program, equals unity. A similar argument has been previously used in other models . Here, an analytical derivation is presented. This provides a theoretical basis for the concept of fictive temperature, and clarifies the accuracy of the results following from the qualitative argument. We have compared them with the values of the residual properties obtained numerically as well as with those following from asymptotic analysis calculations.

Finally, heating processes have also been analyzed. More specifically, we have considered heating processes following a continuous cooling of the system down to very low temperatures. In the description of heating processes, the so-called normal curve plays a fundamental role. The hysteresis effects observed when the system is first cooled and afterwards reheated appear because of the trend of the system to approach the normal solution, along the heating evolution. This is analogous to previous results found in models for structural glasses  and for granular systems . Thus, we think it would be worth investigating whether a similar curve does exist for real complex systems showing this kind of hysteretic behavior.

Although a particular simple model with hierarchically constrained dynamics has been considered in this paper, most of the physical ideas developed seem very general. The validity of Hilbert’s expansion for high temperatures, the natural appearance of the concept of fictive temperature, leading a physically appealing description of the laboratory glass transition, the existence of the normal solution and its fundamental role in explaining the hysteresis effects, are results that are not restricted to the present model. In fact, similar results appear in a quite general class of systems described by master equations . Of course, the details of the dynamical behavior depend on the specific model we are dealing with. On the other hand, the general picture developed here, where the heating and cooling experiments in glasses appear as purely kinetic and relatively simple to understand, might not be valid for real and much more complex systems.

###### Acknowledgements.

This research has been partially supported by the Dirección General de Investigación Científica y Técnica (Spain) through Grant No. PB98-1124.


format: text

text:
# Squeezed light from spin squeezed atoms

## Abstract

We propose to produce pulses of strongly squeezed light by Raman scattering of a strong laser pulse on a spin squeezed atomic sample. We prove that the emission is restricted to a single field mode which perfectly inherits the quantum correlations of the atomic system.

Squeezed light and entangled beams of light can be used to probe matter and to study mechanical motion with better resolution than classical light Fabre and Giacobino (1992). Entangled photon sources can be used for lithography with a resolution below the optical wavelength Kok et al. (2000); Bjork et al. (2001), and the active field of quantum information profits from the development of non-classical light sources Mattle et al. (1996); Bouwmeester et al. (1997); Furusawa et al. (1998). Non-linear crystals in optical parameter oscillators and laser diodes with suitable feed-back have been the work horses in most experiments on non-classical light. As a figure of merit for the degree of non-classicality of these light sources, one may refer to the noise suppression observed in direct or homodyne photon detection measurements. Compared to classical sources it has so far been possible to reduce the noise (variance) by about one order of magnitude. In order to make a significant difference in practical applications, further noise reduction is really necessary.

It has been proposed that the large non-linearity and low absorption in resonant Raman systems can lead to ideal four-wave mixing and substantial squeezing Hemmer et al. (1995); Lukin et al. (1998, 1999). As an alternative approach we propose in this Letter to use spin squeezed ensembles of atoms as sources of squeezed light. Atoms can be entangled in such a way that the fluctuations in occupancies of different internal states are significantly suppressed. This phenomenon is refered to as spin squeezing, because a two level atom can be described formally as a spin 1/2 particle, and the interest in spin squeezed states arose already a long time ago in connection with ultra-precise spectroscopy and atomic clocks Wineland et al. (1994). Squeezing of spins was originally believed to be very complicated, but recent proposals based on quantum non-demolition measurements of atomic populations Kuzmich et al. (2000), on coherent interactions in Bose-Einstein condensates Sørensen et al. (2001); Poulsen and Mølmer (2001), and on interactions between laser excited atoms Bouchoule and Mølmer (2001) have changed this impression and suggested that really significant spin squeezing is achievable. The main purpose of the Letter is to demonstrate that the atomic quantum correlations can be perfectly transferred to the field. This result is readily obtained within a simplified model where both atoms and field are described by single harmonic oscillators Jing et al. (2001) but we show that correlations in the atoms can be mapped perfectly on the field also in an a priori multimode situation.

The emission of light is treated by a simple generalization of the theory of stimulated Raman scattering Raymer and Mostowski (1981); Raymer et al. (1985); Carman et al. (1970) to and this part of our proposal can be analyzed without specifying the model for spin squeezing. Our ensemble of two-state atoms is assumed to be strongly elongated and it is treated in a 1D approximation. It is illuminated by a strong laser field $`E_s`$ propagating along the z-axis of the system. This opens up a channel for an atom in the $`b`$ state to go to the $`a`$ state by absorbing a photon of frequency $`\omega _s`$ and emitting a photon of frequency $`\omega _q\omega _s+\omega _{ba}`$. As a result, a field at this frequency builds up and propagates through the sample. This process is described by the following coupled set of equations for the atomic and field operators

$$\begin{array}{cc}\hfill \frac{}{t}& \left[\widehat{\psi }_b^{}(z,t)\widehat{\psi }_a(z,t)\right]=i\kappa _1^{}E_s^{}(z,t)\times \hfill \\ & \times \left(\widehat{\psi }_a^{}(z,t)\widehat{\psi }_a(z,t)\widehat{\psi }_b^{}(z,t)\widehat{\psi }_b(z,t)\right)\widehat{E}_q(z,t)\hfill \end{array}$$
(1)
$$\left(\frac{}{z}+\frac{1}{c}\frac{}{t}\right)\widehat{E}_q(z,t)=i\kappa _2\widehat{\psi }_b^{}(z,t)\widehat{\psi }_a(z,t)E_s(z,t).$$
(2)
where $`\kappa _1=_i\mu _{ai}\mu _{bi}/(\mathrm{}^2\mathrm{\Delta }_i)`$ and $`\kappa _2=2\pi \mathrm{}\omega _q\kappa _1/c`$. $`\mu _{ji}`$ are dipole moments of the atomic transitions and $`\mathrm{\Delta }_i`$ are the (large) detunings with respect to intermediate levels, see Fig. 1. $`\widehat{\psi }_a(z,t),\widehat{\psi }_a^{}(z,t)`$ and $`\widehat{\psi }_b(z,t),\widehat{\psi }_b^{}(z,t)`$ are annihilation and creation operators for atoms in states $`a`$ and $`b`$; $`\widehat{\psi }_b^{}(z,t)\widehat{\psi }_a(z,t)`$ is the positive frequency part of the atomic dipole operator, taking into account the atomic density at position $`z`$ in the ensemble. In Eq. (1) we have assumed that there is no dephasing of the $`ab`$ coherence and we have assumed the validity of the slowly varying envelope approximation for the emitted field in Eq. (2).

We restrict our analysis to the case where the atoms are almost entirely in the $`a`$ state when the $`E_s`$ field is applied. This implies that the population difference appearing in Eq. (1) can be replaced by the density of atoms $`n(z)`$. This density we represent as a c-number throughout the duration of the output coupling. This allows us to define a dipole operator ”per atom” by $`\widehat{\psi }_b^{}(z,t)\widehat{\psi }_a(z,t)=n(z)\widehat{Q}(z,t)`$ and we obtain the linear operator equations

$`{\displaystyle \frac{}{t}}\widehat{Q}(z,t)`$ $`=`$ $`i\kappa _1^{}E_s^{}(z,t)\widehat{E}_q(z,t)`$ (3)
$`\left({\displaystyle \frac{}{z}}+{\displaystyle \frac{1}{c}}{\displaystyle \frac{}{t}}\right)\widehat{E}_q(z,t)`$ $`=`$ $`i\kappa _2n(z)E_s(z,t)\widehat{Q}(z,t).`$ (4)
If we generalize the analysis in Raymer and Mostowski (1981); Raymer et al. (1985); Carman et al. (1970) to inhomogeneous media, we can solve Eqs. (3) and (4) analytically in terms of the input field $`\widehat{E}_q(0,t)`$ at the entrance face of the sample at $`z=0`$ and the initial position dependent atomic polarization $`\widehat{Q}(z,0)`$. In particular, for the field operator we get

$$\begin{array}{c}\widehat{E}_q(z,\tau )=\widehat{E}_q(0,\tau )i\kappa _2E_s(\tau )_0^z\psi _b^{}(z^{},0)\psi _a(z^{},0)\times \hfill \\ \hfill \times J_0\left(2\sqrt{a(\tau )_z^{}^zn(z^{\prime \prime })𝑑z^{\prime \prime }}\right)dz^{}\end{array}$$
(5)
where the new time coordinate is $`\tau tz/c`$, $`J_0()`$ is a Bessel function of the first kind, and $`a(\tau )=\kappa _1^{}\kappa _2_0^\tau \left|E_s(\tau ^{})\right|^2𝑑\tau ^{}`$.

The expression (5) must be evaluated at the position $`z=L`$ of a detector outside the atomic sample. In the absence of atoms, the field equals the incident quantum vacuum field $`\widehat{E}_q(0,\tau )`$. The atomic sample is able to replace the vacuum with an entirely different field. In order to analyze the quantum properties of $`\widehat{E}_q`$ it is convenient to imagine a time integrated homodyne dection at the detector. By choosing the temporal form of the strong local oscillator field in this detection we select a certain spatio-temporal mode of the field represented by the field operator

$$\widehat{a}=\sqrt{\frac{c}{2\pi \mathrm{}\omega _q}}_0^{\mathrm{}}^{}(\tau )\widehat{E}_q(L,\tau )𝑑\tau $$
(6)
where $`||^2𝑑\tau =1`$. In choosing $`(\tau )`$ we should seek to ensure that $`\widehat{a}`$ is a mapping of the precise collective operator of the atomic sample that is known to be squeezed. Mathematically, it is easy to show that in order to probe $`h(z^{})\widehat{Q}(z^{},0)𝑑z^{}`$ we should choose $`(\tau )`$ as the normalized solution of Eqs. (3,4) with the initial condition $`n(z^{})\widehat{Q}(z^{},0)`$ replaced by $`h(z^{})`$. Physically, this reflects that the question of mode matching coincides with the problem of identifying the classical field radiated by a classical dipole distribution. In particular, if the uniform integral of atomic operators $`\psi _b^{}(z^{},0)\psi _a(z^{},0)𝑑z^{}`$ is squeezed, mode matching is accomplished by taking $`(\tau )`$ to be the (normalized) solution to the classical Raman scattering problem. In that case we get the mapping

$$\widehat{a}=\frac{1}{\sqrt{N}}J_{}$$
(7)
where $`\widehat{J}_{}\widehat{\psi }_b^{}(z,t)\widehat{\psi }_a(z,t)𝑑z`$, and $`N`$ is the total number of atoms.

It is a natural concern, whether a 3D analysis will preserve the possibility to couple to only a single spatio-temporal field mode. The issue has been addressed in the classical case, where a diagonalization of the first order coherence function (Karhunen-Loeve transformation) of the field indeed shows that a single mode dominates the output, provided the Fresnel number of the (active part of the) atomic sample is of the order of unity Raymer et al. (1985). Choosing the experimental parameters accordingly we thus expect this result to apply also for the quantum field in three dimensions.

We recall that several schemes now exist for the generation of non-classical states of atomic spins. The above analytical treatment is general and independent of the method of spin squeezing, and as shown by Eq.(7) (valid only if the majority of atoms occupy the state $`a`$), the atomic state is transfered perfectly to the light field. Mean values and variances for the field observables are therefore explicitly known, and, e.g., the orders of magnitude squeezing derived in Sørensen et al. (2001); Poulsen and Mølmer (2001) apply to the emitted light pulse. It may be useful, however, to relate the emitted quantum field directly to the dynamical variables in the spin squeezing process. This may provide further insight in the origin of squeezing; it may provide a useful tool for the application of the squeezed light as input to another quantum system Gardiner (1993); Gardiner and Parkins (1994); and, it may be useful to analyze processes where the spin squeezing and the light emission occurs simultaneously.

As an example, consider spin squeezing by collisional interactions in a two-component Bose Einstein condensate Sørensen et al. (2001); Poulsen and Mølmer (2001): First, a Bose Einstein condensate is formed in only one of the internal states $`a`$. By a short resonant Raman pulse, all atoms are transfered to an equal superposition of $`a`$ and $`b`$. We assume that the interaction strengths between the atoms are not all equal, $`g_{aa}=g_{bb}g_{ab}`$, so that the interaction terms in the fully quantized interaction Hamiltonian

$$\begin{array}{c}\widehat{H}=\hfill \\ \hfill d^3r\{\underset{i=a,b}{}[\widehat{\psi }_i^{}(\stackrel{}{r})\widehat{h}_i\widehat{\psi }(\stackrel{}{r})+\frac{g_{ii}}{2}\widehat{\psi }_i^{}(\stackrel{}{r})\widehat{\psi }_i^{}(\stackrel{}{r})\widehat{\psi }_i(\stackrel{}{r})\widehat{\psi }_i(\stackrel{}{r})]\\ \hfill +g_{ab}\widehat{\psi }_a^{}(\stackrel{}{r})\widehat{\psi }_b^{}(\stackrel{}{r})\widehat{\psi }_b(\stackrel{}{r})\widehat{\psi }_a(\stackrel{}{r})\}\end{array}$$
(8)
cause a Kerr-like phase evolution of amplitudes on states with different numbers of atoms in the two internal states. This results in squeezing of an appropriate collective spin variable as pointed out by Sørensen et al. Sørensen et al. (2001).

As in our earlier work on this problem Poulsen and Mølmer (2001) we use the positive $`P`$ method Gardiner (1991); Steel et al. (1998); Carusotto et al. (2001) to describe the squeezing process. This means that averages of all normal ordered operator products can be calculated as ensemble averages by the substitution of atomic field operators $`\widehat{𝝍}(\widehat{\psi }_a,\widehat{\psi }_a^{},\widehat{\psi }_b,\widehat{\psi }_b^{})`$ by pairs of two-component “wavefunctions”, $`𝝍=(\psi _{a1},\psi _{a2},\psi _{b1},\psi _{b2})`$. The dynamics of the wave functions is given by four coupled and noisy “Gross-Pitaevski” equations, see details in Ref. Poulsen and Mølmer (2001). Starting from a coherent initial state of all atoms in an equal superposition of internal states $`a`$ and $`b`$, we numerically simulate solutions of these equations to obtain an ensemble of $`𝝍`$’s describing exactly (up to sampling errors) the quantum correlations of the system.

Once a sizable spin squeezing is obtained we want to transfer this special quantum state to a light pulse, and to be able to use the results of the above analysis we first bring the internal state of the atoms close to the $`a`$ state, i.e., with a new resonant Raman pulse we rotate the collective spin close to the north pole of the Bloch sphere. This rotation is applied to the individual sets of ”wavefunction” realisations of the simulation ($`i=1,2`$):

$$\left(\begin{array}{c}\psi _{ai}\\ \psi _{bi}\end{array}\right)\left(\begin{array}{c}\psi _{ai}^{}\\ \psi _{bi}^{}\end{array}\right)=\left(\begin{array}{cc}\mathrm{cos}\frac{\theta }{2}& \mathrm{sin}\frac{\theta }{2}\\ \mathrm{sin}\frac{\theta }{2}& \mathrm{cos}\frac{\theta }{2}\end{array}\right)\left(\begin{array}{c}\psi _{ai}\\ \psi _{bi}\end{array}\right)$$
(9)
with $`\theta \pi /2`$.

The atom field interaction in the coupled equations (1) and (2) leads to a natural positive P representation of the field $`\widehat{E}_q(z,t)`$: Replace in Eq. (5) $`\widehat{E}_q(z,t)`$ by a c-number field $`E_{q1}(z,t)`$ and $`\psi _b^{}(z^{},0)\psi _a(z^{},0)`$ by $`\psi _{b2}(z^{},0)\psi _{a1}(z^{},0)`$, and make a similar replacement in the hermitian conjugate equation where $`\psi _a^{}(z^{},0)\psi _b(z^{},0)`$ is replaced by $`\psi _{a2}(z^{},0)\psi _{b1}(z^{},0)`$ to yield the c-number field $`E_{q2}(z,t)`$. In both cases, the incident vacuum field can be represented by a zero, since the positive P representation yields normally ordered expectation values as simple products. From our original ensemble of quadruples of wavefunctions $`𝝍`$ describing the atoms just before we turn on $`E_s`$ we obtain in this way an ensemble of pairs $`(E_{q1}(z,\tau ),E_{q2}(z,\tau ))`$ describing the generated light field $`\widehat{E}_q(z,\tau )`$. With this ensemble any normally ordered field expectation value can be found.

As an example we have investigated 2000 atoms of mass $`m`$ and with 1D interaction strengths $`(g_{aa},g_{ab},g_{bb}=(1.0,0.5,1.0)\times 5\times 10^3\mathrm{}\mathrm{\Omega }l_0`$ where $`\mathrm{\Omega }`$ is the frequency of the harmonic trap and $`l_0=\sqrt{\mathrm{}/m\mathrm{\Omega }}`$ is the associated characteristic length. The atoms are spin squeezed by collisional interactions for a time $`t=3.0\mathrm{\Omega }^1`$. They are subsequently driven towards the state $`a`$, and hereafter they are illuminated with the $`E_s`$ light which builds up a maximum strength of $`E_{max}=10^2\sqrt{2\pi \mathrm{}\omega /l_0}`$ in a time of roughly $`t_{rise}=100l_0/c`$. The matter-light coupling is chosen to be $`\kappa _1=10^3c/2\pi \mathrm{}\omega `$.

As we know the atoms to be spin squeezed we also know one mode of the $`\widehat{E}_q`$ field which will be squeezed: The one corresponding to the simple uniform integral of $`\psi _b^{}(z^{},0)\psi _a(z^{},0)𝑑z^{}`$ as described in the discussion after Eq. (6). Even with the simplifications leading to Eqs. (3) and (4) the atomic state could in principle radiate into many other modes and do so with varying quantum statistics. To check whether such other modes are present in this case we calculate the first order correlation function of the field $`\widehat{E}_q^{}(L,\tau ^{})\widehat{E}_q(L,\tau )=\overline{E_{q2}(L,\tau _2)E_{q1}(L,\tau _1)}`$ ($`\overline{\mathrm{}}`$ indicates averaging over the positive P ensemble). When it is diagonalized we find that almost all population is in fact in the expected mode $`(\tau )`$ which is plotted in Fig. 1. In other scenarios it might be difficult to calculate beforehand exactly which collective atomic operator is squeezed and then an analysis like this would be necessary in order to pick the local oscillator field for the homodyne detection.

Having confirmed that only one mode is populated we now turn to the quantum character of the field. It depends on how we choose $`\theta `$ in Eq. (9): For $`\theta =\pi /2`$ it will approximate squeezed vacuum, for $`\theta `$ slightly different from $`\pi /2`$ it will approximate a squeezed coherent state. As described in Poulsen and Mølmer (2001) the spin squeezing ellipse is at an angle to the coordinate axes. This carries over to the light-field which is squeezed in an appropriate quadrature component $`X_\varphi =(\widehat{a}e^{i\varphi }+\widehat{a}^{}e^{i\varphi })/\sqrt{2}`$. With the above parameters, we find from the positive P simulations that the minimum variance is $`0.04\pm 0.01`$ corresponding to a reduction by a factor of more than 12 from the standard quantum limit. To illustrate the positive P results in Fig. 2 we show both a histogram and a scatter plot obtained from the positive P ensemble of pairs $`(a_1,a_2)`$ representing $`(\widehat{a},\widehat{a}^{})`$ for the mode of the field. $`a_1`$ and $`a_2`$, on average, are complex conjugate quantities, so that the expectation values of hermitian field operators $`X_\varphi =\overline{a_1e^{i\varphi }+a_2e^{i\varphi }}`$ is real for all $`\varphi `$. To represent our simulated results we have made a histogram for the values of $`x=X_0`$ and $`p=X_{\pi /2}`$, indicated by the real part of $`(a_1+a_2)/\sqrt{2}`$ and $`(ia_1ia_2)/\sqrt{2}`$. The histogram in Fig. 2 has been obtained from $`10^5`$ independent realizations of the noisy Gross-Pitaevskii equation for the atomic spin squeezing, while the scatter plot contains only 6800 points each representing a single realization.

The perfect output coupling of atomic correlations is the main result of the paper, and we imagine that many other atomic quantum states may be taken as starting point for non-classical light generation, e.g., EPR-correlated separated atomic ensembles, and that the mechanism may be used also for reliable interspecies teleportation. A squeezed light pulse may also be applied at the ”dark” input port of a beam splitter, to cause the splitting of a strong pulse (with the appropriate mode function $`(t)`$), into two twin pulses with perfectly matched photon statistics. It was recently demonstrated that a classical light pulse can be brought to a complete halt in an atomic sample and that it can be subsequently released Liu et al. (2001). Our study suggests that a quantum light pulse can be similarly stored and retrieved, and that one may perhaps introduce a known or unknown signal field in a condensed sample, manipulate its quantum state by collisional interactions and release it with the original mode function preserved. We emphasize that the spin squeezing mechanism is not crucial for the squeezed light output, and for example the QND-atomic detection and spin squeezing scheme of Kuzmich et al Kuzmich et al. (2000) may be an interesting possibility. Some schemes for atomic spin squeezing may even be compatible with a transmission or production of light with a group velocity so low, that it may constitute a cw source of squeezed light.