text:
# Topological Atom Laser

## I Introduction

One of the most important technological developments to arise to date from dilute gas Bose-Einstein condensation is the atom laser Mewes et al. (1997); Anderson and Kasevich (1998); Hagley et al. (1999); Bloch et al. (1999); Martin et al. (1999). Analogously to its optical counter-part, the atom laser produces highly coherent, directed matter waves, and properties of this output beam, such as temporal and spatial coherence, and beam divergence have been characterized experimentally Köhl et al. (2001); Trippenbach et al. (2000); Coq et al. (2001). The typical atomic laser configuration uses a confining magnetic potential to act as a cavity for the laser mode, which is occupied by a Bose Einstein condensate. Radio-frequency or optical Raman transitions are then used to coherently transfer atoms into untrapped hyperfine states which can propagate freely away from the remaining trapped atoms Ballagh and C.M.Savage (2000). Usually the condensate is taken to be in the ground motional state of its confining potential, though it may possess thermal excitations Choi et al. (2000).

In this paper we consider the matter wave output from a condensate in a topologically excited state. To be definite, we take the condensate to be in a central vortex state, though our formalism is sufficiently general to apply to a broader range of topological states including non-stationary states such as vortex arrays \[e.g., see Madison et al. (2000)\]. We present a $`T=0`$K Gross-Pitaevskii treatment of the topological atom laser dynamics, and give numerical results that characterize the behavior over a wide parameter regime, for two different scenarios: pulsed and continuous scattering of the matter field with the output beam direction transverse to the vortex core, and continuous scattering in the direction parallel to the vortex core. We obtain a linearized analytic solution for the output matter wave, and use it to determine the resonance condition for the Raman process and to characterize the phase properties of the atom laser.

## II Formalism

### II.1 Raman Coupling

In this paper we consider a Raman output coupling mechanism of the type experimentally demonstrated by Hagley *et al.* Hagley et al. (1999). In that scheme the absorption and emission of optical photons from far detuned laser fields coherently transfers atomic population between hyperfine states and imparts sufficient momentum for the output coupled atoms to recoil relative to the trapped condensate.

The ideal situation for an atom laser is to have only two hyperfine states coupled: (1) a weak field seeking state which we label $`|1`$ which is magnetically trapped by the confining potential into the laser modes; (2) an untrapped state $`|2`$ in which atoms can freely propagate to form a directed output matter wave.

The states $`|1`$ and $`|2`$ have an energy difference $`E_2E_1=\mathrm{}\omega _{21}`$, and are coupled through an intermediate state $`|e`$ by two laser fields, with frequency $`\omega _1`$ and wavevector $`𝐤_1`$, and frequency $`\omega _2`$ and wavevector $`𝐤_2`$ respectively. The atom makes a transition from state $`|1`$ to state $`|2`$ by absorbing an $`\omega _1`$ photon, and emitting an $`\omega _2`$ photon, thus transferring kinetic energy

$$\mathrm{}\omega \mathrm{}(\omega _1\omega _2)\mathrm{}\omega _{21}$$
(1)
and momentum $`\mathrm{}𝐤=\mathrm{}(𝐤_1𝐤_2)`$ to the atom’s center of mass motion. The laser fields are assumed to be sufficiently far detuned from the individual transitions that the effects of spontaneous emission can be ignored, and the laser detunings can be represented by a single value $`\mathrm{\Delta }`$. We can then adiabatically eliminate the intermediate state $`|e`$, and represent the laser field coupling between states $`|1`$ and $`|2`$ by a two photon Rabi frequency $`V=\mathrm{\Omega }_1\mathrm{\Omega }_2^{}/2\mathrm{\Delta }`$, where $`\mathrm{\Omega }_j`$ is the single photon Rabi frequency associated with the $`\omega _j`$ radiation field.

The condensate atoms will be scattered into the output state only if the Raman energy transfer is approximately resonant i.e. $`\omega \omega _k`$ where $`\mathrm{}\omega _k=\mathrm{}^2k^2/2m`$ is the recoil energy due to the momentum transfer to the atoms. We note that the quadratic dependence of the recoil energy on $`k`$ makes subsequent Raman transitions into other Zeeman sublevels nonresonant. Thus provided $`V`$ is small compared to $`\omega _k`$, any subsequent coupling will be negligible, and the system can be regarded as being effectively two state.

### II.2 Matter Wave Evolution

The condensate atoms in internal state $`|1`$ have center of mass wavefunction $`\psi _1(𝐫)`$ and are trapped in the cavity resonator mode by a harmonic potential $`V_T(𝐫)`$. Correspondingly, the atoms in state $`|2`$ have wavefunction $`\psi _2(𝐫)`$. The stationary modes of the trap (or resonator) satisfy the Gross-Pitaevskii equation

$$\mathrm{}\mu \psi _1=\frac{\mathrm{}^2}{2m}^2\psi _1+V_T(𝐫)\psi _1+w|\psi _1|^2\psi _1,$$
(2)
where $`\mathrm{}\mu `$ is the chemical potential and the collisional interaction strength $`w`$ is given by $`4\pi \mathrm{}^2a_{11}N_0/m`$ with $`a_{11}`$ the s-wave scattering length for state $`|1`$ and $`N_0`$ the number of condensate atoms.

For a topological atom laser we require that the trapped condensate is in a topological state, whose properties will be transferred to the output matter wave. For simplicity we consider the case where the initial condition for the condensate is a solution of Eq. (2) that possesses a phase winding: in particular we will assume a central vortex state with a single quantum of circulation.

The time evolution of the trapped condensate and the outcoupled matter fields at $`T=0`$K is given by the coupled Gross-Pitaevskii equations Edwards et al. (1999)

$`i\mathrm{}{\displaystyle \frac{\psi _1}{t}}`$ $`=`$ $`{\displaystyle \frac{\mathrm{}^2}{2m}}^2\psi _1+V_T(𝐫)\psi _1+w(|\psi _1|^2+|\psi _2|^2)\psi _1`$ (5)
$`+{\displaystyle \frac{\mathrm{}V(t)^{}}{2}}e^{i(𝐤𝐫\omega t)}\psi _2,`$
$`i\mathrm{}{\displaystyle \frac{\psi _2}{t}}`$ $`=`$ $`{\displaystyle \frac{\mathrm{}^2}{2m}}^2\psi _2+w(|\psi _1|^2+|\psi _2|^2)\psi _2`$ (9)
$`+{\displaystyle \frac{\mathrm{}V(t)}{2}}e^{i(𝐤𝐫\omega t)}\psi _1,`$
We have taken all scattering lengths between and within states to be degenerate, and we note that the normalization condition for the wavefunctions is $`𝑑𝐫\left[|\psi _1|^2+|\psi _2|^2\right]=1`$. In this paper we solve Eqs. (5) and (9) numerically in two and three spatial dimensions and also obtain an analytic solution for the output wave.

## III Linearized Analytic Solution

We will consider the regime where the depletion of the trapped condensate is negligible on the time scale of the output coupling. This condition will allow us to consider a perturbative solution for the output beam amplitude, valid to first order in $`V`$. At this order, we can ignore depletion of the trapped condensate, and write $`\psi _1(𝐫,t)=\sqrt{\rho (𝐫)}\mathrm{exp}\left(iS_0(𝐫)i\mu t\right)`$ where $`\rho (𝐫)`$($`=|\psi _1(𝐫,0)|^2`$) and $`S_0(𝐫)`$ are the density and phase profiles of the initial trapped state respectively. We will also assume that we can neglect the collisional interaction of atoms in state $`|2`$ on themselves, and on atoms in state $`|1`$ . This requires that the density of the scattered atoms is small compared to the trapped atom density in the region of overlap, which is consistent with the negligible depletion assumption. Transforming $`\psi _2`$ to an interaction picture (to remove the free particle-like motion) defined as $`\psi _2(𝐫,t)=\overline{\psi }_2(𝐫,t)\mathrm{exp}(i𝐤𝐫i\omega t)`$, Eq. (9) can now be written in its linearized form as

$`i\mathrm{}{\displaystyle \frac{\overline{\psi }_2}{t}}`$ $`=`$ $`\left[{\displaystyle \frac{\mathrm{}^2}{2m}}^2\mathrm{}\delta i\mathrm{}𝐊+w\rho (𝐫)\right]\overline{\psi }_2`$ (11)
$`+{\displaystyle \frac{\mathrm{}V(t)}{2}}\sqrt{\rho (𝐫)}e^{iS_0(𝐫)i\mu t},`$
where $`\delta =\omega \omega _k`$ is the Raman detuning from free particle resonance, $`𝐊=\mathrm{}𝐤/m`$ is the recoil velocity, and in accordance with the linearized treatment we have ignored nonlinear terms involving $`\overline{\psi }_2`$.

If the velocity spread in the trapped condensate is small compared to $`K`$ then we can ignore the effects of diffusion (about the central momentum $`\mathrm{}𝐤`$) of the scattered wave packet by neglecting the Laplacian term in Eq. (11). Noting that the local velocity of the trapped condensate at position $`𝐑`$ is

$$𝐯(𝐑)=\frac{\mathrm{}_𝐑S_0(𝐑)}{m},$$
(12)
this approximation is valid where $`𝐯(𝐑)^2𝐊^2`$. For typical experimental parameters this condition will be well satisfied everywhere on the trapped condensate except close to the vortex core where the velocity field diverges. Making this approximation we obtain the formal solution for the output coupled wave packet

$`\overline{\psi }_2(𝐫,t)={\displaystyle \frac{i}{2}}{\displaystyle _0^t}𝑑se^{i\mathrm{\Theta }(𝐫,t,s)}V(s)\sqrt{\rho (𝐫+𝐊(st))},`$ (13)
where

$`\mathrm{\Theta }(𝐫,t,s)`$ $`=`$ $`S_0(𝐫+𝐊(st))\mu s+(ts)\delta `$ (15)
$`{\displaystyle \frac{w}{\mathrm{}}}{\displaystyle _s^t}𝑑s^{}\rho (𝐫+𝐊(s^{}t)).`$
In previous work we have presented a similar approximate solution for Bragg scattering (where only a single internal state is involved) and have verified that it provides a good representation of the behavior of the Gross-Pitaevskii equation over a wide parameter regime Blakie and Ballagh (2001).

The interpretation of Eq. (13) is similar to the Bragg case. As the scattered packet moves across the trapped condensate, the amplitude $`\psi _2`$ at a given point (stationary in the frame moving with velocity $`𝐊`$ ) is built up from the sum of contributions coupled in from the trapped condensate at successive points along the trajectory $`𝐑=𝐫+𝐊(st)`$. The contribution coupled into $`\psi _2`$ at time $`s`$ from position $`𝐑`$, has moved to position $`𝐫`$ at time $`t`$ and and has acquired a net phase of $`\mathrm{\Theta }(𝐫,t,s)`$. If the phase term $`\mathrm{\Theta }`$ varies sufficiently slowly along a given trajectory, then the cumulative contributions of matter scattered from the trapped condensate will interfere constructively. Maximal scattering occurs when $`\delta `$ is chosen so that there is a point of stationary phase along the trajectory (i.e. $`d\mathrm{\Theta }(𝐫,t,s)/ds=0`$), which gives the generalized Raman resonance condition

$`\delta \left[_𝐑S_0(𝐑)𝐊\mu +w\rho (𝐑)/\mathrm{}\right]_{𝐑=𝐫+𝐊(st)}.`$ (16)
It should be stressed that Eq. (16) is a local condition, which will apply only at certain positions in the condensate. This leads to the important property that the Raman scattering process can be spatially selective. The resonance condition Eq. (16) differs from the Bragg resonance condition we previously derived Edwards et al. (1999), and gives rise to the following physical interpretation. Immediately after the Raman ejection process, an atom at $`𝐑`$ will have velocity $`𝐊+𝐯(𝐑)`$ \[see Eq. (12)\]. The atom will also have center of mass energy $`\mathrm{}(\omega +\mu ),`$ arising from the chemical potential released from the trapped condensate ($`\mathrm{}\mu `$), together with an additional energy $`\mathrm{}\omega `$ gained from the radiation field. This energy is divided between kinetic energy $`m(𝐊+𝐯(𝐑))^2/2`$ and potential energy $`w\rho (𝐑)`$ arising from interaction with the atoms remaining in the trapped condensate. Energy conservation during the Raman process is therefore expressed as

$$\mathrm{}(\omega +\mu )=\frac{m}{2}(𝐊+𝐯(𝐑))^2+w\rho (𝐑).$$
(17)
One can easily see that Eq. (17) is equivalent to Eq. (16) provided we neglect the term in $`𝐯(𝐑)^2`$, consistent with the validity condition $`𝐯(𝐑)^2𝐊^2`$ used in deriving the analytic solution. It is also convenient below to interpret the term $`S_0𝐊`$ as the Doppler shift due to the local velocity of the condensate. Two particular cases of the resonance condition are worth examining:

* For a non-interacting ground state (i.e. $`S_0`$ $`=0,w=0`$) the resonance condition (16) is $`\delta =\mu .`$ This is in contrast to the result for Bragg scattering for this case, where resonance is at $`\delta =0,`$ with the difference due to the fact that the Raman process ejects the atom from the trapped initial state, so releasing the energy of that state.
* For an interacting ground state (i.e. $`S_0=0`$) the resonance condition is $`\delta =w\rho (𝐑)/\mathrm{}\mu `$. By density weighting the contributions from all $`𝐑`$ across the condensate (assuming a Thomas-Fermi profile) we find that the maximum amount of atoms will be scattered when $`\delta =3\mu /7`$.

## IV Pulsed Output Coupling

Here we consider a configuration of Raman fields which are pulsed in time with a scattering direction (i.e. $`\widehat{𝐊}`$) perpendicular to the vortex core. With this choice, one axis of the plane containing the characteristic spatial and phase variations of the vortex is parallel to the direction of the beam propagation. For convenience we present 2D simulations of this geometry with the Raman momentum difference vector along the $`y`$-direction, i.e. $`𝐤=k\widehat{𝐲}`$. The initial trapped state is assumed to be of the form $`\psi _1(x,y)=R(r)e^{i\varphi }`$, where $`r`$ and $`\varphi `$ are the radial coordinate and polar angle respectively.

In Fig. 2 we show the output matter wave (i.e. $`\overline{\psi }_2`$) after the application of four Raman pulses. The length of each pulse $`\tau _p`$ is chosen to be sufficiently short that the atoms hardly move during the pulse, i.e. $`\tau _pR_c/K`$ (where $`R_c`$ is the condensate radius), while the time between pulses is sufficiently long that the scattered atoms can traverse the trapped condensate before the next pulse is applied. The Raman intensity envelope $`V(t)`$ is shown in Fig. 2(a), and the resulting output coupled density and phase are shown in Figs. 2(b) and (c). It is clear that this procedure produces copies of the trapped vortex state in the output state, including a $`2\pi `$ phase circulation, a result that we confirm analytically below. We note that because $`\overline{\psi }_2`$ is defined in an interaction picture with respect to the free particle center-of-mass motion, the rapid phase gradient due to the mean wave packet velocity $`𝐊`$ has been removed, allowing the phase information in Fig. 2(c) to be more easily discerned. For this reason all the phase plots we show in this paper will be given for the interaction wavefunction.

The Fourier frequency width of each short pulse in our simulation is $`210\omega _T`$, which is much larger than the typical frequencies associated with meanfield and Doppler shifts. In this transient temporal regime, the resonance condition Eq. (16) is inapplicable, however the frequency spread is still sufficiently narrow that there is no significant coupling into other magnetic sublevels. As long as the central frequency is close to resonance ($`\omega \omega _k`$) this temporal-limited frequency spread causes the Raman coupling to be resonant everywhere on the trapped condensate facilitating the complete copying of the wavefunction into the output component. The analytic solution Eq. (13) contains the vortex output result of Fig. 2 as we can show by modeling an individual Raman pulse at time $`t_j`$ as a delta function $`V_0\delta (tt_j)`$. For simplicity we shall choose $`V_0`$ to be real, and then using Eq. (13) we obtain for $`t>t_j`$ (in the non-interaction picture)

$$\psi _2(𝐫,t)=\frac{i}{2}V_0\sqrt{\rho (𝐫𝐑_j)}\mathrm{exp}(i𝐤𝐫)\mathrm{exp}(i\theta )$$
(18)
where $`𝐑_j=𝐊(tt_j),`$ and

$`\theta `$ $`=`$ $`\omega t+S_0(𝐫𝐑_j)\mu t_j+(tt_j)\delta `$ (19)
$`{\displaystyle \frac{w}{\mathrm{}}}{\displaystyle _{t_j}^t}𝑑s^{}\rho (𝐫+𝐊(s^{}t)).`$
The density distribution of the output field in Eq. (18) is the same as the initial vortex state, but displaced by $`𝐑_j`$ (and multiplied by $`V_0/2`$ ). The phase $`\theta `$ of the field can be written

$`\theta `$ $`=`$ $`\left[S_0(𝐫𝐑_j)\mu t_j\right]\omega t_j\omega _k\left(tt_j\right)`$ (20)
$`{\displaystyle \frac{w}{\mathrm{}}}{\displaystyle _{t_j}^t}𝑑s^{}\rho (𝐫+𝐊(s^{}t))`$
which has the following interpretation. The term in square brackets represents the phase of the trapped vortex at time $`t_j`$ , which is transferred to the output vortex along with the phase of the Raman coupling ($`\omega t_j`$) at the time of the pulse. In the absence of collisional interactions the vortex would behave as a free particle of momentum $`𝐤`$, and energy $`\omega _k`$ (we have neglected diffusive kinetic energy) and its phase would increment as $`\omega _k\left(tt_j\right).`$ The integral term in Eq. (20) represents the correction to the free particle evolution during the time the output field is transiting the trapped vortex. Notice that the integrand cuts off once $`|𝐫+𝐊(s^{}t)|>R_c,`$ so that this contribution to the phase becomes constant once the output vortex has separated from the initial vortex. The rotation of the line of zero phase between successive vortices in Fig. 2(c) is thus described by the first 3 terms of Eq. (20).

## V Continuous Output Coupling

In Fig. 3 we show the output matter wave following the application of a long duration, low intensity Raman pulse to a vortex state - which we will refer to here as continuous output coupling. This result displays typical behavior of the scattered solution in the continuous regime: the matter wave is preferentially scattered from a selected spatial region of the trapped condensate and has an asymmetric density distribution transverse to the direction along which it propagates. The shape of the output density profile depends on the detuning $`\delta `$, and by suitably adjusting this parameter, matter can be selectively coupled out from either side of the trapped vortex. Here the frequency width of the finite duration Raman pulse is sufficiently narrow that the generalized resonance condition (16) is applicable. From this resonance condition, the spatial asymmetry of the output density profile can be understood as arising from the Doppler term $`_𝐑S_0(𝐑)𝐊`$. A unit circulation central vortex state with angular momentum number $`m=\pm 1`$ has a phase profile of the form $`S_0=\pm \varphi `$, where $`\varphi `$ is the azimuthal angle. The Doppler effect gives a local shift of the resonant frequency; on the side of the vortex where the current flow runs parallel to the scattering direction the resonant frequency is shifted upward, while on the other side the flow is anti-parallel and the frequency is shifted downward. By choosing $`\delta `$ above or below the central resonance value, it is possible to selectively scatter from either side of the vortex.

To investigate the properties of the continuously output matter wave in more detail we consider the steady state region of the output coupled matter wave. This region, which we indicate pictorially in Fig. 3(a), is where the output beam has constant transverse density profile, that we shall refer to as the *steady state beam profile*. From the formal solution of Eq. (13) it can be shown that for two points $`𝐫`$ and $`𝐫+\alpha 𝐊`$, both within the steady state beam region, then the values of the wavefunction at these points are related by

$$\overline{\psi }_2(𝐫+\alpha 𝐊,t)=\overline{\psi }_2(𝐫,t)e^{i(\delta +\mu )\alpha }.$$
(21)
This result indicates that the steady-state velocity of the output matter wave is $`𝐊+𝐕_0,`$ where we define

$$𝐕_0=\frac{\mathrm{}}{m}\frac{\delta +\mu }{K}\widehat{𝐊}.$$
(22)
Remarkably, this is uniform across the transverse profile, a result we can understand using Eq. (16), which shows that at the position where the Raman transfer takes place, the trapped condensate has a local velocity component parallel to $`𝐊`$ given by

$$𝐯(𝐑)_{}=\frac{\mathrm{}}{m}\frac{\left(\delta +\mu w\rho (𝐑)\right)}{K}\widehat{𝐊}.$$
(23)
The velocity of the atom (in direction $`\widehat{𝐊}`$) immediately after the Raman transfer is therefore $`𝐊+𝐯(𝐑)_{}.`$ We showed in section III that this local ejection velocity gives a kinetic energy consistent with an atom residing in a potential of $`w\rho (𝐑),`$ (due to the trapped condensate atoms). Once the atom has moved away from the condensate, this potential energy is converted to kinetic energy, consistent with an atom of velocity $`𝐊+𝐕_0.`$ We note that the transverse components of the local ejection velocity cancel out in the final scattered wave, due to the symmetry of the vortex velocity distribution. For a general straight-line trajectory across the vortex, there are two positions of local resonance with the same density $`\rho (𝐑),`$ but with opposite transverse local velocities. It is also worth emphasizing that the control available over the output matter wave velocity is in contrast with the case of optical lasers, where the group velocity of the light cannot be altered without changing the properties of medium in which it propagates.

We have numerically verified that Eq. (21) is in good agreement with numerical simulations of the Gross-Pitaevskii Equation over a wide regime. We note that changing the Raman detuning $`\delta `$ to control the phase gradient of the output coupled beam also affects the efficiency at which the atoms are scattered, and to characterize this we show the scattered population, i.e. $`n_2=𝑑𝐫|\overline{\psi }_2(𝐫,t)|^2`$ in Fig. 4. This result shows that in this case the frequency width of the Raman coupling (from $`\delta 30\omega _T`$ to $`30\omega _T`$) is sufficient to allow the output velocity to be changed over a range of $`1.0x_0\omega _T`$ without significant attenuation of the output matter wave. In these results the Doppler effect is the dominant broadening mechanism, so that the approximate width of the Raman transition, $`\mathrm{\Delta }\omega `$, will be

$$\mathrm{\Delta }\omega =\frac{k\mathrm{\Delta }p}{m},$$
(24)
where $`\mathrm{\Delta }p`$ is the momentum width of the condensate in the scattering direction. For the results in Fig. 4 Eq. (24) gives the full width at half maximum $`\mathrm{\Delta }\omega =60.8\omega _T`$, where $`\mathrm{\Delta }p`$ is numerically evaluated from the initial vortex eigenstate (the initial state for $`\psi _1`$).

Additionally, it should be noted that the scattered population ($`n_2`$) is a convenient experimental observable for measuring the linear response of the condensate to the perturbing Raman potential. This technique for probing the condensate is similar to Bragg spectroscopy, which was first used on condensates by Stenger *et al.* Stenger et al. (1999), except that for the Raman case the condensate is scattered into a different internal state. This technique has several advantages over Bragg spectroscopy for conducting high precision measurements of condensate properties which we will explore elsewhere.

Finally we consider continuous output coupling of a central vortex in a direction parallel to the vortex core. Because this arrangement has the superfluid current flow perpendicular to $`𝐊`$, the Doppler term and hence the output coupling is not spatially selective<sup>1</sup><sup>1</sup>1There is a residual resonance shift due to the inhomogeneous meanfield, but we have found that typically this does not give rise to significant spatial selectivity if the condensate chemical potential is of order or smaller in size that the Doppler width (see Eq. (24)).. This means that the vortex profile will be reproduced as the transverse profile of the atom laser output, so that the atom laser will have a similar intensity and phase profile to that of a TEM$`{}_{01}{}^{}{}_{}{}^{}`$ mode of an optical laser, commonly known as the donut mode Siegman (1986).

In Figs. 5 (a) and (b) we show phase and density iso-surfaces of a vortex output-coupled in the manner described above, for realistic experimental parameters. The noticeable beam divergence in the stable beam region, which extends from about $`z=10\mu `$m to $`z=60\mu `$m, is enhanced by the centrifugal forces of the rotating condensate beyond that expected simply from repulsive wave packet spreading

The helicity of the phase iso-surface arises from the addition of the phase gradient along the propagation direction (due to the Raman output coupling) to the vortex-like phase profile the matter wave has in the transverse direction. The phase gradient of the matter wave and hence the phase helicity is controlled by Raman detuning, as given in Eq. (21). We have numerically verified that for a case with the same parameters as in Fig. 5, except using $`\delta \mu `$, the helicity disappears as the phase iso-surfaces are approximately vertical. Also for $`\delta <\mu `$ the helicity changes sign corresponding to the phase iso-surfaces spiraling in the opposite sense.

The helical structure of the phase could be investigated by an interference experiment, such as superimposing the topological atom laser beam with a co-propagating plane phase atom laser.

## VI Conclusion

In this paper we have introduced the idea of a topological atom laser and considered its behavior using numerical and analytic approaches over a broad parameter regime. We have developed a linearized solution for the output coupled matter wave from which we have determined a spatially dependent resonance condition for Raman scattering and characterized how the output beam phase properties relate to the Raman detuning and the condensate chemical potential.

## VII Acknowledgments

This work was supported by the US Office of Naval Research, the Advanced Research and Development Activity. PBB and RJB acknowledge support from the Marsden Fund of New Zealand under contracts PVT902 and PVT202.


format: text

text:
# 1 Introduction

## 1 Introduction

The techniques of Shannon information theory are steadily gaining application in the field of biology, especially bioinformatics. Information theory has been used to study the information content of biomolecules (Gatlin, 1972; Schneider et al., 1986; Pavesi et al., 1997; Adami and Cerf, 2000), as well as the mode and manner in which this information is acquired (Adami et al., 2000, Schneider, 2000), but not as a tool to isolate the selective pressures that adapting populations are subject to. This is surprising as information theory describes not only the difference between entropy and information, but also the transfer and storage of information.

Recently Adami and Cerf (2000) argued that a genomic complexity can be defined rigorously within standard information theory as the information the genome of an organism contains about its environment. Adami et al. (2000) then demonstrated how this complexity must generally increase in simple evolving systems. Here, we focus on the selective pressures of molecular evolution and how they contribute to the evolution of robust and complex structures. Thus, information theory can be a tool that ties together evolutionary biology and molecular biology.

From the point of view of information theory, it is convenient to view Darwinian evolution on a molecular level as a collection of information transmission channels, subject to a number of constraints. In these channels, the organism’s genomes code for the information (a message) to be transmitted from progenitor to offspring, and are subject to noise due to an imperfect replication process. Information theory is concerned with analyzing the properties of such channels, how much information can be transmitted, and how the rate of perfect information transmission of such a channel can be maximized. This information transmission rate can be obtained by multiplying the channel’s symbol transmission rate (the rate at which symbols can be sent from the source to the destination without taking into account the error rate) and its capacity (the maximal information that can be encoded per symbol while being able to guarantee faithful transmission at a given level of noise). A major assumption of the information-theoretic treatment of evolution is that selection acts to maximize the probability that the full genetic information is faithfully transmitted to the next generation. This maximization involves a number of factors that we isolate below, and guides us in elucidating the pressures evolution exerts on the genome itself.

We shall identify three pressures acting on evolving genomes and recast them in an information-theoretic language. This allows us to formulate the forces acting during selection from a unified point of view, while making unique predictions that go beyond the standard population genetics treatment of evolution. The first pressure results in the compression of a message that is not perfectly coded, and (in the absence of any other pressures) tends to reduce message length. The second pressure favors longer messages by increasing the amount of information transmitted per usage of the channel. The third pressure forces populations to code the information in a fault-tolerant manner (when mutation rates are high). While the first two pressures are well-known within standard theory (Bürger, 2000; Baake & Gabriel, 2000, Baake & Wagner, 2001), the third one is novel, and its theoretical treatment has only been elucidated recently (van Nimwegen, 1999; Wilke et al., 2001, Wilke, 2001b, Wilke 2001c). The pressure to code information in a fault-tolerant manner implies that codes should evolve that are robust to deleterious mutations, or, in other words, that the fraction of mutations that are neutral with respect to fitness should increase.

The existence of this neutrality selection has implications for the evolution of evolvability, an effect otherwise known as canalization (Waddington, 1942; Mather, 1953; Gibson & Wagner, 2000). It appears to be active in viral (Burch & Chao, 2000), as well as digital (Wilke et al., 2001) populations, and may even be responsible for the evolution of “robustness genes” such as Hsp90 (Rutherford & Lindquist, 1998; Wagner et al., 1999; Queitsch et al., 2001)

## 2 Selective Pressures

Standard population genetics views the process of adaptive evolution (in the absence of frequency-dependent selection) as the maximization of a suitably defined fitness. Under simple conditions (a population of self-replicating molecules in a single niche, without co-evolution) and ignoring the effects of mutation, a genotype will dominate if its basic Malthusian parameter—the growth-rate of a population seeded by this genotype—is highest. When analyzed from an information-theoretic point of view, we see that the growth rate is in fact affected by three distinct parameters: the length of the message, the maximum rate (in units of bits/sec) at which that message can be transmitted, and the message’s robustness to mutations. We shall take a look at each one of them separately, and demonstrate their effect on an evolving population of digital organisms.

### 2.1 Compression Selection

No matter how information is transmitted, there is a cost for passing along each symbol. This general truth is obvious in biology where transmitting a symbol involves copying a nucleotide, and this cost can be minimized simply by transmitting fewer symbols (a shorter message). This results in a selective pressure that we call compression selection, and illustrate it in a simple experimental system: a population of digital organisms adapting to a virtual world whose complexity and structure is determined by the experimenter (Adami, 1998; Lenski et al, 1999; Adami et al., 2000; Wilke & Adami, 2002). The experiment demonstrating compression selection is identical in spirit to those of Spiegelman (Mills et al., 1967) involving RNA replicating in vitro.

Among a fixed-size population of self-replicators (see Methods), a single limiting resource is allocated evenly to the organisms regardless of phenotype or genotype, such that the only trait selected for is self-replication. In this scenario, compression selection will dominate the other evolutionary pressures on the population. Since it takes time to copy each symbol, there is a pressure to shrink the genome to the minimal length that can still self-replicate. Taking organisms evolved in a complex environment (see Methods) and introducing them into this simple one, we can demonstrate this selective pressure in action. Figure 1(A) shows the adaptive drop in average replication time (CPU time itself being the single resource in the system) as sections of the genome that are meaningless in the simple environment are stripped away. In Figure 1(B) we witness the corresponding drop in genome length. Of course, this is completely analogous to the seven-fold decrease in sequence length observed in the evolution of RNA in a simple in vitro environment (Mills et al., 1967). Naturally, the other experiments described below would be much more difficult (if not impossible) to carry out with RNA.

Although compression selection is fundamental and unavoidable, it is also expected to be weak in systems where the cost associated with transmitting the message on a per-symbol basis is low. However, in simple molecular systems such as the one described above, as well as in periods following the invasion of a new niche or an environmental change, the effects on the genome can be significant due to the removal of genes that have lost their usefulness.

### 2.2 Transmission Selection

Most resources in extant biological systems are spent building “survival machines” as Dawkins (1976) eloquently put it. If an organism obtains an increased share of resources (including expected lifetime), its opportunity for replication is correspondingly increased. This is often considered first by biologists when thinking of interesting examples of adaptation. In terms of information theory, these “survival machines” increase the symbol transmission rate of their channel. Indeed, in a population of self-replicators with a single limiting resource, an organism with twice the affinity to that resource compared to its competitors will produce twice as many offspring, and will soon dominate. In other words, by taking better advantage of its environment, an organism can increase the ”bandwidth” of the channel through which it transmits its genetic message to future generations.

To demonstrate this selective pressure, we modified the previous environment that favored a short sequence length, by allowing the single limiting resource to be actively competed for. For our digital organisms, this meant that the resource (CPU cycles) was now distributed according to phenotype, in this case to those phenotypes that performed selected computations (see Methods). Additionally, in these experiments we scaled the resource allocation with genome length to neutralize the effects of compression selection. This allocation method ensures that organisms with identical computational abilities are equally fit in a mutation-free environment, regardless of genome length. While a non-zero mutation rate creates an additional pressure to reduce the target area for mutations, in a complex environment we nevertheless observe a tendency for organisms to increase sequence length in order to store information about their environment. It is this information that modifies the organisms’ phenotype, allowing them to better exploit the environment (i.e., it increases the capacity of their channels, through which they pass their genomic information to their offspring.)

We performed trials (sets II-IV, see Table 1) starting from a simple, self-replicating ancestor in three environments of differing complexity: simple where no computational tasks are rewarded (as in our compression selection example), medium where the completions of 10 different logical operations are each rewarded by increased resource affinity (in turn, increasing the channel’s symbol transmission rate), and complex where a total of 78 distinct operations are rewarded. In this setup, we expect that transmission selection will cause populations in more complex environments to incorporate more information about those environments into their genomes, resulting in relatively more fit organisms.

Figure 2(A) shows the trajectory of resource affinity divided by replication time (the fitness $`w`$) as function of time. This fitness is a Malthusian parameter measuring rate of replication. Replication time does increase as computational tasks are acquired, but this is more than balanced by the improved resource affinity, resulting in an increase of the overall fitness. Figure 2(B) demonstrates that genome length increases in relation to increasing environmental complexity, as organisms learn to exploit their environment by acquiring information that increases the channel’s symbol transmission rate.

### 2.3 Neutrality Selection

Up to this point, we have been careful to exclude the effects of mutation. Mutations are easily incorporated in an information-theoretical analysis as noise affecting the information transmission channel. This inclusion reveals a third selective pressure not easily apparent at the population level, but affecting the rate at which complexity grows in adapting populations with implications for the evolution of evolvability.

For illustrative purposes, suppose that all mutations are lethal, and as a consequence directly lower the number of expected offspring by the probability of a mutation occurring anywhere in the genome. In this case, the probability of faithfully transmitting genetic information to the next generation, which we have assumed is maximized by evolution, must account for the noise affecting the transmission process. However, just like there are many ways to express a concept in English words, there are many messages (sequences) that contain the same meaning (information), and thus in reality not all mutations are deleterious. If a mutation modifies the genetic sequence without affecting the message, but does make it more robust to further mutations, such a sequence will have an enhanced probability of successful transmission, and therefore possesses a selective advantage. This is most easily thought of as a neutrality selection that acts to increase the probability that mutations are otherwise selectively neutral. While the possibility of such a pressure has been mentioned before (Clarke, 1970; Modiano et al., 1981), its theoretical basis has only been rigorously studied for the case that all mutations are either fully neutral or strongly deleterious (van Nimwegen et al. 1999, Wilke 2001b). Here, we show that this pressure is more general, and exists even when the effects of mutations are distributed on a continuous scale.

As the replication process is not perfect, mutated genotypes often do not contribute to the strength of a species. Denote by $`F`$ (the fidelity of replication) the probability that a birth process gives rise to a genetically identical offspring. If all mutations are deleterious, the true rate of growth (fitness) of a genotype is fidelity times raw replication rate,

$$w_F=Fw,$$
(1)
which is the number of genetically identical offspring produced per unit time. In the language of information theory, this is the number of perfect copies of a message received per unit time.

In reality, of course, not all mutations are detrimental. From an information-theoretic perspective, a sequence whose information content is not altered is as good as an un-mutated copy: the neutrality of a code is akin to a limited error-protection. Therefore, we are interested in the probability that a replication process gives rise to either a true copy, or one that differs from the wild-type by neutral mutations only.

In our system (and macro-molecular evolution in general) the fidelity of the copy process can be written in terms of the error rate per monomer copied, $`R`$, and the length of the genome $`\mathrm{}`$, as $`F=(1R)^{\mathrm{}}`$. Denote by $`\nu `$ the probability that a mutation does not affect the replication rate, $`w`$, of an organism’s offspring. Then, the neutral fidelity

$$F_\nu =(1R+\nu R)^{\mathrm{}}$$
(2)
is the probability of either replicating accurately or else with only neutral mutations. Thus, the combination

$$w_\nu =F_\nu w$$
(3)
is the replication rate of the actual genetic information (in its original or a neutrally modified form). Because $`F_\nu `$ is an increasing function of $`\nu `$, a pressure to increase $`w_\nu `$ implies that there is a selective pressure maximizing $`\nu `$. This is neutrality selection. Note, however, that evolution acts to increase the average $`w_\nu `$ in the population (van Nimwegen et al. 1999; Wilke 2001c). Thus, $`\nu `$-selection implies that selection at high mutation rates acts on groups of mutually neutral genotypes, rather than on the individuals themselves.

Before discussing evidence for neutrality selection in experiments with digital organisms, let us point out that the above formulae can be recast in terms of standard population biology of clonal or randomly mating sexual populations, if $`w`$ is understood as the fitness contribution of a locus or gene.

Equilibrium (meaning mutation-selection balance) demands that the true rate of growth of the wild-type $`Fw_0`$ equals the average growth rate $`\overline{w}`$, or, if there are neutral mutations, $`F_\nu w_0=\overline{w}`$. (Note that the neutrality $`\nu `$ going into the calculation of $`F_\nu `$ is the average neutrality of a population grown from the wild-type (Wilke 2001c), which we approximate in the experiments below by the average neutrality of all one point mutants of the wild-type, see Methods).

We can then calculate the mutational load on the population (the fraction of the population that dies each generation because of deleterious mutations) as

$$L=1\frac{\overline{w}}{w_0}=1F_\nu .$$
(4)
Now, using Eq. (2) and the fact that $`R<<1`$, we obtain

$$L1e^{R\mathrm{}(1\nu )}.$$
(5)
We see that a population can reduce its mutational load, at constant mutation rate, by increasing the neutrality (fraction of all possible one-point mutations that are not deleterious to the organism) $`\nu `$ of its genomes. This has been noted independently by van Nimwegen et al. (1999) and by Wilke (2001b) for the case of strictly neutral evolution.

A decrease of mutation rate (e.g., by the development of error correction) would provide another means of lowering the mutational load. However, a lower mutation rate would also result in reduced genetic variability. According to Fisher’s Fundamental Theorem, a lower genetic variability correlates to lower adaptive potential. As a consequence, neutrality selection may be thought of as an indirect selection for evolvability.

We demonstrate the influence of neutrality selection using Avida by directly measuring the neutrality of the organisms, $`\nu `$, (see Methods) as a function of mutation rate (i.e., for populations adapting under different fixed mutation rate). We examined three sets of trials subjected to different mutation rates (sets V-VII), all in identical complex environments (equalizing the pressure of transmission selection) and constrained to remain at a genome length of 100 instructions (preventing length changes, and hence nullifying compression selection, see Table 1).

Higher mutation rates result in a stronger pressure to increase neutrality. This trend is observed in the measured neutralities of populations adapting in environments at different mutation rates, see Figure 3. Additionally, the interplay between transmission and neutrality selection is apparent from the slow decrease in neutrality for all three sets due to information that increases resource affinity, and is incorporated into the genomes. Since length is fixed, the genome sequence is being “filled up” with information in the form of adaptive mutations that become fixed in the genome, and neutrality must decrease. In other words, as the population gains information about the environment, sites that were once neutral will now cause a loss of information were they to be mutated.

## 3 Methods

Evolutionary experiments were performed on digital organisms (Adami, 1998; Lenski et al., 1999; Adami et al., 2000) implemented with version 1.3 of the Avida software. This software and its documentation can be obtained from http://dllab.caltech.edu/avida. The Avida system controls populations of self-replicating computer programs in a complex and noisy environment within a computer’s memory. Self-replicating programs may be viewed, within limits, as the computational analogs of catalytically active RNA sequences that serve as the template of their own self-reproduction. Populations are propagated for 50,000 updates (an arbitrary unit of time: every update represents the execution of an average of 30 instructions per individual in the population) from the same single ancestor genotype of length 20 instructions (except for set I, see Table 1) to the carrying capacity of 3,600 individuals. Sequence length is variable except for trials V-VII. Mutations occur during replication, at the rate of $`R`$ mutations per instruction copied per organism. Default settings were used in all experiments unless otherwise indicated.

Computations can be carried out by evolved programs if they develop sequences of code (“genes”) that perform logical, bitwise, computations on random numbers provided in their environment. Such genes evolve spontaneously if the performance is rewarded with bonus CPU cycles (their unit of energy). The complexity of the environment can be controlled by changing the number of logical operations whose performance is rewarded. The development of the necessary code to perform such computations is, in the digital world, the analog of the evolution of a sequence or gene enabling the organism to catalyze exothermic reactions leading to faster replication. Outside materials (in this case binary numbers) are taken in and processed for an energy advantage to the organism (see Wilke & Adami, 2002).

From an information theoretic viewpoint, the rewards for these computations directly increase the symbol transmission rate of the channel, because longer messages can be sent per usage of the channel. In our simplest world (set II), no such computations are rewarded, and thus no genes beyond replication develop. In the intermediate world (set III), 10 computations are rewarded, and different combinations of these functions are acquired in the different replicas of that set. In the third and most complex set, 78 different operations are potentially rewarded.

Neutrality was measured by obtaining the fitness of all possible one-point mutants of the most abundant genotype, which allowed us to classify mutations as either deleterious, neutral, or beneficial. Neutral mutations are those for which the fitness has changed less than $`1/N`$, where $`N`$ is the population size (3600 organisms, for all experiments presented here). The neutrality is the fraction of all mutations that are either neutral or beneficial (as these are those that do not decrease the information content).

## 4 Conclusion

In natural systems, neutral mutations often occur on a non-critical nucleotide (synonymous substitutions), or else in code that is never expressed. That most molecular evolutionary change is neutral with respect to the phenotype has been widely accepted since Kimura’s seminal work on the subject (Kimura, 1983). Explaining why so much evolutionary change is neutral has been difficult, and has complicated the discussion of increasing genomic complexity (even though genome length has generally been increasing over the course of evolutionary history, Cavalier-Smith, 1985), since the correlation between genome length and organismic complexity appears weak—the so-called C-value paradox.

Here, we point out that we can think of the quantity maximized by the process of evolution as the probability of faithfully passing on the information encoded in the genome on to the next generation, which is simply the Malthusian parameter ($`w`$) multiplied by the neutral fidelity (2) of the organism. This implies that a population can minimize its mutational load by increasing the probability that a mutation is neutral, or, equivalently, decreasing the probability that a mutation is deleterious. Experiments with digital organisms confirm this view. In particular, an analysis of an aggregate of approximately 10 million generations of the evolution of digital organisms has never revealed a single case in which the effective fitness, Eq. (3), has decreased.

The selection of neutrality has a number of interesting consequences when translated to population biology. First, it might provide a convincing mechanism to explain the actual amount of neutrality in the genomes of higher organisms. Furthermore, because a higher neutrality implies both larger fitness variations and the potential for more information acquisition, neutrality selection implies a selective pressure for the evolution of evolvability (see, e.g., Wagner & Altenberg, 1996; Gerhart & Kirschner, 1997, and Ancel & Fontana, 2000) that is molecular at its origin, and provides added support for the idea that robustness is actively selected for in the evolution of development (Gibson & Wagner, 2000). That it may also have left traces in the genomes of higher organisms in the form of proteins that confer neutrality to other proteins (Rutherford & Lindquist, 1998; Wagner et al., 1999; Queitsch et al., 2000) is a tantalizing suggestion that awaits confirmation.

Acknowledgements This work was supported by the National Science Foundation under grants PHY-9723972 and DEB-9981397. Access to a Beowulf system was provided by the Center for Advanced Computation Research at the California Institute of Technology. Part of this work was performed at the Jet Propulsion Laboratory, California Institute of Technology, under contract with the National Aeronautics and Space Administration. We thank Richard Lenski and Claus Wilke for valuable suggestions.


text:
# Trap Activation Energy and Transport Parameters of HgI2 Crystals for Bubble-Plasma Diagnostics.

111mmiller@jinr.ru

, , ,

## Abstract

In recent data on neutron induced acoustic cavitation in deuterium–containing liquids obtained by neutron measurements it was shown that very high temperatures could arise in some special cases. To study temperature of so–called bubble plasma it is desirable to have micro–detectors of X–rays, which can be prepared on the basis of room–temperature semiconductor detectors, in particular on mercuric iodide ($`\alpha `$–HgI2) crystals. Having in view this aim, the properties of gel–grown ($`\alpha `$–HgI2) crystals was studied by means of isothermal currents, and trap parameters was estimated. Results are promising for special aim of preparing X–ray detectors with moderate energy resolution needed in bubble–plasma diagnostic, though improving of crystal growing technology is necessary222Two authors (M.M.$`\&`$A.S.) are grateful to Russian Foundation for Basic Research for financial support of this work, Pr. #02–02–16397–a..

PACS: 29.40.Wk; 52.70.La

Keywords: X-ray and gamma–ray measurements; semiconductor detectors; mercuric iodide; plasma diagnostics; cavitation

Recent indication on high temperature plasma arising in collapse of cavitation micro–bubbles is of high importance . The data was obtained by neutron measurements in deuterium-containing liquid, and they have only quality character in what regards plasma temperature. For obtaining the value of temperature it is desirable to measure the $`X`$–rays from the bubbles, for which aim it is necessary to have micro–detectors. Such type of instruments can be prepared on the basis of room–temperature semiconductor detectors, in particular on mercuric iodide ($`\alpha `$–HgI2) crystals.

Mercuric iodide ($`\alpha `$–HgI2) is a wide bandgap ($`E_g=2.14`$ MeV) high atomic number ($`Z=62`$) semiconductor, which has attracted considerable interest as a low noise, room temperature $`X`$–ray and gamma ray detector. The detection properties are strongly influenced by energy levels in the forbidden gap arising from chemical impurities and/or structural defects. In the present investigation we carried out measuring isothermal currents on gel–grown HgI2 platelets, in order to estimate trap parameters. The deep energy level $`\mathrm{\Delta }E_t`$ can be determined using the following equation :

$$\mathrm{\Delta }E_t=kT\mathrm{ln}(N_c\sigma _t\nu _{th}\tau _D),$$
(1)
where
$`k`$ – Boltzman’s constant,
$`T`$ – absolute temperature,
$`N_c`$ – effective density of states in the conduction band,
$`\sigma _t`$ – trapping cross section,
$`\nu _{th}`$ – thermal velocity of electrons,
$`\tau _D`$ – detrapping time constant.

Detrapping time and corresponding activation energy of traps are summarized in the Table 1. Moreover, the estimate of mobility–lifetime product of the electrons $`(\mu \tau )_e`$ was performed by measuring the pulse height from alpha–particles 239Pu ($`E_\alpha 5.15`$ MeV). We estimated $`(\mu \tau )_e(1÷5)10^6`$ cm2/V. This value is about one order of magnitude lower than in the best detector grade materials. In future, it would be of interest to study the nature of defects observed and to improve the technology of growing crystals. Taking into consideration not extremely high energy resolution necessary for the aims of diagnostic micro–bubble plasma the above results can be considered as encouraging.


text:
# The iron line in MCG–6-30-15 from XMM-Newton: evidence for gravitational light-bending?

## 1 Introduction

The X-ray spectrum of the bright Seyfert 1 galaxy MCG–6-30-15 ($`z=0.00775`$) has a broad emission feature stretching from below 4 keV to about 7 keV. The shape of this feature, first clearly resolved with ASCA by Tanaka et al. (1995), is skewed and has a peak at about 6.4 keV. This profile is consistent with that predicted from iron fluorescence from an accretion disc inclined at 30 deg and extending down to within about 6 gravitational radii ($`6r_\mathrm{g}=6GM/c^2`$) of a black hole (Fabian et al. 1989; Laor 1991). In part of the ASCA observation the line extended below 4 keV (Iwasawa et al. 1996) which means that the emission originates at radii less than $`6r_\mathrm{g}`$, possibly due to the black hole spinning. XMM-Newton has now observed MCG–6-30-15 twice (Wilms et al. 2001; Fabian et al. 2002a) and in both cases the line extended down to about 3 keV.

The X-ray continuum emission of MCG–6-30-15 is highly variable (see Vaughan, Fabian & Nandra 2002 for a recent analysis). If the observed continuum drives the iron fluorescence then the line flux should respond to variations in the incident continuum on timescales comparable to the light-crossing, or hydrodynamical time of the inner accretion disc (Fabian et al 1989; Stella 1990; Matt & Perola 1992; Reynolds et al 1999). This timescale ($`100M_6`$ s for reflection from within $`10r_\mathrm{g}`$ around a black hole of mass $`10^6M_6`$ $`\mathrm{M}_{}`$) is short enough that a single, long observation can span many light-crossing times. This has motivated observational efforts to find variations in the line flux (e.g. Iwasawa et al. 1996, 1999; Reynolds 2000; Vaughan & Edelson 2001; Shih, Iwasawa & Fabian 2002). These analyses indicated that the iron line in MCG–6-30-15 was indeed variable on timescales of $`10^4`$ s, but that the amplitude of the variations was considerably less than expected and the variations were not correlated with the observed continuum. The data also showed a strong correlation between the photon index of the best-fitting power-law and the flux of the source.

In this letter we examine the later, long XMM-Newton observation in more detail using a simple two-component model to explain the observed spectral variability. Specifically, the data are compared to the model proposed by Shih et al. (2002) comprising a highly variable power-law component plus a constant (or less variable) harder component carrying the iron line. This gives a good fit to the data, with the harder, line-carrying component dominating lowest flux states of the observation. On the assumption that the variable component of the spectrum is just a power-law, the variability enables us to determine the attenuation at low energies due to both Galactic absorption and the warm absorber in MCG–6-30-15. It also enables us to demonstrate that there is no subtle additional absorption influencing the shape of the extensive low-energy “red” wing to the iron line.

## 2 Analysis

Fabian et al. (2002a) discuss details of the observation and basic data reduction. For the present paper the data were reprocessed entirely with the latest software (SAS v5.3.3) but the procedure followed that discussed in the previous paper (except that this time the source extraction radius was 35 arcsec). For both the EPIC pn and MOS data, standard spectral redistribution matrices were used and ancillary responses were generated using ARFGEN v1.48.10. The source spectra were fitted using XSPEC v11.1 (Arnaud 1996) and the quoted errors on the derived model parameters correspond to a 90 per cent confidence level for one interesting parameter (i.e. a $`\mathrm{\Delta }\chi ^2=2.7`$ criterion), unless otherwise stated, and fit parameters (specifically line and edge energies) are quoted for the rest frame of the source.

### 2.1 EPIC spectra

Fig. 1 shows the time-averaged EPIC pn and MOS spectra for the entire dataset compared to a simple power-law model fitted between 2.5–3.0 keV and 7.5–10.0 keV. The emission feature peaking at 6.4 keV is prominent, as is the absorption below 2 keV. However, as is also clear from this plot, the detailed spectral features revealed by the two detectors are subtly different. As the instruments were operated simultaneously these differences represent differences in the spectral calibration between the MOS and pn. In particular, the pn shows a somewhat steeper spectral slope ($`\mathrm{\Gamma }=1.97\pm 0.08`$ for the pn compared to $`\mathrm{\Gamma }=1.86\pm 0.15`$ for the MOS), a larger excess of emission in the 3.5–5.5 keV band, and obvious sharp features around $`1.8`$ keV and $`2.2`$ keV. These last two features are most likely associated with the instrumental edges of Si-K at 1.84 keV and Au-M at 2.3 keV. The two spectra also diverge strongly below $`0.6`$ keV.

Detailed spectral fits to the 2.5–10.0 keV spectrum were presented in Fabian et al. (2002a) based on the MOS data only. These gave a more conservative estimate of the strength and breadth of the broad iron line than the pn data. The greater excess in the pn data in the 3.5–5.5 keV band leads to an apparently stronger line with a significantly stronger red wing. Fitting Model 3111Model 3 consists of a power-law continuum plus a pexrav reflection component ($`R=1`$) and an iron line, both relativistically blurred using a broken power-law emissivity function. of Fabian et al. (2002a) to the 2.5–10 keV pn spectrum gave a good fit ($`\chi ^2=1310.9`$ for $`1492`$ degrees of freedom, $`dof`$) with the following parameters: continuum slope $`\mathrm{\Gamma }=2.02\pm 0.01`$, disc inclination $`i=27.5\pm 0.9`$ deg, inner radius $`r_{\mathrm{in}}=1.95\pm 0.03r_\mathrm{g}`$, break radius $`r_{\mathrm{br}}=6.1_{0.5}^{+0.8}r_\mathrm{g}`$, inner emissivity index $`q_{\mathrm{in}}=5.5\pm 0.3`$, outer emissivity index $`q_{\mathrm{out}}=2.7\pm 0.1`$ and a line equivalent with of $`EW=685_{32}^{+50}`$ eV. The equivalent width of the possible $`6.9`$ keV emission line is $`EW=17\pm 3`$ eV.

It should be noted that while the above calibration problems remain with the EPIC spectral responses, it is not clear whether the MOS or the pn data give a more accurate description of the shape of the broad iron line. However, while the absolute values of best-fitting model parameters are sensitive to detector-dependent calibration, the spectral variability characteristics should be, to first order at least, independent of these problems. Therefore, the pn data are used for the purposes of the spectral variability analysis below as they have the highest signal-to-noise ratio. It should be noted that while the form of the spectral variability should be calibration-independent, the exact details of any spectral models will remain sensitive to the calibration uncertainties outlined above.

### 2.2 Spectral variability

#### 2.2.1 The difference spectrum

In this section the two-component model discussed in Shih et al. (2002) is applied to spectra obtained from the XMM-Newton observation. If the spectrum can be broken into two distinct emission components then the observed spectrum during any given time interval, accounting for the effects of absorption, can be written as $`S=A\times (N_1s_1+N_2s_2)`$, where $`s_1`$ and $`s_2`$ represent the spectra of the two emission components, $`N_1`$ and $`N_2`$ are their normalisations, respectively, and $`A`$ gives the absorption profile.

Shih et al. (2002) found that such a model adequately described the spectral variability observed in MCG–6-30-15 by ASCA, with one component having a constant flux while the other varied only in normalisation. If the spectral variability can indeed be described in terms of a variable emission component and another component that remains constant (or, at least, varies very little on the timescales of interest) then it becomes possible to examine the absorption profile using the spectral variability. The difference between two spectra that differ significantly in flux (different $`N_1`$ but identical $`N_2`$) is simply $`A\times (N_{1,\mathrm{high}}N_{1,\mathrm{low}})s_1`$. The difference spectrum is thus the spectrum of the variable component only, modified by absorption. If the variable spectral component is a power-law, as suggested by Shih et al. (2002), then the difference between two spectra of different fluxes should be a simple power-law modified by (Galactic and intrinsic warm) absorption.

In order to examine this possibility the XMM-Newton observation was divided into $`32\times 10`$ ks intervals (Fig. 2) from each of which spectra were obtained. The lowest-flux spectrum is from interval 303:k. A high-flux spectrum was extracted from the first 80 ksec of data from this revolution (303:a–h). The average count rate during the high-flux interval is $`32.6`$ ct s-1 and during the low-flux interval the count rate is $`13.5`$ ct s-1 (after correcting for the live-time of the pn camera in small-window mode). Subtracting the two spectra (i.e. high$``$low), accounting for the difference in exposure times, gives the difference spectrum shown in Fig. 3.

#### 2.2.2 A simple absorption model

Across the 3–10 keV band the difference spectrum is well fitted by a power-law, with photon index $`\mathrm{\Gamma }=2.11\pm 0.05`$. The strong, broad iron line, while clear in both the high- and low-flux spectra when examined separately, is not apparent in the difference spectrum. The iron line flux must be very similar in these two spectra, even though the total (i.e. continuum) flux differs by a factor of about 2 above 7 keV. This then implies a very strong line equivalent width during the low-flux interval.

Extrapolating the power-law model to 0.5 keV reveals the attenuation spectrum of the source (Fig. 3). The warm absorber is clear as the large drop above $`0.7`$ keV which recovers by about $`2`$ keV, and the Galactic absorption mostly accounts for the decrease at lower energies. These features can be accounted for using a very simple absorption model (similar to that used to model ASCA data of the source; Otani et al 1996) comprising a power-law modified by absorption from neutral gas of column density $`N_\mathrm{H}=3.4\pm 0.4\times 10^{20}`$ cm-2, explained by the Galactic column density ($`4.06\times 10^{20}`$ cm-2; Elvis, Lockman & Wilkes 1989), together with two warm absorption edges (at $`E_1=0.72\pm 0.005`$ keV of depth $`\tau _2=0.61\pm 0.03`$ and at $`E_2=0.86\pm 0.01`$ keV with depth $`\tau _2=0.19\pm 0.04`$). The obvious candidates for such edges are OVII and OVIII, respectively.

We note that the above simple absorption model may be adequate for fitting data at CCD resolution above 1 keV but cannot be the whole explanation for data below that energy, where absorption by lines and Unresolved Transition Arrays are expected. The important point for the present work is that there are no major absorption features associated with iron-K, i.e. above 7 keV, nor with Si, S etc between 2–4 keV. The broad iron emission line is not strongly affected by any absorption features. Very similar results are obtained if the absolute minimum, 11 ks, spectrum from within intervals 303:k and l is used. Note that these results depend on the absorption not varying, which is consistent with the rms variability spectrum (Fabian et al 2002a) showing no edges.

#### 2.2.3 The low-flux spectrum

The absorption model derived above was then applied to the spectrum obtained from the low-flux interval. Fitting this spectrum with a power-law, and including the absorption as described above, reveals the very strong iron emission line and an excess of emission below $`1.5`$ keV (Fig. 4).

The iron line can be modelled by a LAOR profile (Laor 1991) with a high equivalent width ($`EW=1.4`$ keV) extending into an inner radius of $`1.8r_\mathrm{g}`$ with an emissivity profile index of $`3.3\pm 0.3`$. The spectrum can also be fitted reasonably well with the constant density ionised disc model of Ross, Young & Fabian (1999), when relativistic blurring is included (again using a broken emissivity law). The model used included three times solar abundances and the best-fitting ionisation parameter was $`\xi =10`$, with the reflection component dominating over the power-law ($`R=5.5`$). We hereafter refer to this spectral model as the Reflection-Dominated Component (RDC) .

#### 2.2.4 Fitting the whole light curve

We now extend the study of the simple, absorbed, two-component model to the whole dataset. A model consisting of the simple absorption model from Section 2.2.2 applied to the best-fitting RDC and a Power-Law Component (PLC), with all parameters apart from the normalizations and PLC photon index fixed, was then fitted to all 32 spectra over the energy band 1–10 keV.

The results are shown in Fig. 5. All fits gave acceptable values of $`\chi ^2`$. The photon index $`\mathrm{\Gamma }`$ of the PLC is seen to lie between 2 and 2.3 for most of the spectra. The deviant results occur when the PLC normalization is low and disappear if the fit is made over the 2–10 keV band. It is possible that there are small changes in the absorption when the flux is low. The RDC normalization is approximately constant and shows variations of up to $`\pm 25`$ per cent, with the variations showing no obvious correlation with other parameters.

## 3 Interpretation

We have shown that the spectral variability of MCG–6-30-15, on timescales of 10 ks, is accounted for by an almost constant reflection-dominated component together with a highly variable power-law component, which has an almost constant photon index, supporting the model of Shih et al (2001). We now attempt to interpret the source using a two component RDC plus PLC model. An alternative interpretation, in which the line is progressively ionized and weakened as the source brightens, is presented and discussed by Ballantyne, Vaughan & Fabian (2002).

Explaining the constant emission component provides us with a significant challenge. It appears to be mostly due to reflection. However it is not simply due to reflection of the observed power-law component since that repeatedly varies by factors of two or more on short timescales (but not shorter than the inferred emission region size, if it is indeed only a few $`r_\mathrm{g}`$). In other words the RDC and PLC appear to be separate. Since however both show the effects of the warm absorber they must originate in a similar location to within the dimensions of the warm absorber (probably less than 1 pc). As the iron line in the RDC indicates emission peaking at only a few gravitational radii we shall henceforth assume that this is indeed where that component originates. We then have to explain why so little of the continuum which causes the line is observed and why the power-law continuum we do see does not seem to have any associated reflection.

We are driven to consider that the intrinsic X-ray emission is anisotropic, so that the reflection component is a smoothed average of many variable power-law components, most of which do not illuminate our line of sight. This could be due to source geometry, relativistic beaming of the emitters, or to the general relativistic bending of light expected from an emission region so close to the black hole. We now address each of these in turn.

### 3.1 Geometrical effects

As discussed by Fabian et al. (2002b) in the case of 1H 0707–495, the surface of the disc could be corrugated. If the continuum is emitted close to the surface within the corrugations then it could show large apparent variability as part of it disappears within a corrugation. The reflected component could be strong, since we are only seeing part of the intrinsic PLC (and multiple reflection may operate; Ross et al. 2002), and be approximately constant.

### 3.2 Intrinsic beaming

If the emitting region, or the emitting electrons, are relativistic, then the emission can be beamed. We may see variations as beams from different regions or populations of electrons sweep our line of sight. Reflection from the disc however responds to the average of many beams sweeping over it. Note that if the particles are relativistic then the reflection can be strong due to the intrinsic anisotropy of the inverse-Compton process (Ghisellini et al. 1991).

### 3.3 General relativistic effects

If we accept the evidence from the line shape the the emission peaks at say $`2r_\mathrm{g}`$, then the returning radiation (Cunninham 1975) will be very important (see also Martocchia et al. 2002; Dabrowski et al. 1997). An observer on the disc at that radius would see much of the Sky covered by disc, due to gravitational light bending. The reflection is therefore strong and dominated by the average PLC emission from the opposite side of the disc. The PLC which we see will be dominated by the approaching side, perhaps associated with the plunge region just within the marginally stable orbit (Krolik 1999, Agol & Krolik 2000). This will be a much smaller region, also Doppler boosted, than that to which the reflection responds. If its height varies with time then so will the flux of the PLC, whereas the RDC will see a much more constant illumination and so vary little (see Dabrowski & Lasenby 2001, Fig. 12, for the expected variation in equivalent width, dominated by changes in height of the PLC). The appearance of the PLC and RDC are therefore essentially decoupled, as required.

## 4 Discussion

The iron line in MCG–6-30-15 appears to indicate an extreme disc around a spinning black hole. We have proposed several explanations for the apparent constancy of the RDC and distinguishing between them is not simple. General relativistic effects ought however to be important if much of the emission originates at about $`2r_\mathrm{g}`$. They do provide a simple explanation in which small changes in the height of the power-law emission region lead to large changes in the observed flux as a consequence of gravitational light bending; much of the flux is bent onto the disk giving a strong, almost constant, reflection component. This model will be pursued in detail elsewhere.

If instead the emission is anisotropic due to relativistic beaming, then the PLC might extend as a power-law to high energies, perhaps above 500 keV as seen in the hard tail in the soft state of Cygnus X-1.

In summary, we find that the spectrum of MCG–6-30-15 is well represented by the sum of an almost constant reflection-dominated component and a highly variable power-law component. We emphasise that this is only an approximation to the spectral variability as it does not account for the observed changes in the power-spectral density with energy, nor the hard lags (Vaughan et al 2002). Acting on both is a warm absorber, which above 1 keV is relatively simple and dominated by highly-ionized oxygen. Provided that the emission predominantly originates from about $`2r_\mathrm{g}`$, then the strengths of the components and their behaviour are accounted from by strong gravtiational light-bending. The black hole in this case is presumably spinning rapidly.

Strong broad iron lines in active galactic nuclei may require the central black hole to have a high spin in order that returning radiation creates a strong reflection component. If the rapid variability in MCG–6-30-15 is in part due to strong light-bending amplifying small changes in source height, then the similarity with the timing behaviour of the soft state of Cygnus X-1 (Vaughan, Fabian & Nandra 2002) suggests that it too may be rapidly spinning. Further extrapolation to the variability of many Narrow-Line Seyfert 1 galaxies (e.g. Boller et al 1997) is also interesting.

## Acknowledgements

Based on observations obtained with XMM-Newton, an ESA science mission with instruments and contributions directly funded by ESA Member States and the USA (NASA). We thank David Ballantyne, Russell Goyder, Kazushi Iwasawa, Anthony Lasenby and Andy Young for discussions. ACF thanks the Royal Society for support.


text:
# The order parameter–entropy relation in some universal classes: experimental evidence

## I Introduction

The observation of fluctuation-driven critical phenomena in structural phase transitions is subject of long and persistent controversies. Although power law behaviours near phase transitions can be “universal”, the temperature range where this asymptotic behaviour can be observed within experimental resolution (the so-called asymptotic or critical region) is not. It depends in general on the particular features of each system and can even vary for different physical quantities. It is no wonder, therefore, that the size of the critical region in structural phase transitions has been a matter of discussion for several decades. In general, one should expect it to be smaller than in the pure order-disorder systems, the spin hamiltonians, that represent the corresponding universality classes. However, the Ginzburg criterion has been equally invoked to predict typical critical regions of size too small to be of any experimental significance,Ginzburg et al. (1987) as to justify critical regions of the order of 0.1 in reduced temperature ($`T/T_c`$)Müller and Berlinger (1971) (an interval that would greatly surpass that of the actual spin-hamiltonian representatives of some of the universality classes). In addition, structural defects may produce strong deviations from the universal asymptotic laws predicted for ideally perfect crystals.Levanyuk et al. (1979); Scott (1981); Marro et al. (1986)

Authors expecting insignificantly small critical regions have argued repeatedly that the experimental observation of power laws around the phase transition with non-mean-field exponents may have a rather simple explanation, e.g. in terms of the Landau theory when sixth-order or higher order terms are taken into account in the Landau expansionGinzburg et al. (1987). For instance, a 2-4-6 Landau potential:

$$\mathrm{\Delta }G=\frac{1}{2}AϵQ^2+\frac{1}{4}BQ^4+\frac{1}{6}CQ^6$$
(1)
with temperature independent parameters $`A,B,C`$ and $`ϵ=1T/T_c`$ yields indeed a mean field critical behaviour for the order parameter ($`Qϵ^{1/2}`$) with a mean-field exponent $`\beta =1/2`$, but the asymptotic behaviour is limited to $`ϵB^2/4AC`$. This temperature interval can be beyond any observation if the fourth order coefficient is small enough, a typical situation happening if the system is close to a tricritical point. In $`ϵ`$-intervals of the order of 0.01-0.1, or even larger, the behaviour of the order parameter can then be fitted to power laws with “effective” exponents which are intermediate between the mean field value 1/2 and the one corresponding to mean field tricritical behaviour (1/4).

More generally, it became clear that any experimental observation of expected asymptotic critical behaviour for the order parameter requires further confirmation through cross-examination of its consistency with the behaviour of other physical parameters near the transition point. The correlation of the order parameter with the transition excess entropy has been recently used for this purposeSalje et al. (1998). The excess entropy $`\mathrm{\Delta }s`$ is expected to be proportional to the square of the order parameter $`Q`$ not only in the mean field asymptotic regime ($`\beta =1/2,\alpha =0`$), but also for any behaviour described by a Landau potential like equation (1). The relation would be true as long as the terms higher than quadratic in the Landau potential are temperature independent and the quadratic term has a linear temperature dependence. This particular relation between entropy and order parameter is in principle broken in the asymptotic regime for all other universality classes. Its experimental check can be then an effective way to distinguish between genuine fluctuation-driven critical behaviour and pure Landau behaviour. Indeed, this proportionality relation has recently been shown for SrTiO3, in contradiction with a common belief that this compound exhibits in a significant temperature interval asymptotic critical behaviour corresponding to the Heisenberg–3D or the Ising–3D universality classSalje et al. (1998); Chrosch and Salje (1998); Hayward and Salje (1999). The analysis of $`\mathrm{\Delta }s`$, instead of checking directly the expected critical behaviour of the excess specific heat $`\mathrm{\Delta }c`$, has the fundamental advantage of using an integrated quantity with a smooth temperature dependence that is independent of experimental statistical errors in the measurement of $`\mathrm{\Delta }c`$.

In the present paper, the approach of comparing $`\mathrm{\Delta }s`$ and $`Q^2`$ for analysing the thermal behaviour at phase transitions is further investigated and applied to several materials. The expected behaviour of $`\mathrm{\Delta }s/Q^2`$ for the different universality classes and the Landau case, is discussed and compared with that derived from experimental data for Rb2CoF4 (antiferromagnetic phase transition), Rb2ZnCl4 (normal-incommensurate phase transition) and SrTiO3 (antiferrodistortive transition).

## II  Entropy-order parameter relation in the critical region

The asymptotic entropy-order parameter relation for the different universality classes follows from the critical laws for specific heat and order parameter. It can be expressed as

$$\frac{\mathrm{\Delta }s}{Q^2}ϵ^\kappa \mathrm{\Delta }sQ^{2\kappa ^{}}$$
(2)
with

$$\kappa =\{\begin{array}{cc}1\alpha 2\beta \hfill & \alpha >0\hfill \\ 12\beta \hfill & \alpha 0\hfill \end{array}$$
(3)
and

$$\kappa ^{}=\{\begin{array}{cc}\frac{1\alpha }{2\beta }\hfill & \alpha >0\hfill \\ \frac{1}{2\beta }\hfill & \alpha 0\hfill \end{array}$$
(4)
In order to get these expressions the normal asymptotic behaviour of the order parameter, characterized by the exponent $`\beta `$, is assumed while that of the excess specific heat is supposed to be:

$$\mathrm{\Delta }c\frac{ϵ^\alpha 1}{\alpha }$$
(5)
Note that the specific heat is not reduced to a pure power law as it includes the cases of negative $`\alpha `$ values and the logarithmic divergence for $`\alpha 0`$.Heller (1967); Griffiths (1967)

The resulting exponents $`\kappa `$ and $`\kappa ^{}`$ for the different universality classes and the Landau behaviour are listed in table 1. Non-classical exponents $`\kappa `$ and $`\kappa ^{}`$ differ significantly from the mean-field/Landau trivial values associated with the above mentioned proportionality law. Nevertheless, these deviations may be rather difficult to detect in a typical $`\mathrm{\Delta }s`$ vs. $`Q^2`$ plot. In figure 1, we present these plots for the different universality classes assuming that, in an $`ϵ`$ interval of 0.1 and within experimental resolution, both the order parameter and the excess specific heat follow their asymptotic critical behaviour (see equation (5) for that of the excess specific heat).

The excess entropy in figure 1 has been calculated by numerical integration of the expression:

$$\mathrm{\Delta }s=\underset{T_c}{\overset{T}{}}\frac{\mathrm{\Delta }c}{T^{}}dT^{}=\underset{0}{\overset{ϵ}{}}\frac{\mathrm{\Delta }c}{1ϵ^{}}dϵ^{}$$
(6)
this integration leads generally to smooth curves $`\mathrm{\Delta }s(T)`$ even if $`\mathrm{\Delta }c`$ contains significant statistical noise.

Classical Landau behaviour is represented in figure 1 by a straight line through the origin. The Ising–3D, Heisenberg–3D and XY–3D classes yield curves different but close to a straight line. Therefore, in practice, the deviation from the Landau behaviour is difficult to assess quantitatively through these plots. Only in the case of the Ising–2D class the curve, being a fourth order parabola, clearly deviates from the proportionality law. The different relationships between $`\mathrm{\Delta }s`$ and $`Q^2`$ become more obvious in a log-log plot of $`\mathrm{\Delta }s/Q^2`$ vs $`ϵ`$. The plots corresponding to the curves of figure 1 are presented in figure 2. The deviation from the horizontal Landau behaviour is then obvious for all universality classes. The larger deviation is seen for the Ising–2D case, while all other universality classes lead again to more similar curves. It also becomes evident in the figure that although the curves have been obtained assuming perfect asymptotic behaviour for $`\mathrm{\Delta }c`$ and $`Q`$ up to $`ϵ=0.1`$, the magnitude $`\mathrm{\Delta }s/Q^2`$ is far from being in the asymptotic regime represented by the power law of equation (3) even at $`ϵ`$ values so low as 0.0005. The curves deviate rather smoothly from straight lines and in fact their approximate slopes are far from the asymptotic values given by equation (3) and listed in table 1. Indeed, the “critical” region for the excess entropy is much smaller than for $`\mathrm{\Delta }c`$, since the integration represented by equation (6) yields the sum of several competing terms with similar but different powers of $`ϵ`$. Hence, according to figure 2, the asymptotic behaviour for $`\mathrm{\Delta }s/Q^2`$ represented by equation (2) is not experimentally accessible in any case. Nevertheless, the curves in figure 2 can be considered “universal” in the sense that they should be satisfied for any system so long as its order parameter and specific heat are following the asymptotic critical behaviour of the corresponding universality class. In the $`ϵ`$ interval \[0.001-0.01\] the curves can be approximated to straight lines, their slopes correspond then to “universal” effective exponents $`\kappa _e`$ for the approximate law $`\mathrm{\Delta }s/Q^2ϵ^{\kappa _{e_T}}`$ restricted to this temperature interval. The exponents $`\kappa _{e_T}`$ obtained from the linear fit of the curves in figure 2 in the mentioned temperature interval are listed in table 1.

## III Experimental cases

We consider now under this approach three very different systems: strontium titanate (SrTiO3), rubidium tetrafluorocobaltate (Rb2CoF4) and rubidium tetrachlorozincate (Rb2ZnCl4).

### III.1 SrTiO3

The antiferrodistortive phase transition of SrTiO3 has been widely considered one of the few examples of a structural phase transition with a clear observable critical region. Müller *et al.*Müller and Berlinger (1971) could fit their EPR data assuming a power law for the order parameter with $`\beta =0.35`$ in an $`ϵ`$ interval of the order of 0.1. This behaviour was then identified as the expected asymptotic behaviour of the Heisenberg–3D class. Recently, the excess entropy was determined and compared with the order parameter data; a mean field linear correlation was observedSalje et al. (1998), thus the critical region might be beyond observation. For comparison with earlier specific heat data see Ref. Gallardo et al. (2002).

For the sake of comparison with the two following systems, the $`\mathrm{\Delta }s`$-$`Q^2`$ relation obtained in Ref. Salje et al. (1998) for SrTiO3 is reproduced in figure 3 in a format similar to next figures.

### III.2 Rb2CoF4

This compound exhibits an antiferromagnetic phase transition at $`107\mathrm{K}`$. This system is believed to belong to the Ising-2d universality class. The structure of Rb2CoF4 is isomorphic to that of K2NiF4 which is a well known two-dimensional anti-ferromagnetic system. The two-dimensional character is related to the presence of coupled anti-ferromagnetic planes separated by planes with weak interplanar effective interactions. Specific heat Ikeda et al. (1976) and the reduced sublattice magnetization (order parameter) obtained from neutron scattering experimentsSamuelsen (1973) have been reported to be consistent with the critical behaviour predicted for this universality class. In figure 4 we present the correlation of the excess entropy entropy $`\mathrm{\Delta }s`$ and the square of the order parameter (the reduced sublattice magnetization $`\sigma `$) for this material. The entropy has been calculated upon integration of the specific heat measured by Ikeda *et al.*Ikeda et al. (1976). The magnetization was determined from neutron scattering experiments by SamuelsenSamuelsen (1973). The figure shows clearly how the entropy fails to be proportional to the square of the order parameter.

### III.3 Rb2ZnCl4

Rb2ZnCl4 presents an orthorhombic (Pnma) phaseCummins (1990) above $`T_c305\mathrm{K}`$. At this temperature it transforms through a continuous phase transition into an incommensurately modulated crystal with a modulation wave vector close to one third of $`𝐚^{}.`$ Further below, at $`195\mathrm{K}`$ the modulation wave vector locks into the exact value $`1/3𝐚^{}`$, and the incommensurate phase transforms into an orthorrhombic commensurate ferroelectric phase with a three-fold unit cell along the modulation. This compound belongs to the family of type A2BX4, where this incommensurate instability is rather common, and has been extensively studied in order to test the theoretical predictions on incommensurate phase transitions. In 1978 Bruce and CowleyBruce and Cowley (1978) argued that this transition should belong to the XY–3D universality class. Since then, many workers have related their results to what is expected in the critical zone for this universality class. For instance, Walisch *et al.*Walisch et al. (1989) accurately determined from quadrupole perturbed NMR experiments the critical exponent of the order parameter with full agreement with the theoretical value for the XY–3D class. According to these results the critical zone would reach up to $`ϵ`$ values of the order of 0.07. Diffraction synchrotron measurementsZinkin et al. (1996) confirmed these results and a critical region of the order of 0.1. Specific heat was also reportedChen (1990) to agree expected critical behaviour in a range up to $`ϵ0.2`$. More recently, Haga *et al.*Haga et al. (1995) made a more comprehensive study of the specific heat and fitted their data to an expression which includes corrections up to second order of the asymptotic law. The fitted $`ϵ`$ interval reached in this case up to $`ϵ=0.1`$. The existence of observable large critical regions in incommensurate phase transitions, in comparison with normal structural phase transitions, was predicted and justified.Patashinskii and Prokovskii (1974); Prokovskii (1979) However, other authors maintain that there is no theoretical ground for expecting in incommensurate systems larger critical regions than in other structural phase transitionsGinzburg et al. (1987); Sandler and Aleksandrov (1983) and have interpreted for instance other NMR experiments in this compounds within a Landau approachFajdiga et al. (1992).

For the present study, a new independent measurement of the specific heat of Rb2ZnCl4 was performed in a conduction calorimeter. The calorimeter has been previously described del Cerro et al. (1987) and has been used successfully for the measurements of thermal properties of single-crystalsMartín et al. (1997). A single crystal of $`1.078\mathrm{g}`$ ($`2.850\mathrm{mmol}`$) was measured. The sample was grown from aqueous solution and was annealed at the high symmetry phase and then cooled down at a low rate of $`0.4\mathrm{Kh}^1`$ while the heat capacity of the sample was recorded. The heat capacity agrees previously reported dataHaga et al. (1995); López-Echarri et al. (1980); Atake et al. (1983) and they show a $`\lambda `$-shaped peak at about $`304\mathrm{K}`$.

The background contribution to the specific heat was determined by fitting a smooth polynomial function to data above $`T_c`$ and data well below the lock-in phase transition located at $`T_i=195\mathrm{K}`$. By subtracting the background contribution, the excess specific heat was determined and the excess entropy versus temperature was obtained through integration. This latter could then be compared with the temperature evolution of the width $`\mathrm{\Delta }\nu `$ of the NMR frequency distribution measured by Walisch *et al.*Walisch et al. (1989), which is a quantity proportional to the order parameter. Due to the sharpness of the specific heat anomaly we have assumed that $`T_c`$ corresponds to the temperature of maximum specific heat. The temperature origin of the NMR data was slightly modified to make $`T_c`$ coincide. Excess entropy was set to zero at $`ϵ=0`$. Finally, an interpolation function was used to get the values of the NMR frequency distribution width $`\mathrm{\Delta }\nu `$ at the temperatures values were entropy was determined. Figure 5 presents the correlation between both quantities in a $`\mathrm{\Delta }s`$ versus $`\mathrm{\Delta }\nu ^2`$ graph. The deviation from a mean field linear behaviour is rather small, but as shown in figure 1, XY–3D asymptotic behaviour implies in this representation small changes with respect to the Landau mean field behaviour.

## IV Discussion

### IV.1 Comparative analysis of experimental data

Figure 6 depicts the log-log plot of $`\mathrm{\Delta }s/Q^2`$ versus $`ϵ`$ for the three compounds discussed in section III. One can clearly see the deviation of both Rb2ZnCl4 and Rb2CoF4 from the Landau mean field behaviour, in contrast with the case of SrTiO3. In accordance with the simulations in section II, the curves in general can be approximated to straight lines when reduced to the interval of smaller $`ϵ`$ values. Their slopes yield “experimental effective exponents” $`\kappa _e`$, listed in table 2, that can be compared with the theoretical ones listed in table 1. It can be stated that in the case of Rb2ZnCl4 the experimental effective exponent $`\kappa _e`$ value coincides with the theoretical value of the XY–3D class considering its standard deviation. It should be noted that the values of $`\kappa _{e_T}`$ are similar for XY–3D and Ising–3D universality classes. However, the value $`\beta `$ that fits the order parameter behaviour is consistent with only the XY–3D class. The $`\kappa _e`$ value for Rb2CoF4 is close to $`\kappa _{e_T}`$ expected for Ising–2D class, with a deviation of the order of 5%; while in SrTiO3 not only the experimental $`\kappa _e`$ clearly disagrees with the theoretical value of the Heisenberg–3D class, but its zero value indicates a classical Landau behaviour. It should be stressed that a mere deviation from the linearity between $`\mathrm{\Delta }s`$ and $`Q^2`$ should not be taken as a proof of fluctuation-driven critical behaviour, since a generalized Landau potential with temperature dependent quartic and higher order terms could explain such behaviour. It is the quantitative agreement of this deviation (measured by the effective exponent $`\kappa _e`$) with that expected for the corresponding universality class ($`\kappa _{e_T}`$), together with the agreement of the fitted exponent $`\beta `$ that proofs the consistency of an interpretation in terms of asymptotic critical behaviour for both the order parameter and the specific heat.

### IV.2 Crossover in Rb2ZnCl4

The curve of Rb2ZnCl4 in figure 6 strongly deviates from its low $`ϵ`$ linear behaviour for $`ϵ>0.01`$ and changes then rapidly into a near horizontal straight line. This suggests a crossover around this temperature from the XY–3D critical behaviour to a mean-field/Landau regime. However, this interpretation would be somehow in contradiction with the fact that the observed order parameter behaviour, used in fact for the calculation of this curve, does not suffer any kind of crossover at $`ϵ0.01`$ and follows the $`\beta `$ power law with excellent agreement up to much larger $`ϵ`$ values, of the order of 0.1. In any case, the curve clearly deviates for $`ϵ>0.1`$ from the theoretical curve for the XY–3D class in figure 2 and this can then only be caused by the deviation beyond this $`ϵ`$ value of the experimental $`\mathrm{\Delta }c`$ from the theoretical asymptotic behaviour given by equation (5) and used to produce figure 2. The reason for such deviation could be intrinsic, i.e. the critical region for the specific heat would be then of the order 0.01, while for the order parameter is about 0.1 in reduced temperature, but it is also possible that the critical region were sample dependent. In order to get a more comprehensive analysis of this phase transition we have also analyzed the measurement done by Haga *et al.*Haga et al. (1995) where the excess specific heat was fitted over $`30\mathrm{K}`$, i.e. 0.1 in reduced temperature, using a complex expression involving corrections up to second order of the pure asymptotic regime:

$$\mathrm{\Delta }c=Aϵ^\alpha \left(1+D_1ϵ^{\mathrm{\Delta }_1}+D_2ϵ\right)+B_c$$
(7)
where $`\alpha =0.007`$, $`A=8.153,D_1=0.040,D_2=0.026`$, $`B_c=8.201,\mathrm{\Delta }_1=0.524`$; where $`A,D_1,D_2,B_c`$ are expresed in arbitrary units. This curve is compared in figure 7 with that obtained from a mere critical law as given in equation (5), which essentially corresponds to the neglect of the coefficients $`D_1`$ and $`D_2`$ in equation (7).

The difference between both curves become obvious beyond $`ϵ=0.01`$, showing that for lower temperatures corrections to the asymptotic behaviour are significant. For comparison, we have computed again the curve of $`\mathrm{\Delta }s/Q^2`$ using as specific heat data the values given by the function equation (7) as proposed by Haga. The results are presented in figure 8.

The agreement between both calculations is excellent. The approximate linear behaviour for $`ϵ`$ values smaller than 0.01 is reproduced and its slope gives a value of around 0.15 similar to that obtained from our specific heat analysis. On the other hand, a small shift along the $`ϵ`$-axis of the crossover region can be observed and the behaviour for larger $`ϵ`$ values slightly differs. These differences are most probably due to the choice of the baseline for obtaining the “excess” specific heat. The influence of the baseline is slight at small values of $`ϵ`$ but becomes very important when $`ϵ`$ increases. The results presented for Rb2ZnCl4 in figure 6 correspond to a baseline described by a polynomial of minimal order fitting the specific data well above and below the phase transition. This polynomial can be significantly changed, by shifting the low-temperature fitting region to higher o lower temperatures, while keeping a physically consistent baseline. The effect of this change of baseline is shown in figure 9. It is slight for low $`ϵ`$ values, where the approximate linear behaviour remains unchanged, while the crossover to what one might consider a Landau regime (the horizontal part of the curve in figure 6) becomes less pronounced. As shown in the figure, the change is such that the calculated curve coincides with the theoretical one of the XY-3D class in a much larger $`ϵ`$ interval. Hence, a large uncertainty exists with respect to the actual $`ϵ`$-interval where a crossover to the Landau behaviour does take place. But, the behaviour of the $`\mathrm{\Delta }s/Q^2(ϵ)`$ relationship for $`ϵ<0.01`$ confirming a XY-3D critical regime is robust with respect to modifications of the baseline.

###### Acknowledgements.

This work was supported by Spanish *Ministerio de Ciencia y Tecnología* project number PB98-115.