text:
# The plasma-insulator transition of spin-polarized Hydrogen

## Abstract

A mixed classical-quantum density functional theory is used to calculate pair correlations and the free energy of a spin-polarized Hydrogen plasma. A transition to an atomic insulator phase is estimated to occur around $`r_s=2.5`$ at $`T=10^4K`$, and a pressure $`P0.5Mbar`$. Spin polarization is imposed to prevent the formation of $`H_2`$ molecules.

Département de Physique des Matériaux (UMR 5586 du CNRS),
Université Claude Bernard-Lyon1, 69622 Villeurbanne Cedex, France
Department of Chemistry, University of Cambridge
Lensfield Road, Cambridge CB2 1EW, UK

PACS numbers: 31.15Ew, 61.20.-p, 64.70.-p

Although Hydrogen is generally considered to be the simplest of elements, its expected metallization under pressure  has proved a rather elusive transition. It is now accepted that the behaviour of solid and fluid molecular Hydrogen ($`H_2`$) may be very different at high pressures. Despite considerable experimental efforts with static, room temperature compression of solid $`H_2`$ in diamond anvils  beyond pressures of 2 Mbar, there is still no compelling evidence for a metallic state . The situation is somewhat more favourable for fluid $`H_2`$, since shock compression to 1.4Mbar, and a temperature of about 3000K, led to measurements of metallic resistivities . However, theoretical interpretation is hampered by the absence of a clear-cut scenario; in particular it is not clear whether molecular dissociation precedes ionization or conversely . The presence of several species, $`H_2`$, $`H_2^+`$, H, $`H^+`$ and electrons at a “plasma phase transition”  complicates a theoretical analysis considerably; to gain a clearer picture of pressure-induced ionization, it may be instructive to consider a model system, which would not involve molecular dissociation.

The model system considered in this letter is spin-polarized Hydrogen. If all electron spins are assumed to be polarized by a strong external magnetic field $`𝐁`$, only triplet pair states $`{}_{}{}^{3}\mathrm{\Sigma }`$ can be formed, preventing the binding into $`H_2`$ molecules. The low pressure phase will be made up of H atoms and the only possible scenario upon compression will be the ionization of atoms to form an electron - proton plasma, which is expected to be crystalline at low temperature and fluid at higher temperatures. A rough estimate of the magnetic field needed to spin - polarize the electrons is obtained by equating the magnetic coupling energy $`\mu _BB`$ (where $`\mu _B`$ is the magnetic moment of an electron) to the difference between the triplet and singlet H-H potential energy functions, calculated at the equilibrium distance of the $`H_2`$ molecule ; this leads to $`B10^5Tesla`$. This value exceeds the highest magnetic fields achievable in a laboratory by 3 orders of magnitude, but is well within the range of astrophysical situations. The present calculation neglects possible orbital effects due to a strong applied magnetic field; such effects are expected to be small and their inclusion would lead to a much more involved and less transparent calculation. We prefer to think of our model system as a plasma which has been prepared in a spin-polarized state, and is assumed to remain such even when B is switched off. The subsequent calculation will be restricted to fluid Hydrogen.

The thermodynamic properties of the low pressure atomic $`H`$ phase may be easily calculated from the known triplet pair potential  by standard methods of the theory of classical fluids. We calculated the atom-atom pair distribution function $`g(r)`$ from the HNC integral equation , and deduced from it the equation of state via the virial and compressibility routes. The resulting excess free energies per atom are plotted in Fig.4 as a function of the usual density parameter $`r_s=a/a_0`$, along the isotherm $`T=10^4K`$; here $`a_0`$ is the Bohr radius, and $`a=[3/(4\pi n)]^{1/3}`$, where $`n`$ is the number of $`H`$ atoms per unit volume. There is a thermodynamic inconsistency, typical of HNC theory, but the small difference between the “virial” and “compressibility” free energies will have no influence on our conclusions. To allow for a meaningful comparison with the free energy calculated for the high pressure plasma phase, the free energies shown in Fig.4 contain an electron binding energy contribution of $`0.5a.u.`$ . It is implicitly assumed that this binding energy, valid for isolated atoms (i.e. in the limit $`r_s\mathrm{}`$) does not change upon compression up to $`r_s=2.5`$, due to overlap and distortion of the individual electron $`1s`$ orbitals.

A statistical description of the high pressure phase is more challenging. The key parameter characterizing the spin-polarized electron component is its Fermi energy $`ϵ_F=2.923/r_s^2a.u.`$ ;the corresponding Fermi temperature $`T_F\mathrm{9.2\hspace{0.17em}10}^5/r_s^2K`$. Along the isotherm $`T=10^4K`$ considered in the present calculations, the electrons may be considered to be completely degenerate (i.e. in their ground state) up to $`r_s3`$. The degeneracy temperature of the protons is $`2000`$ times smaller, so that for $`T=10^4K`$, the latter may be considered as being essentially classical, down to $`r_s0.5`$. The proton component is characterized by the Coulomb coupling constant $`\mathrm{\Gamma }=e^2/(ak_BT)=31.56/r_s`$ along the above isotherm, showing that classical Coulomb correlations are expected to be strong over the density range $`1r_s3`$ considered in this paper. Note that while $`\mathrm{\Gamma }`$ decreases as $`r_s`$ increases, the corresponding electron Coulomb coupling constant $`\gamma =e^2/(aϵ_F)=0.342r_s`$ increases.

In the ultra-high density regime $`r_s1`$, the electron kinetic energy dominates, and the proton and electron components decouple in first approximation (“two-fluid” model); the weak proton-electron coupling may be treated by linear response theory , suitably adapted to the spin-polarized case. Within linear response, the free energy per atom (ion-electron pair) splits into three terms: the ground-state energy of the uniform, spin-polarized electron gas (“jellium”), $`ϵ_e`$, the free energy of protons in a uniform neutralizing background (the so-called “one-component plasma” or OCP), $`f_{OCP}`$, and the first order correction due to linear screening of the Coulomb interactions by the electron gas, $`\mathrm{\Delta }f`$:

$$f=\frac{F(\mathrm{\Gamma },r_s)}{N}=ϵ_e+f_{OCP}+\mathrm{\Delta }f$$
(1)
where $`ϵ_e(r_s)`$ is taken to be the sum of kinetic ($`1.754/r_{s}^{}{}_{}{}^{2}`$), exchange ($`0.5772/r_s`$) and correlation  contributions; $`f_{OCP}`$ is given by an accurate fit to Monte Carlo simulations of the OCP ; $`\mathrm{\Delta }f`$ follows from first order thermodynamic perturbation theory :

$$\mathrm{\Delta }f=\frac{1}{2(2\pi )^3}S_{OCP}(k)\widehat{w}(k)𝑑𝐤$$
(2)
where $`S_{OCP}(k)`$ is the static structure factor of the OCP (which plays the role of reference system). According to linear response theory, $`\widehat{w}(k)`$ is the difference between screened and bare ion-ion pair potentials:

$$\widehat{w}(k)=\frac{4\pi e^2}{k^2}\left[\frac{1}{ϵ(k)}1\right]$$
(3)
where $`ϵ(k)`$ is the dielectric function of the electron gas which we calculated within the RPA from the Lindhard susceptibility of a gas of spin-polarized, non-interacting electrons, supplemented by a local exchange and correlation correction . All necessary ingredients for the calculation of $`\mathrm{\Delta }f`$ may be found in , and the resulting free energy curve is shown in Fig.4. Although linear response cannot, a priori, be expected to be quantitatively accurate for $`r_s>1`$, it provides a rough estimate of the plasma to atomic phase transition, from the intersection of the free energy curves, which is seen to occur at $`r_s1.9`$. The corresponding transition pressure would be $`2.3Mbar`$.

However, as $`r_s`$ increases, the ion-electron coupling becomes stronger, and the non-linear response of the electron component to the “external” potential field provided by the protons, is expected to lower the free energy of the plasma phase. To explore the non-linear regime we have adapted the HNC-DFT formulation of our earlier work on (unpolarized) metallic H  to the spin-polarized case. Within this formulation, proton-proton and proton-electron correlations are treated at the HNC level, which is expected to be a good approximation for the long-range Coulomb interactions, while the energy of the inhomogeneous electron gas follows from the density functional ($`E=Nϵ_e`$):

$$E\left[\rho (𝐫)\right]=E_K\left[\rho (𝐫)\right]+E_H\left[\rho (𝐫)\right]+E_X\left[\rho (𝐫)\right]+E_C\left[\rho (𝐫)\right]$$
(4)
where $`\rho (𝐫)`$ denotes the local electron density, and $`E_K`$, $`E_H`$, $`E_X`$ and $`E_C`$ are the kinetic, Hartree, exchange and correlation contributions. For $`E_K`$ we adopted the Thomas-Fermi approximation, corrected by a square gradient term:

$$E_K\left[\rho (𝐫)\right]=C_K[\rho (𝐫)]^{5/3}𝑑𝐫+\frac{\lambda }{8}\frac{|\rho (𝐫)|^2}{\rho (𝐫)}𝑑𝐫$$
(5)
where $`C_K=3(6\pi ^2)^{2/3}/10a.u.`$, while the choice of $`1/9<\lambda <1`$ will be specified below. The mean field Hartree term is of the usual form:

$$E_H\left[\rho (𝐫)\right]=\frac{1}{2}𝑑𝐫𝑑𝐫^{}\frac{\mathrm{\Delta }\rho (𝐫)\mathrm{\Delta }\rho (𝐫^{})}{|𝐫𝐫^{}|}$$
(6)
where $`\mathrm{\Delta }\rho (𝐫)=\rho (𝐫)n`$, while:

$$E_X[\rho (𝐫)]=C_X[\rho (𝐫)]^{4/3}𝑑𝐫$$
(7)
with $`C_X=3(6/\pi )^{1/3}/4a.u.`$ . The correlation contribution $`E_C[\rho (𝐫)]`$ (within the LDA) can be found in . This functional yields an explicit form for the electron-electron direct correlation function $`c_{22}(r)`$ (henceforth the indices 1 and 2 will refer to protons and electrons respectively) . The remaining direct and total correlation functions $`c_{11}(r)`$, $`c_{12}(r)`$, $`h_{11}(r)`$ and $`h_{12}(r)`$ are calculated by a numerical resolution of the HNC closure equations and the quantum version of the Ornstein-Zernike (OZ) relations , which form a closed set of coupled non-linear integral equations for the four functions.

Solutions were obtained by a standard iterative procedure along the isotherm $`T=10^4K`$ and for density parameters in the range $`0.5r_s2.5`$, corresponding to more than one-hundred-fold compression of the lowest density state ($`r_s=2.5`$), which would correspond to $`0.17gr/cm^3`$. The temperature is roughly equal to that expected inside Saturn, and comparable to temperatures reached in shock compression experiments on the NOVA laser facility at Livermore . The iterative solutions were first obtained at the highest densities ($`r_s=0.5`$ and $`1`$), where linear response theory provides reasonably accurate intial input. The prefactor $`\lambda `$ in the square gradient correction to the electron kinetic energy functional (5) was adjusted to provide the best match between the HNC-DFT result for the local radial density of electrons around a proton, $`r^2g_{12}(r)=r^2[1+h_{12}(r)]`$, and its linear response prediction, at the highest densities ($`r_s=0.5`$ and $`1`$), where linear response should be most accurate; $`g_{12}(r)`$ turns out to be rather sensitive to $`\lambda `$, and the best agreement is achieved for $`\lambda =0.18`$, which is close to the value $`1/5`$ frequently advocated in electronic structure calculations for atoms .

Results for the local radial density $`r^2g_{12}(r)`$ are shown in Fig.1 for several values of $`r_s`$. As expected, electrons pile up increasingly at small $`r`$ as $`r_s`$ increases, and a shoulder is seen to develop around $`r/a0.4`$. For comparison, the linear response prediction is shown at $`r_s=1`$, while at the lowest density ($`r_s=2.5`$), an estimate of $`g_{12}(r)`$ in the atomic phase is obtained by adding to the electron density in a $`H`$ atom (namely $`r^2\rho (r)=r^2\mathrm{exp}(2r)/\pi a.u.`$) the convolution of the latter with the atom-atom pair distribution function $`g(r)`$. HNC-DFT results for the proton-proton pair distribution function are shown in Fig.2 for three densities. As expected, proton-proton correlations are seen at first to weaken, as the density decreases, due to enhanced electron screening. However at the lowest density ($`r_s=2.5`$), weakly damped oscillations build up at long distances, which may be indicative of an incipient instability of the proton-electron plasma. The atom-atom $`g(r)`$ at the same density agrees reasonably well with $`g_{11}(r)`$ up to the first peak, but it does not exhibit the long-range correlations in the latter. The proton-proton structure factors $`S_{11}(k)`$ are plotted in Fig.3 for several values of $`r_s`$. A considerable qualitative change is again seen to occur at the lowest density ($`r_s=2.5`$), where the main peak is shifted to larger $`k`$, while a significant peak builds up at $`k=0`$. Such enhanced “small angle scattering” is reminiscent of the behaviour observed in simple fluids near a spinodal (subcritical) instability. In fact we were unable to obtain convergence of the HNC-DFT integral equations for $`r_s>2.5`$, which hints at an instability of the electron-proton plasma at lower densities. This strongly suggests a transition to the insulating atomic phase, but the simple density functional used in this work cannot properly describe the recombination of protons and electrons into bound (atomic) states .

In order to confirm this scenario, the free energy of the plasma phase should be compared to that of the atomic phase. This is easily achieved within the high density linear response regime, as shown earlier. However the calculation of the free energy in the non-linear regime appropriate for lower densities ($`r_s>1`$) is less straightforward . In fact the present HNC-DFT formulation provides only one direct link with thermodynamics, namely via the compressibility relation :

$$\underset{k0}{lim}S_{11}(k)=\underset{k0}{lim}S_{12}(k)=nk_BT\chi __T$$
(8)
where $`\chi __T`$ denotes the isothermal compressibility of the plasma. From the calculated values of $`\chi __T`$, the free energy of the plasma follows by thermodynamic integration, starting from a reference state (e.g. $`r_s=1`$) for which the linear response estimate is expected to be accurate. The resulting “compressibility” free energy curve is plotted in Fig.4. Somewhat unexpectedly it lies above the linear response prediction. An alternative route to the free energy is via the virial relation for the pressure; only an approximate virial expression is known within the present HNC-DFT formulation , and the resulting “virial” free energy curve is also shown in Fig.4. It falls well below the “compressibility” free energy curve, thus illustrating the well known thermodynamic inconsistency of the HNC closure for Coulombic fluids . Any reasonable extrapolation of the two free energy curves would miss the low density limit $`0.5a.u.`$ by as much as twenty per cent. We suggest instead an estimate of the free energy of the plasma phase by taking the average of the “compressibility” and “virial” values, despite the lack of foundamental a priori justification for doing this. A short extrapolation of the resulting curve is likely to intersect or smoothly join on to the free energy of the atomic phase just beyond $`r_s=2.5`$. An intersection would correspond to a first-order phase transition, reminiscent of the “plasma phase transition” of Saumon and Chabrier . However, due to the uncertainty on the thermodynamics of the plasma phase, a continuous transition cannot be ruled out. The transition pressure $`P`$ would be of the order of $`0.5Mbar`$, well below the current experimental and theoretical estimates for the transition of fluid molecular Hydrogen to a conducting state .

In summary the structure and thermodynamic results derived from an HNC-DFT theory of the spin-polarized proton-electron plasma strongly suggest that this plasma will recombine into an insulating atomic phase at $`r_s2.5`$, for a temperature $`T=10^4K`$. We are presently exploring the behaviour of the system at lower temperatures.


format: text

text:
# Vibrational properties of the one-component 𝜎 phase

## I Introduction

The local atomic order in disordered condensed materials is well defined and governs many physical properties . Quite often, for a disordered marerial, it is possible to find a corresponding crystal with similar local and even intermediate-range order which give rise to similarities in many structural and dynamical features of these two solids. Such a crystal can be regarded as a reference crystalline structure (crystalline counterpart) for the corresponding disordered substance. In some cases, the reference structure can be uniquely defined. The simplest examples are toy structural models with force-constant and/or mass disorder. In these toy models, the atoms occupy their equilibrium positions at the sites of a crystalline lattice (e.g. simple cubic), which can be considered to be a reference one (see e.g. ). Another related example is a binary substitutional alloy, the reference system for which is a periodic point lattice with one of the two atomic species placed at the lattice sites . The disorder in such models does not influence the equilibrium positions of the atoms arranged in an ideal crystalline lattice. This makes possible the use of approximate analytical approaches (e.g. the coherent potential approximation ) to treat the vibrational properties of the models, provided that the vibrational properties of the counterpart crystal are known.

In amorphous solids, or glasses, the atoms do not occupy the sites of a crystalline lattice, which results in positional disorder. For these materials, a choice of a reference structure becomes problematic. Good counterparts can usually be found among the crystalline polymorphs having the same (or similar) chemical composition as the corresponding glass. For example, $`\alpha `$-cristobalite appears to be a good crystalline counterpart for vitreous silica .

The main purpose of this paper is to investigate numerically the vibrational properties of a one-component $`\sigma `$-phase  crystal which is conjectured to be a good crystalline counterpart for a one-component glass with icosahedral local order (IC glass) .

The motivation for this choice of a crystalline counterpart of the IC glass is the following. The computational model of the IC glass is based on a simple empirical pair interatomic potential  resembling the effective interionic potentials conjectured for simple metallic glass-forming alloys . The use of the same potential allows us to construct models of bcc and $`\sigma `$-phase crystals that are stable for a wide range of thermodynamical parameters . Of these two crystalline structures, the $`\sigma `$ phase is expected to be a good reference structure for the IC glass because of the following reasons: The supercooled IC liquid (where the interactions between atoms are described by the same potential ) undergoes a transition either to the IC glass or to a dodecagonal quasicrystal  depending on the quench rate . This quasicrystal has similar local structural properties with the IC glass . However, the absence of global periodic order in the quasicrystalline phase makes the analysis of its vibrational properties a task of comparable complexity to that for the glass itself. The $`\sigma `$ phase is one of the closest low-order crystalline approximants  for this dodecagonal quasicrystal , which means that these two (crystalline and quasicrystalline) structures are built up from the same structural units. This implies that the IC glass and the $`\sigma `$ phase, being both tetrahedrally closed-packed structures , are nearly isomorphous in terms of local order.

Knowledge of the vibrational properties of the $`\sigma `$-phase crystal allows for a direct comparison with those of the IC glass. The apparent similarity in the vibrational densities of states of these two structures gives stronger support for the choice of this crystalline counterpart for the IC glass.

The $`\sigma `$-phase structure, used in our computations, has been obtained by means of molecular dynamics simulation with the use of an interatomic pair potential . The vibrational properties have been investigated by using both a the normal-mode analysis and by computing the spectra of appropriate time-correlation functions.

The paper is arranged as follows: The $`\sigma `$-phase structure is described in Sec. II. The model and technical details of the calculations are presented in Sec. III. In Sec. IV we present the results of the simulations. Some concluding remarks are contained in Sec. V.

## II The $`\sigma `$ phase

In this Section, we review the known structural and dynamical properties of the $`\sigma `$ phase.

### A Structure

The $`\sigma `$ phase belongs to an important class of tetrahedrally close-packed crystallographic structures , viz. the Frank-Kasper phases . The first coordination shells of the constituent atoms in these structures form triangulated (Frank-Kasper) polyhedra composed entirely of slightly distorted tetrahedra. The four possible coordination numbers ($`Z`$) in these structures are $`Z=12,14,15`$ and $`16`$. The least distorted tetrahedra are found in icosahedra ($`Z12`$ polyhedra). Structures of small clusters of atoms interacting via pairwise central potentials favor icosahedral order  as having the lowest energy. The prototype $`\sigma `$ phase structures are $`\beta `$-U  and Cr48Fe52 . There are $`30`$ atoms per tetragonal unit cell ($`tP30`$) with $`c/a0.52`$, where $`c`$ and $`a`$ are the dimensions of the cell (lattice parameters). The space group of this phase is $`P4_2/mnm`$. There are $`10`$ atoms with the coordination number $`12`$, $`Z12`$ or icosahedra, $`16`$ $`Z14`$ atoms and $`4`$ $`Z15`$ atoms. The $`72^{}`$ disclination lines  form a network (a major skeleton, in the parlance of Frank and Kasper ), where rows of $`Z14`$ atoms parallel to the tetragonal $`c`$-axis thread planar networks of $`Z14`$ and $`Z15`$ atoms. A projection of the $`\sigma `$-phase structure down the $`c`$-axis is shown in Fig. 1.

The Frank-Kasper phases share their significant geometrical property of icosahedral local order with simple metallic glasses . Some liquid alloys which form Frank-Kasper phases have a tendency to freeze into metastable amorphous structures (metallic glasses) when quenched sufficiently rapidly . It is now generally well accepted that, at least in the case of metallic alloys of simple constitution, glass formation is caused by the incompatibility of local icosahedral coordination with the translational symmetry in Euclidean space (geometrical frustration) . There exist statistical mechanical arguments in favor of this scenario of glass formation  based on a Landau free-energy analysis. The average coordination number, $`\overline{Z}`$, in the Frank-Kasper phases ($`13.333\overline{Z}13.5`$) is very close to that of a sphere-packed “ideal glass” model  ($`\overline{Z}_{\mathrm{ideal}}=13.4`$). In a sense, this “ideal glass” could be regarded as a Frank-Kasper phase with an infinitely large unit cell . Thus the class of Frank-Kasper phases is a natural choice for reference crystalline structures for metallic glasses of simple constitution.

From a structural point of view, the $`\sigma `$ phase can be also regarded as a crystalline low-order approximant for dodecagonal quasicrystals . Such quasicrystals , morphologically close to Frank-Kasper phases, represent an alternative class of noncrystallographic structures which combine icosahedral local order with non-translational long-range order manifested by infinitely sharp diffraction peaks.

### B Dynamics

Similarities in the local structure of metallic glasses and Frank-Kasper phases are reflected in the dynamical properties of these materials.

The available data about the vibrational dynamics of Frank-Kasper phases is limited to some of the Laves phases, a subclass of the Frank-Kasper phases . For instance, the similarity between the phonon-dispersion relations of the Mg70Zn30 glass  and those of the Laves phase MgZn2  was emphasized by Hafner .

Another interesting aspect of the dynamics of Frank-Kasper phases is related to the appearance of soft vibrational modes in these materials. Such a soft low-frequency optic mode at the $`\mathrm{\Gamma }`$-point (the origin of the reciprocal lattice) has been found numerically in the same Laves phase MgZn2 . The frequency of this mode decreases with increasing pressure (accompanied by volume compression) and eventually becomes negative, indicating a structural phase transition . The authors of Ref.  applied a group-theoretical analysis and demonstrated that the polarization vector of the soft optic mode in MgZn2 is determined by the structure symmetry and is independent of interatomic interactions. This suggests that the soft-mode character of some vibrations is a generic property of Frank-Kasper phases. However, no soft-mode behavior was observed in the isostructural CaMg2  where the mass ratio of the constituent elements is reversed with respect to MgZn2. Hafner  suggested that the soft modes in MgZn2 should be attributed to the relatively large mass of the Zn atoms. This is an example of the chemical composition of the materials introducing considerable difficulties into the analysis of the interplay of the structure and the dynamics.

A numerical simulation of a one-component Frank-Kasper phase allows us to eliminate this uncertainty. We investigate the behavior of the lowest-frequency optic modes in the one-component $`\sigma `$ phase with variable pressure in Sec. IV B 3.

## III Methods of computation

We have constructed a thermodynamically stable structural model of the one-component $`\sigma `$ phase by means of classical molecular dynamics simulation. In this case, success strongly depends on the choice of the interatomic potential for the model. For example, the Lennard-Jones potential, widely used in creating simple models of liquids , glasses , and crystals , is not suitable for this purpose, because the $`\sigma `$ phase is not stable for it in the range of thermodynamical parameters investigated below (the stable phase is the fcc lattice). Instead, we use a pair interatomic potential suggested in Ref.  and show that it is possible to construct a one-component $`\sigma `$ phase which is stable for a wide range of thermodynamical parameters .

### A Model

As a mathematical model for the study of atomic dynamics in the crystalline $`\sigma `$-phase structure, we consider a classical system of $`N`$ identical particles interacting via a spherically symmetric pair potential.

The pair potential used in this study  is designed to favor icosahedral local order. The main repulsive part of this potential is identical to that of the Lennard-Jones potential $`u_{\mathrm{LJ}}(r)=4ϵ[(\sigma /r)^{12}(\sigma /r)^6]`$; therefore all the quantities in these simulations are expressed in reduced Lennard-Jones units , i.e. with $`\sigma `$, $`ϵ`$, and $`\tau _0=(m\sigma ^2/ϵ)^{1/2}`$ chosen as length, energy, and time units. To convert the reduced units to physical units one can refer to argon ($`m`$ = 39.948 a.m.u.) by choosing the Lennard-Jones parameters $`\sigma =3.4`$ Å and $`ϵ/k_B=120`$ K. In this case, our frequency unit $`\nu _0=\tau _0^1`$ corresponds to 0.4648 THz. The analytical expression defining the potential is given in Ref. . This potential resembles those for simple glass-forming metallic alloys  with only the first of the Friedel oscillations being retained (see Fig. 2).

We use conventional molecular dynamics simulations  in which the Newtonian equations of motion are solved using a finite-difference algorithm with time step equal to 0.01 while the particles are enclosed in a simulation box of volume $`V`$ with periodic boundary conditions. In this case, the total energy $`E`$ is a constant of motion and time averages obtained in the course of simulations approximate the ensemble averages in the microcanonical (constant-$`NVE`$) statistical ensemble. Wavevectors,

$$𝐐=n_x𝐐_{x,0}+n_y𝐐_{y,0}+n_z𝐐_{z,0}$$
(1)
with $`n_x`$, $`n_y`$ and $`n_z`$ being integers, compatible with periodic boundary conditions are multiples of the three fundamental wavevectors $`𝐐_{x,0}=\frac{2\pi }{L_x}(1,0,0)`$, $`𝐐_{y,0}=\frac{2\pi }{L_y}(0,1,0)`$ and $`𝐐_{z,0}=\frac{2\pi }{L_z}(0,0,1)`$, where $`L_x`$, $`L_y`$, and $`L_z`$ are the dimensions of the (tetragonal) simulation box. In order to have a sufficient number of allowed wavevectors within the Brillouin zone, the sample dimensions must be sufficiently large. The time-correlation functions resulted from the molecular dynamics simulation reported below were obtained for a system of 20580 particles ($`7\times 7\times 14`$ unit cells with $`30`$ atoms per unit cell). This relatively large system size also gives sufficient statistical accuracy . Where necessary, we performed simulations in other ensembles by modifying the equations of motion .

The arrangement of atoms in a unit cell of the model $`\sigma `$-phase structure used in our computer simulations is shown in fig. 1(b). The optimal (with respect to an energy minimization) $`c/a`$ ratio was taken to be equal to 0.5273. To determine the optimal model structure, it was sufficient to perform molecular dynamics simulations with the number of particles in the system equal to $`N=1620`$ ($`3\times 3\times 6`$ unit cells). Details of the preparation of this atomic configuration are given in section IV A.

### B Time-correlation functions

A straightforward method to analyze the vibrational dynamics in a molecular dynamics model is to imitate inelastic neutron scattering experiments by calculating the dynamical structure factor $`S(𝐐,\omega )`$ , proportional to the neutron scattering cross-section , which is the spectrum of the density autocorrelation function:

$$F(𝐐,t)=<\rho (𝐐,t)\rho (𝐐,0)>$$
(2)
where

$$\rho (𝐐,t)=\underset{k=1}{\overset{N}{}}\mathrm{exp}(i𝐐𝐫_k(t))$$
(3)
is the Fourier transform of the local particle density , $`N`$ is the number of particles in the system, $`𝐫_k(t)`$ is the position vector of particle $`k`$, and the wavevector $`𝐐`$ takes on the values according to Eq. (1). A longitudinal phonon is associated with a maximum in $`S(𝐐,\omega )`$ at a fixed $`𝐐`$. In order to get information about the transverse modes from $`S(𝐐,\omega )`$, one has to select wavevectors outside the first Brillouin zone .

In a more convenient way, the vibrational modes can be studied using the current autocorrelation function :

$$C_𝐞(𝐐,t)=\frac{Q^2}{N}<j_𝐞(𝐐,t)j_𝐞(𝐐,0)>$$
(4)
where

$$j_𝐞(𝐐,t)=\underset{k=1}{\overset{N}{}}(𝐞𝐯_k(t))\mathrm{exp}[i𝐐𝐫_k(t)]$$
(5)
is the Fourier transform of the local current, $`𝐯_k(t)`$ is the velocity of particle $`k`$, and $`𝐞`$ is the unit polarization vector.

Note that for the longitudinal polarization, $`𝐞𝐐`$, Eq. (4) can be obtained from (2) by double differentiation with respect to time. For the transverse-current correlation functions, the polarization vector must be chosen consistently with the lattice symmetry.

At a temperature $`T`$, the vibrational density of states $`g(\omega )`$ can be calculated as the Fourier transform of the normalized velocity autocorrelation function :

$$Z(t)=\frac{1}{3NT}\underset{k=1}{\overset{N}{}}𝐯_k(t)𝐯_k(0)$$
(6)
We computed the time-correlation functions using the overlapped data collection technique . The number of overlapped measurements used for statistical averaging was about 10000. The time origins of the measurements were separated by 0.2 r.u. (20 time steps).

In order to reduce the finite-time truncation effects in the spectra of the time-correlation functions, we used a Gaussian window function with the half-width equal to 3 r.u.

### C Normal-mode analysis

To calculate the dispersion relations in the harmonic approximation, we used the standard method based on diagonalization of the Fourier transformed dynamical matrix . From the known dispersion relations, $`\omega _s(𝐐)`$, $`s=1,2,\mathrm{},90`$, the vibrational density of states can be computed by integration over the first Brillouin zone according to

$$g(\omega )=\frac{v}{r(2\pi )^3}\underset{s=1}{\overset{r}{}}_{\mathrm{BZ}}\delta (\omega \omega _s(𝐐))𝑑𝐐$$
(7)
where the sum is over all $`r`$ dispersion branches and $`v`$ stands for the volume of the unit cell. In the computations, a Gaussian function with small but finite width is substituted for the $`\delta `$-function. The half-width of the Gaussian function in our computations was equal to about 0.05. Formally, the wavevector $`𝐐`$ in Eq. (7) is a continuous variable, but in the simulations the integral was estimated by a sum over a uniform rectangular grid of $`100\times 100\times 100`$ points in the first Brillouin zone.

## IV Results

### A Optimization of the structure in the $`\sigma `$ phase

In this subsection, we describe the method of construction of the thermodynamically stable model of the $`\sigma `$ phase and analyse the range of its stability.

There are several ways to obtain numerical values of the atomic coordinates in the $`\sigma `$ phase. One way is to use those available for Cr48Fe52 . Alternatively, a unit cell of the $`\sigma `$-phase structure can be constructed either by manipulating the kagomé tiling according to the algorithm given by Frank and Kasper in ref.  or by stacking the square and triangular basic elements of the dodecagonal quasicrystal model  into the $`3^2,4,3,4`$ square-triangle net  (see fig. 1(a)). The arrangements of the atoms resulting from these constructions do not correspond exactly, although the difference is rather small – the root-mean-square distance between the corresponding atoms in different configurations is of the order of a few percent of the root-mean-square distance between different atoms in the same configuration. In either case, the resulting structure is an approximate one in the sense that the atomic positions do not correspond to a minimal potential-energy configuration for a given interaction potential. To obtain the true structure corresponding to the potential, the approximate configuration must be relaxed by a molecular dynamics program. Moreover the $`c/a`$ ratio is slightly different for different natural $`\sigma `$-phase crystals, which means that this ratio is not uniquely defined.

The atomic configuration of the $`\sigma `$ phase used in this study was prepared as follows. A sample of $`N=1620`$ particles ($`3\times 3\times 6`$ unit cells) was constructed by filling a tetragonal box of appropriate dimensions with $`\sigma `$-phase unit cells. We used the unit-cell atomic configuration suggested by Gähler  with $`c/a=0.5176`$. The number density, $`\rho =N/V`$, $`V=L_xL_yL_z`$, of this atomic configuration was $`\rho =1.0048`$. This configuration was then used to provide the initial atomic coordinates for variable-shape $`NST`$ (constant number of atoms, pressure tensor, and temperature) molecular dynamics , the run performed at $`S=0`$ and $`T=0`$. This procedure is equivalent to a potential-energy minimization by the steepest descent method under the condition of independent pressure balance in each spatial dimension. The variable-shape $`NST`$ run resulted in an ideal crystalline structure for which the fractional coordinates of the atoms in all unit cells were identical within the precision of the calculations. The structure of the $`\sigma `$ phase thus obtained is characterized by the minimum potential energy per atom, $`U_{\mathrm{min}}=2.5899`$, with respect to variations of thermodynamical parameters. In order to check this, we have performed similar $`NST`$ runs at different pressures and indeed found that the energy is minimum at zero pressure (see Fig. 3a,b). The density for the optimal structure has been found to be $`\rho =0.8771`$ and the ratio of the lattice parameters is $`c/a=0.5273`$, close to that of $`\beta `$-U ($`c/a=0.5257`$, at $`T=720^{}`$C) . The potential-energy minimum for the bcc structure at the same density with the same potential was $`U_{\mathrm{min}}=2.6148`$. At zero pressure, the density of the bcc structure is $`\rho =0.8604`$ and the minimum potential energy per atom is $`U_{\mathrm{min}}=2.6357`$, i.e. in both cases lower than for the $`\sigma `$ phase. It was shown earlier  that the potential energy per atom for the $`\sigma `$ phase becomes lower than that for the bcc structure at the same density as the temperature increases. This is consistent with the fact that natural crystalline $`\sigma `$ phases are stable only at high temperatures. They undergo a solid-solid phase transition to a simpler crystalline phase as the temperature decreases. In our simulations, however, the system was stable in the range of temperatures $`T0.9`$ for as long as $`t_{\mathrm{run}}=5000`$.

We have also investigated the thermodynamical stability of the $`\sigma `$ phase under variable pressure. We have found that the $`\sigma `$ phase is stable for pressures in the range $`5P12`$. At high pressures $`P12`$, a structural transformation occurs, resulting in the fcc structure. The phase diagram of the IC potential is not known at present. We can expect that at densities greater than the triple-point density for the Lennard-Jones system, $`\rho 0.85`$ , the solid-fluid coexistence curve for the IC system is close to that for the LJ system. We can only estimate that the melting temperature at $`\rho =0.8771`$ is about $`0.8T0.9`$ from the fact that the $`\sigma `$-phase crystalline structure is stable at $`T=0.8`$ on the time scale of our computations. No diffusion was observed at temperatures up to $`T=0.8`$. At $`T=0.9`$ the system stayed in a metastable superheated state for about $`t_{\mathrm{run}}=5000`$, after which it melted.

We used the density and the $`c/a`$ ratio obtained from the $`NST`$ run to perform $`NVE`$ (constant number of atoms, volume, and total energy) molecular dynamics runs starting from the three configurations mentioned above and scaling the velocities to zero at each time step, which is also equivalent to a potential-energy minimization by the steepest-descent method. The same was done for one instantaneous configuration corresponding to the temperature $`T=0.8`$. The configurations resulting from this procedure were identical, which is an indication that there is a well-defined potential-energy minimum corresponding to a unique crystallographic arrangement of atoms within the $`\sigma `$-phase unit cell. This structure, scaled so that $`a=c=1`$, is shown in fig. 1(b). The atomic layers with $`z`$ close to 0.25 and 0.75, which are not closely packed, show a small but significant puckering – an effect present in the $`\beta `$-U  and Cr48Fe52  structures.

### B Vibrational dynamics

Above, we have discussed the similarities in the local structure of the $`\sigma `$ phase and the IC glass. These similarities are expected to cause the vibrational spectra in these two materials also to be similar. In order to check this assumption, in this subsection we investigate the vibrational properties of the $`\sigma `$ phase and compare the vibrational spectra of this crystal and the IC glass.

If the vibrational spectra are similar, the $`\sigma `$ phase can be considered to be a good crystalline reference structure for the IC glass. One consequence of these similarities is that we can use the data about the vibrational properticies of the $`\sigma `$ phase crystal to explain the nature of the vibrational excitations in the corresponding amourphous structure.

#### 1 Phonon dispersion in the $`\sigma `$ phase

The $`\sigma `$ phase has $`30`$ atoms per unit cell which result in $`3`$ acoustic and $`87`$ optic branches, as shown in Fig. 4. The vibrational density of states (VDOS) obtained by integration over the first Brillouin zone (see Eq. (9)) is shown on the right-hand side of Fig. 4. The linear dispersion of the acoustic branches in the low-frequency range ($`\omega 4`$) results in the Debye law for the VDOS, $`g(\omega )=3\omega ^2/\omega _\mathrm{D}^3`$, with the Debye frequency equal to $`\omega _\mathrm{D}23.89`$. The Debye frequency has been estimated from a fit of the initial part of the VDOS by a parabolic function. The optic branches are densely distributed above the acoustic part. There are no large gaps in the spectrum, which is a consequence of tight binding and the mutual penetration of the basic structural units (Frank-Kasper polyhedra) in the $`\sigma `$ phase. In other words, there are no isolated structural units, like molecules in molecular crystals  and crystalline fullerens , or tetrahedra in silica , the vibrations of which form separate optic bands. At some of the zone boundaries (e.g. the X point; see Fig. 4) the dispersion curves do not show zero derivatives. This is because the space group $`P4_2/mnm`$ of the $`\sigma `$ phase is nonsymmorphic, i.e. it contains nonpoint symmetry elements involving fractional translations .

#### 2 Comparison with the IC glass

An informative characteristic of the vibrational dynamics in the IC glass which can be compared with the $`\sigma `$ phase is its VDOS  (see Fig. 5) which can be easily obtained from the velocity autocorrelation function Eq. (6). The VDOS for the bcc lattice is also presented for comparison in the figure (the dashed line). We can clearly see that the frequency range of the whole spectrum is the same for the $`\sigma `$ phase and the IC glass but differs for the bcc lattice. The shape of the IC-glass spectrum mainly reproduces the basic features of the $`\sigma `$ phase spectrum and can be imagined as a superposition of broadened (by disorder) crystalline peaks. This is a consequence of the presence of a large number of optic modes in the vibrational spectrum of the sigma-phase structure located in the same frequency region as the whole spectrum of the IC glass. Therefore, the similarities in the VDOS of the $`\sigma `$ phase and the IC glass strongly support the choice of the $`\sigma `$ phase as a crystalline counterpart.

#### 3 Soft modes in the $`\sigma `$ phase

As was mentioned in Sec. II B, an interesting feature of atomic dynamics in the Frank-Kasper phases is related to the appearance of low-frequency soft modes. We have investigated whether a soft mode appears in our model of the $`\sigma `$-phase structure. For this purpose, we followed the evolution of the vibrational spectrum with variable pressure (see Fig. 6 and, indeed, found that one of the lowest-frequency optic modes (doubly degenerate) in the $`\mathrm{\Gamma }`$-point shows soft-mode behavior. The frequency of this mode decreases both with decreasing and increasing pressure (see Fig. 7). The decrease of the mode frequency at negative pressures is not surprising and reflects the softening of the whole vibrational spectrum (see Fig. 6a). However, with increasing pressure, the whole spectrum is shifted to higher frequencies (see Fig. 6c), while the frequency of the soft mode (a small peak around $`\omega =3.5`$ in Fig. 6c) moves in the opposite direction, approaching zero and thus indicating a structural instability (structural phase transition to the fcc lattice) at a critical pressure $`P_{}12.5`$ (see Fig. 7). Around this value of the pressure the structure of the $`\sigma `$ phase becomes extremely unstable and an investigation of the details of atomic motion requires a thorough analysis. We hope to address this point in another study.

#### 4 Anharmonicity in the $`\sigma `$ phase

One of the interesting questions concerning the vibrational dynamics of the $`\sigma `$ phase is related to the range of applicability of the harmonic approximation for the lattice vibrations. We are able to anwer this question by investigating the vibrational spectrum using the velocity autocorrelation function with increasing temperature and comparing it with the results of the normal-mode analysis (harmonic approximation).

To assess the degree of temperature-induced anharmonicity, we computed the dispersion relations for the symmetry direction $`[001]`$ ($`𝐐𝐜`$, $`\mathrm{\Gamma }`$Z in Fig. 4) at different temperatures by using both these techniques (molecular dynamics and normal-mode analysis) and compared the results. These are shown in Fig. 8 for several low- and high-frequency dispersion branches for two temperatures $`T=0.01`$ and $`T=0.8`$. At intermediate frequencies, the density of dispersion branches is so high that a comparison between the results of the two methods of calculation of dispersion relations is hardly possible, mainly because of the finite width of the respective peaks in the spectra of the current autocorrelation functions (see sec. III B). Due to the fact that the $`\sigma `$-phase space group $`P4_2/mnm`$ is nonsymmorphic, i.e. it contains nonpoint symmetry elements involving fractional translations , the phonon-dispersion relations, derived from the peak positions in $`C_l(𝐐,\omega )`$ and $`C_t(𝐐,\omega )`$, appear in the extended zone scheme . The optic modes cannot be measured in the vicinity of the origin of the first Brillouin zone $`Q=0`$ ($`Q=|𝐐|`$), because this long-wavelength limit corresponds to motion of the system as a whole, forbidden by the periodic boundary conditions. Information about these modes is available at the boundaries ($`Q=2\pi /c,6\pi /c`$) and at the origin ($`Q=4\pi /c`$) of the second extended Brillouin zone. The molecular dynamics results for the dispersion relations at $`T=0.01`$ are adapted from the second extended Brillouin zone. To make possible the comparison with the results obtained in the harmonic approximation, the data in the region $`Q[5\pi /c,6\pi /c]`$ were folded with respect to $`Q=5\pi /c`$ into the region $`Q[4\pi /c,5\pi /c]`$, which corresponds to half of an irreducible zone. From these results we can see that the harmonic approximation works quite well at the low temperature of $`T=0.01`$. At $`T=0.8`$, only the acoustic branches could be resolved without ambiguity. Therefore, for this temperature, we used the data available from the first Brillouin zone.

One important signature of temperature-induced anharmonicity is a softening of the acoustic modes, i.e. a lowering of the acoustic branches with respect to those calculated in the harmonic approximation which occurs as the temperature increases . This effect can be clearly seen for $`T=0.8`$ in Fig. 8. In accordance with this observation, the vibrational density of states for this temperature, shown in Fig. 9, exhibits the presence of excess states with respect to the harmonic approximation. A deviation from the harmonic approximation in $`g(\omega )`$ at low frequencies (see the inset in Fig. 9) starts to be noticeable at a temperature of about $`T=0.2`$. Therefore, we can conclude the the lattice dynamics in the $`\sigma `$ phase is harmonic in a wide range of temperatures $`T0.2`$.

Finally, we would like to note the similarity between the high-temperature VDOS for the $`\sigma `$ phase and the low-temperature VDOS for the IC glass (see Fig. 10). Since $`g(\omega )`$ for the glass is only slightly temperature dependent, we show it only for $`T=0.01`$. The fact that the densities of states for the glass and the high-temperature $`\sigma `$ phase in fig. 10 are remarkably similar clearly indicates that the effect of the thermally-induced dynamical disorder in the crystalline structure on the vibrational spectrum is similar to that of the configurational disorder characteristic of the amorphous structure.

## V Conclusions

In this paper we have studied the structural and vibrational properties of a $`\sigma `$-phase crystal. First, we have shown that it is possible to construct a structural model of a one-component $`\sigma `$ phase by means of molecular dynamics simulations using an appropriate pair potential. This $`\sigma `$-phase structure is stable in a wide range of thermodynamical parameters. Our model of the $`\sigma `$ phase contains only one atomic component. This is important in understanding the role of topological icosahedral order alone on the structural and dynamical properties and avoids the effects arising from the presence of different atomic species.

Second, we have investigated atomic vibrational dynamics of the $`\sigma `$ phase. In particular, we have found the range of applicability of the harmonic approximation in a description of atomic dynamics. We have also demonstrated the existence of soft modes in the $`\sigma `$ phase which leads to a structural phase transformation with increasing pressure.

Third, we have demonstrated that the $`\sigma `$ phase is a good crystalline counterpart of the IC glass. This has been done on the basis of a comparative analysis of the vibrational dynamics (vibrational density of states).

We think that the results on the vibrational properties of the $`\sigma `$ phase discussed above can be used in an analysis of the peculiar vibrational properties of the IC glass (e.g. the Boson peak ). We also believe that the computational data of the vibrational properties of the $`\sigma `$ phase could be of value for metallurgy where this phase has received much detailed attention, chiefly because of the detrimental effect which the formation of this phase has on mechanical properties of certain steels .

## Acknowledgements

S.I.S. and M.D. thank Trinity College for hospitality. We are grateful to H.R. Schober for bringing to our attention Ref. .


text:
# Dynamics of Atom-Mediated Photon-Photon Scattering II: Experiment

## I Introduction

As described in Part I , atom-mediated photon-photon scattering is the microscopic process underlying the optical Kerr nonlinearity in atomic media. The Kerr nonlinearity produces such effects as self-phase modulation, self-focusing and self-defocusing and four-wave mixing. In an atomic medium, resonanant nonlinearities can give rise to very large nonlinear optical effects, suggesting the possibility of nonlinear optical interactions with only a few photons . Unfortunately, the same resonances which could facilitate such experiments make them difficult to analyze . In Part I we showed theoretically that the photon-photon interaction is not intrinsically lossy, and can be fast on the time scale of atomic relaxation. Here we describe an experiment to directly measure the time-duration of the photon-photon interaction in a transparent medium.

In the scattering experiment, two off-resonance laser beams collide in a rubidium vapor cell and scattering products are detected at right angles. The process of phase-matched resonance fluorescence in this geometry has been described as spontaneous four-wave mixing , a description which applies to our off-resonant excitation as well. This geometry has been of interest in quantum optics for generating phase-conjugate reflection . Elegant experiments with a barium atomic beam  showed antibunching in multi-atom resonance fluorescence, but a separation of timescales was not possible since the detuning, linewidth and doppler width were all of comparable magnitude.

## II Setup

A free-running 30mW diode laser at 780nm was temperature stabilized and actively locked to the point of minimum fluorescence between the hyperfine-split resonances of the D2 line of rubidium. Saturation spectroscopy features could be observed using this laser, indicating a linewidth $`\delta \nu <200`$ MHz. This linewidth is small compared with the detuning from the nearest absorption line $`\delta \nu =1.3`$ GHz . Direct observation of the laser output with a fast photodiode (3 dB rolloff at 9 GHz) showed no significant modulation in the frequency band $`100`$ MHz – $`2`$ GHz. The laser beam was was shaped by passage through a single-mode polarization-maintaining fiber, collimated and passed through a scattering cell to a retro-reflection mirror. The beam within the cell was linearly polarized in the vertical direction. The beam waist (at the retroreflection mirror) was $`0.026`$ cm $`\times 0.023`$ cm (intensity FWHM, vertical $`\times `$ horizontal). The center of the cell was 1.9 cm from the retroreflection mirror, thus within a Rayleigh range of the waist. With optimal alignment, the laser could deliver 1.95 mW to the cell, giving a maximal Rabi frequency of $`\mathrm{\Omega }_{\mathrm{Rabi}}2\times 10^9`$s-1, significantly less than the minimal detuning of $`\delta =2\pi \times 1.3`$ GHz $`=8\times 10^9`$s-1. For this reason, we have neglected saturation of the transitions in the analysis.

The retro-reflected beam returned through the fiber and was picked off by a beamsplitter. The single-mode fiber acted as a near-ideal spatial filter and the returned power through the fiber provides a quantitative measure of the mode fidelity on passing through the rubidium cell. With optimal alignment it was possible to achieve a mode fidelity (described below) of 36%.

The cell, an evacuated cuvette filled with natural abundance rubidium vapor, was maintained at a temperature of 330 K to produce a density of about 1.6 $`\times 10^{10}`$ cm-3. Irises near the cell limited the field of view of the detectors. Stray light reaching the detectors was negligible, as were the detectors’ dark count rates of $`<100`$ cps.

With the aide of an auxiliary laser beam, two single-photon counting modules (SPCMs) were positioned to detect photons leaving the detection region in opposite directions. In particular, photons scattered at right-angles to the incident beams and in the direction perpendicular to the drive beam polarization were observed. Each detector had a 500 $`\mu `$m diameter active area and a quantum efficiency of about 70%. The detectors were at a distance of 70 cm from the center of the cell. The effective position of one detector could be scanned in two dimensions by displacing the alignment mirrors with inchworm motors. A time-to-amplitude converter and multichannel analyzer were used to record the time-delay spectrum. The system time response was measured using sub-picosecond pulses at 850nm as an impulse source. The response was well described by a Lorentzian of width 810 ps (FWHM).

Optimal alignment of the laser beam to the input fiber coupler could not always be maintained against thermal drifts in the laboratory. This affected the power of the drive beams in the cell but not their alignment or beam shape. These were preserved by the mode-filtering of the fiber. Since the shape of the correlation function depends on beam shape and laser tuning but not on beam power, this reduction in drive power reduced the data rate but did not introduce errors into the correlation signal.

### A Experimental Results

The time-delay spectrum of a data run of 45 hours is shown in Fig. 2. The detectors were placed to collect back-to-back scattering products to maximize the photon-photon scattering signal. A Gaussian function $`P(t_At_B)`$ fitted to the data has a contrast $`[P(0)P(\mathrm{})]/P(\mathrm{})`$ of $`0.046\pm 0.008`$, a FWHM of $`1.3\pm 0.3`$ ns, and a center of $`0.07\pm 0.11`$ ns. This center position is consistent with zero, as one would expect by the symmetry of the scattering process. For comparison, a reference spectrum is shown. This was taken under the same conditions but with one detector intentonally misaligned by much more than the angular width of the scattering signal.

The angular dependence of the scattering signal was investigated by acquiring time-delay spectra as a function of detector position. To avoid drifts over the week-long acquision, the detector was scanned in a raster pattern, remaining on each point for 300 s before shifting to the next. Points were spaced at 1 mm intervals. Total live acquisition time was 9 hours per point. The aggregate time-spectrum from each location was fitted to a Gaussian function with fixed width and center determined from the data of Fig. 2. The position-dependent contrast C(x,y) is shown in Fig. 3. A negative value for the contrast means that the best fit had a coincidence dip rather than a coincidence peak at zero time. These negative values are not statistically significant. Fitted to a Gaussian function, C(x,y) has a peak of 0.044 $`\pm .010`$ and angular widths (FWHM) of $`1.1\pm 0.7`$ mrad and $`3.7\pm 0.4`$ mrad in the horizontal and vertical directions, respectively.

These angular widths are consistent with the expected coherence of scattering products . Seen from the detector positions, the excitation beam is narrow in the vertical direction, with a Gaussian shape of beam waist $`w_y=0.009`$ cm, but is limited in the horizontal direction only by the apertures, of size $`\mathrm{\Delta }z=0.08`$ cm. Thus we expect angular widths of $`0.9\mathrm{mrad}`$ and $`3.25\mathrm{mrad}`$, where the first describes diffraction of a Gaussian, the second diffraction from a hard aperture.

## III Comparison to theory

The correlation signal predicted by the theory of Part I is shown in Fig. 4. The ideal contrast is 1.53 and the FWHM is 870 ps. The shape of the time correlations is altered by experimental limitations. First, beam distortion in passing through the cell windows reduces the photon-photon scattering signal. Second, finite detector response time and finite detector size act to disperse the signal. None of these effects alters the incoherent scattering background.

Beam distortion is quantified by the fidelity factor introduced in Part I

$$F4\frac{\left|d^3xG(𝐱)H(𝐱)\right|^2}{\left[d^3x\left(|G(𝐱)|^2+|H(𝐱)|^2\right)\right]^2}$$
(1)
The greatest contrast occurs when $`H`$ is the phase-conjugate, or time-reverse of $`G`$, i.e., when $`H(𝐱)=G^{}(𝐱)`$. In this situation $`F=1`$. Under the approximation that the field envelopes obey the paraxial wave equations

$`{\displaystyle \frac{d}{dz}}G`$ $`=`$ $`{\displaystyle \frac{i}{2k}}_{}^2G`$ (2)
$`{\displaystyle \frac{d}{dz}}H`$ $`=`$ $`{\displaystyle \frac{i}{2k}}_{}^2H,`$ (3)
Green’s theorem can be used to show that the volume integral is proportional to the mode-overlap integral

$$d^3xG(𝐱)H(𝐱)=\mathrm{\Delta }z𝑑x𝑑yG(𝐱)H(𝐱),$$
(4)
where the last integration is taken at any fixed $`z`$ and $`\mathrm{\Delta }z`$ is the length of the interaction region. Similarly, the beam powers are invariant under propagation and the mode fidelity can be expressed entirely in terms of surface integrals as

$`F`$ $`=`$ $`4\left|{\displaystyle 𝑑x𝑑yG(𝐱)H(𝐱)}\right|^2`$ (6)
$`\times \left[{\displaystyle 𝑑x𝑑y\left(|G(𝐱)|^2+|H(𝐱)|^2\right)}\right]^2.`$
The overlap of $`G`$ and $`H`$ also determines the efficiency of coupling back into the fiber. This allows us to determine $`F`$. In terms of $`P_{\mathrm{in}}`$, the power leaving the output fiber coupler and $`P_{\mathrm{ret}}`$, the power returned through the fiber after being retro-reflected, this is

$$F=\frac{4}{(1+T^2)^2}\frac{P_{\mathrm{ret}}}{\eta TP_{\mathrm{in}}}.$$
(7)
where $`\eta =0.883`$ is the intensity transmission coefficient of the fiber and coupling lenses and $`T=0.92`$ is the transmission coefficient for a single-pass through a cell window. We find $`F=0.36\pm 0.03`$. The mode fidelity acts twice to reduce contrast, once as the drive beams enter the cell, and again on the photons leaving the cell. This beam distortion has no effect on the incoherent scattering background, thus the visibility is reduced by $`F^2`$.

The finite time response of the detector system acts to disperse the coincidence signal over a larger time window. This reduces the maximum contrast by a factor of 0.27 and increases the temporal width to 1.62 ns. Similarly, the finite detector area reduces the maximum contrast by a factor of 0.81 and spreads the angular correlations by a small amount. The resulting coincidence signal is shown in Fig. 6. Fitted to a Gaussian, the final signal contrast is 0.042 $`\pm 0.007`$, where the uncertainty reflects the uncertainty in $`F`$. This is consistent with the observed contrast of 0.044 $`\pm `$ 0.010.

## IV conclusion

We have measured the temporal and angular correlations in photon-photon scattering mediated by atomic rubidium vapor. We found good agreement between experiment and the perturbative theory presented in Part I. The observed temporal correlations are of the order of one nanosecond, much faster than the system can relax by radiative processes. This is consistent with the prediction that the duration of the photon-photon interaction is determined by the inhomogeneous broadening of the vapor.


text:
# References

EXACT SOLUTIONS OF THE DIRAC EQUATION

FOR MODIFIED COULOMBIC POTENTIALS

ANTONIO SOARES DE CASTROa AND JERROLD FRANKLINb

aDFQ, UNESP, C.P. 205, 12500-000 Guaratinguetá SP, Brasil
bDepartment of Physics, Temple University, Philadelphia, PA 19122-6082
Received February 2000

> Exact solutions are found to the Dirac equation for a combination of Lorentz scalar and vector Coulombic potentials with additional non-Coulombic parts. An appropriate linear combination of Lorentz scalar and vector non-Coulombic potentials, with the scalar part dominating, can be chosen to give exact analytic Dirac wave functions. The method works for the ground state or for the lowest orbital state with $`l=j\frac{1}{2}`$, for any $`j`$.

In a previous letter, simple exact solutions were found for the Dirac equation for the combination of a Lorentz vector Coulomb potential with a linear confining potential that was a particular combination of Lorentz scalar and vector parts. In this letter, we extend the method of Ref. to the more general case of an arbitrary combination of Lorentz scalar and vector Coulombic potentials with a particular combination of Lorentz scalar and vector non-Coulombic potentials, $`S(r)`$ and $`V(r)`$. The non-Coulombic potentials can have arbitrary radial dependence, but must be related by

$$V(r)=\frac{E}{m}S(r),$$
(1)
where $`E`$ is the bound state energy for a particle of mass $`m`$. This is the same relation between the vector and scalar non-Coulombic potentials as was required in Ref.

The Dirac equation we solve is

$$\left[\alpha 𝐩+\beta m\frac{(\lambda +\beta \eta )}{r}+V(r)+\beta S(r)\right]\psi =E\psi ,$$
(2)
where $`𝜶`$ and $`\beta `$ are the usual Dirac matrices. The four component wave function $`\psi `$ can be written in terms of two component spinors $`u`$ and $`v`$ as

$$\psi =N\left(\begin{array}{c}u\\ v\end{array}\right),$$
(3)
with $`N`$ an appropriate normalization constant. The two component spinors $`u`$ and $`v`$ satisfy the equations

$`(\sigma 𝐩)v`$ $`=`$ $`\left[Em+{\displaystyle \frac{(\lambda +\eta )}{r}}V(r)S(r)\right]u`$ (4)
$`(\sigma 𝐩)u`$ $`=`$ $`\left[E+m+{\displaystyle \frac{(\lambda \eta )}{r}}V(r)+S(r)\right]v.`$ (5)
The key step in generating relatively simple exact solutions of the Dirac equation is to choose a particularly simple form for the function $`v(r)`$

$$v(r)=i\gamma (\sigma \widehat{𝐫})u(r),$$
(6)
where $`\gamma `$ is a constant factor, to be determined by the solution to the Dirac equation. This is the form of $`v(r)`$ that was found in Ref. for a Coulomb plus linear confining potential. This ansatz for $`v(r)`$ has also been used as the basis for generating approximate saddle point solutions for the Dirac and Breit equations.

For now, we limit our discussion to spherically symmetric states $`u(r)`$ that depend only on the radial coordinate. We will discuss orbitally excited states later in this paper. Using the form of $`v(r)`$ given by Eq. (6), equations (4) and (5) reduce to two first order ordinary diffential equations for $`u(r)`$

$`{\displaystyle \frac{du}{dr}}`$ $`=`$ $`{\displaystyle \frac{1}{\gamma }}\left[EmV(r)S(r)+{\displaystyle \frac{(\lambda +\eta 2\gamma )}{r}}\right]u(r)`$ (7)
$`{\displaystyle \frac{du}{dr}}`$ $`=`$ $`\gamma \left[E+mV(r)+S(r)+{\displaystyle \frac{(\lambda \eta )}{r}}\right]u(r).`$ (8)
Equations (7) and (8) are two independent equations for the same quantity, so that each term in one equation can be equated with the corresponding term in the other equation having the same radial dependence. This leads to

$$\gamma ^2=\frac{mE}{m+E}=\frac{S(r)+V(r)}{S(r)V(r)}=\frac{\lambda +\eta 2\gamma }{\eta \lambda }.$$
(9)
The relations in Eq. (9) can be rearranged, after some algebra, to give

$$\gamma =\frac{\lambda +\eta }{1+b},$$
(10)
with

$$b=\pm \sqrt{1\lambda ^2+\eta ^2}.$$
(11)
The constant b can have either sign. Although $`b`$ must be positive in the pure Coulombic case, we will see that a negative $`b`$ is possible if the Lorentz scalar potential $`S(r)`$ is more singular at the origin than $`1/r`$. The bound state energy can be written as

$$E=m\left(\frac{1\gamma ^2}{1+\gamma ^2}\right)=m\left(\frac{b\lambda \eta }{\lambda b\eta }\right),$$
(12)
or as

$$E=m\frac{V(r)}{S(r)}.$$
(13)
The wave function $`u(r)`$ can be found by solving differential equation (8) to give

$$u(r)=r^{b1}\mathrm{exp}\left[a\left(r+\frac{1}{m}S(r)𝑑r\right)\right],$$
(14)
where the constant a is given by

$$a=\gamma (m+E)=\pm \sqrt{m^2E^2}.$$
(15)
The constants $`\gamma `$ and $`a`$ could be negative if $`S(r)`$ approaches a large enough negative constant or diverges negatively as $`r`$ becomes infinite. The integral in equation (14) can diverge at the origin or as $`r\mathrm{}`$, or for any finite $`r`$, as long as the quantity in square brackets in Eq. (14) remains negative.

The constraint equation (13) shows that in order for this class of exact solutions to apply, the Lorentz vector and scalar non-Coulombic potentials must have the same radial dependence. As long as this constraint is satisfied, the results in Eqs. (10)-(15) represent a complete exact solution for the wave function and bound state energy of the Dirac Hamiltonian given in equation (2).

We note from Eq. (12) that the energy seems to depend only on the Coulombic coupling constants $`\lambda `$ and $`\eta `$, and does not seem to depend on the non-Coulombic potentials. However, we could alternatively say from Eq. (13) that the energy depends only on the ratio of the non-Coulombic potentials $`V(r)`$ and $`S(r)`$, and does not seem to depend on the Coulombic coupling constants. Then, Eq. (12) could be considered a constraint equation on the Coulombic coupling constants. The actual situation is that the energy does depend on all the potentials but, because of the severe constraints imposed by equations (12) and (13) taken together, the energy can be written in terms of one set of potentials or the other. Although the possibility of this class of exact solutions is limited by the constraints on the potentials, this still permits a wide range of non-Coulombic potentials.

We now consider conditions imposed on the potentials and the wave function parameters by the physical requirements that the potentials be real and the wave function normalizable. We see from Eq. (9) that $`\gamma `$ must be real, and then from Eq. (10) that $`b`$ must be real. This requires the Coulombic potentials to satisfy the condition

$$1\lambda ^2+\eta ^20.$$
(16)
The reality of $`\gamma `$ restricts possible bound state energies to the range

$$m<E<m.$$
(17)
Note that negative energies can occur, but $`E+m`$ cannot be negative. Also, $`E`$ cannot equal $`\pm m`$, because this would lead to an unormalizable wave function. This condition on $`E`$, along with Eq. (13), means that $`V(r)`$ must always be less in magnitude than $`S(r)`$.

We discuss the remaining conditions on the parameters in terms of three sub-classes of solution:

1. The ”normal” class of solutions has $`b`$, $`\gamma `$, and $`a`$ all positive. In this case, we see from Eq. (10) that the Coulombic potentials must satisfy the further condition

$$\lambda +\eta >0.$$
(18)
2. Sub-case 2 has $`b`$ negative, with $`\gamma `$ and $`a`$ still positive. The constant $`b`$ can be negative if the product $`aS(r)`$ is positively divergent at the origin faster than $`1/r`$. Then each of Eqs. (10)-(15) holds just as for positive $`b`$, and the wave function is still normalizable. Sub-case 1 with positive $`b`$ transforms smoothly into the pure Coulombic solution as the non-Coulombic potential tends to zero everywhere. But this is not true for sub-case 2 with $`b`$ negative. This sub-case requires the non-Coulombic potential to be dominant at the origin, and so has no corresponding pure Coulombic limit.
3. Sub-case 3 has a negative $`\gamma `$ and a negative $`a`$, while $`b`$ can have either sign, as discussed in sub-cases 1 and 2 above. A negative $`a`$ is possible if the non-Coulombic potential diverges or approaches a constant as $`r\mathrm{}`$, so that the integral in Eq. (15) diverges faster than $`r`$ at large $`r`$. Since $`a`$ is negative, the potential $`S(r)`$ must be negative at large $`r`$. Then all of Eqs. (10)-(15) hold as for positive $`a`$, and the wave function is still normalizable. This case is highly unusual, because it allows the possibility of a potential that is negative everywhere and diverges negatively at both the origin and infinite $`r`$. We know of no other example in quantum mechanics where a potential that diverges negatively at infinity can lead to a normalizable ground state. The reason this is possible here can be seen from Eqs. (7) and (8). There it is seen that $`S(r)`$ enters the differential equations for $`u(r)`$ only in the combinations $`\gamma S(r)`$ or $`S(r)/\gamma `$. Since these effective potentials are positive, the resulting wave function is normalizable. As with sub-case 2, the case with $`\gamma `$ and $`a`$ negative cannot not approach a pure Coulombic case because the non-Coulombic potential must be dominant at large r.

We now look at some special cases. If the non-Coulombic potentials are absent, then the solutions are for a general linear combination of Lorentz vector and Lorentz scalar Coulombic potentials. If either constant, $`\lambda `$ or $`\eta `$, is zero, we recover the usual solutions of the Dirac equation for a pure scalar or vector Coulombic potential. The Coulombic potentials cannot both be absent (while keeping a non-Coulombic part) because then $`\gamma `$ would be zero and b one, leading to a constant, unnormalizable wave fucntion.

For a power law non-Coulombic potential of the form

$$S(r)=\mu (s+1)r^s,s1,$$
(19)
the wave function will be given by

$$u(r)=r^{b1}exp\left[a\left(r\frac{\mu }{m}r^{s+1}\right)\right].$$
(20)
All three sub-classes are possible for this wave function, depending on the ranges of the parameters $`b`$, $`a`$, $`\mu `$, and $`s`$.

The method does not work for radially excited states, because the simple ansatz of Eq. (6) for $`v(r)`$ does not lead to consistent equations for $`du/dr`$ in that case. But, the method does work for the lowest orbitally excited state for which $`l=j\frac{1}{2}`$. This can be seen by writing the Dirac wave function in terms of radial functions $`f(r)`$ and $`g(r)`$

$$\psi =N\left(\begin{array}{c}f(r)𝒴_{j\frac{1}{2}}^j\\ ig(r)𝒴_{j+\frac{1}{2}}^j\end{array}\right),$$
(21)
where $`𝒴_l^j`$ is a two component angular spinor function corresponding to total angular momentum $`j`$ and orbital momentum $`l`$. The radial functions satisfy the equations

$`\left({\displaystyle \frac{d}{dr}}+{\displaystyle \frac{1\kappa }{r}}\right)f`$ $`=`$ $`\left[E+m+{\displaystyle \frac{(\lambda \eta )}{r}}+S(r)V(r)\right]g`$ (22)
$`\left({\displaystyle \frac{d}{dr}}+{\displaystyle \frac{1+\kappa }{r}}\right)g`$ $`=`$ $`\left[Em+{\displaystyle \frac{(\lambda +\eta )}{r}}S(r)V(r)\right]f,`$ (23)
where $`\kappa =j+\frac{1}{2}`$ is the principal quantum number of the state.

For the radial function $`g(r)`$, we make the ansatz

$$g(r)=\gamma f(r).$$
(24)
Then, equations (22) and (23) reduce to two first order differential equations for $`f(r)`$

$`{\displaystyle \frac{df}{dr}}`$ $`=`$ $`\gamma \left[E+mV(r)+S(r)+{\displaystyle \frac{(\lambda \eta )}{r}}{\displaystyle \frac{(\kappa 1)}{\gamma r}}\right]f(r)`$ (25)
$`{\displaystyle \frac{df}{dr}}`$ $`=`$ $`{\displaystyle \frac{1}{\gamma }}\left[EmV(r)S(r)+{\displaystyle \frac{(\lambda +\eta )}{r}}{\displaystyle \frac{\gamma (\kappa +1)}{r}}\right]f(r).`$ (26)
Equating the corresponding terms having the same radial dependence in these two equations results in

$$\gamma ^2=\frac{mE}{m+E}=\frac{S(r)+V(r)}{S(r)V(r)}=\frac{\lambda +\eta 2\gamma \kappa }{\eta \lambda }.$$
(27)
Equation (27) is the same as Eq. (9) with the replacements

$$\lambda \lambda /\kappa ,\eta \eta /\kappa .$$
(28)
Thus, equations (10)-(13), and (15) hold for the orbitally excited case, with the replacements $`\lambda \lambda /\kappa `$ and $`\eta \eta /\kappa `$. The radial wave function is $`f(r)`$ is given by

$$f(r)=r^{\kappa b1}exp\left[a\left(r+\frac{1}{m}S(r)𝑑r\right)\right].$$
(29)
Note that the constant $`b`$ used in this paper differs from the $`b`$ defined by Eq. (22) of Ref. 1. The $`b`$ in this paper is the $`b`$ in Ref. 1 divided by $`\kappa `$.

In summary, we have presented a class of potentials for which the Dirac equation has relatively simple exact analytical solutions for the ground state of each angular momentum state with $`l=j\frac{1}{2}`$.


text:
# An extended multi-zone model for the MCG-6-30-15 warm absorber

## 1 INTRODUCTION

X-ray absorption by partially ionized, optically thin material, along the line of sight of the central engine, the so called warm absorber, is one of the prominent features in the X-ray spectrum of many AGN . The presence of such gas was first postulated in order to explain the unusual form of the X-ray spectrum of QSO MR 2251-178 . A direct confirmation of the existence of circumnuclear ionized matter came from ASCA, which for the first time was able to resolve the OVII and OVIII absorption edges (at rest energies of 0.74 and 0.87 keV , respectively) in the X-ray spectra of the Seyfert 1 galaxy MCG$``$6-30-15 . Systematic studies of warm absorbers in Seyfert 1 galaxies with ASCA have shown their ubiquity, being detected in half of the sources .

Warm absorbers were usually described by single zone, photoionization equilibrium models. Under this assumption, simple variability patterns were expected: increasing ionization of the matter when the source brightens. This was not found by Otani et al.  during their long-look ASCA observation of the nearby (z=0.008) Seyfert 1 galaxy MCG$``$6-30-15. The source varied with large amplitude on short time scales ($`10^4`$ s or so). The depth of the OVII edge, on the contrary, stayed almost constant, while that of the OVIII edge was anticorrelated with the flux. To explain the behaviour of the OVII and OVIII edges, the authors adopted a multizone model in which the OVII and OVIII edges originate from spatially distinct regions. OVII ions in the region responsible for the OVII edge, the outer absorber, have a long recombination timescale (i.e. weeks or more), whereas the OVIII edge arises from more highly ionized material, the inner absorber, in which most oxygen is fully stripped. The recombination timescale for the OIX ions in the inner absorber is of the order of $`10^4s`$ or less. A decrease in the primary ionizing flux is then accompanied by the recombination of OIX to OVIII, giving the observed variation in the OVIII edge depth.

Orr et al.  raised the possibility of a more complex warm absorber. During their MCG$``$6-30-15 BeppoSAX observation the depth of the OVII edge was marginally consistent with being constant, whereas the optical depth for OVIII, $`\tau (OVIII)`$, exhibited significant variability. The authors claimed that its large value during epoch 1111The epochs in the BeppoSAX observation are chronologically numerated (i.e. number 1 corresponds to the first epoch of the observation, etc.). ($`1.7\pm 0.5`$, $`1\sigma `$ uncertainty) was inconsistent with the values at all other epochs. They also pointed out that the value of $`\tau (OVIII)`$ during epoch 5 (a low luminosity state) did not show the expected anticorrelaction with the ionizing flux .

In both Otani et al.  and Orr et al., $`\tau (OVIII)`$ was plotted versus count rate. In order to compare ASCA and BeppoSAX observations, these two observational results for $`\tau (OVIII)`$ have been plotted versus luminosity222The conversion from count rate to luminosity has been obtained using $`H_o=50`$ $`\mathrm{km}\mathrm{s}^1\mathrm{Mpc}^1`$, $`q_0=0`$ and PIMMS (http://heasarc.gsfc.nasa.gov/Tools/w3pimms.html). in 0.1-10 keV in figure 1 (instead of versus count rate). The conversion from count rate to luminosity is different for each apparatus and therefore, the comparison can not be made using count rate, but luminosity.

The only discrepancy between both sets of data is the result for epoch 1 (labelled 1 in figure 1)333This epoch was at the very beginning of the BeppoSAX observation. In a private communication, A. Orr notes that this high value for the optical depth for OVIII could not be due to instrumental effects, since no similar behaviour, neither at the beginning nor during the observation has been observed in any other BeppoSAX target yet.. Note also that there is no disagreement at all with any other point (Orr et al. mentioned epoch 5 as problematic, but in this plot it appears consistent with the rest of the Otani et al. data).

In this paper we present a simple photoionization model that accounts for the experimental results of both ASCA and BeppoSAX observations. Section 2 describes the code used to model the inner warm absorber. The application to the MCG$``$6-30-15 warm absorber is presented in Section 3. The extension of the multi-zone model for the warm absorber (i. e. the third component) is addressed in Section 4. Finally our conclusions are discussed in Section 5.

## 2 TIME DEPENDENT PHOTOIONIZATION CODE

In order to reproduce both ASCA and BeppoSAX observations, the state of the inner absorber has been modelled using a time dependent photoionization code for oxygen . This code treats the material in the inner absorber as containing only ions of oxygen in an otherwise completely ionized hydrogen plasma. A given oxygen ion is assumed to be in one of the 9 states corresponding to its ionization level. Different excitation levels of a given ionization state are not treated. The total abundance of oxygen is fixed at $`7.41\times 10^4`$ relative to hydrogen. Both ions and electrons are assumed to be in local thermal equilibrium (LTE) with a common temperature T and the plasma is assumed to be strictly optically-thin at all frequencies (i.e. radiative transfer processes are not considered).

A point source of ionizing radiation is placed at a distance $`R`$ from this material with a frequency dependent luminosity $`L_\nu `$ (total luminosity $`L`$). In the calculations presented here, the ionizing continuum is taken to be a power-law with a photon index $`\mathrm{\Gamma }=2`$ extending from $`\nu _{min}`$ to $`\nu _{max}`$: i.e. we take

$$L_\nu =\frac{L}{\mathrm{\Lambda }\nu }$$
(1)
where $`\mathrm{\Lambda }=\mathrm{ln}(\nu _{max}/\nu _{min})`$. The upper and lower cutoff frequencies are chosen such that $`h\nu _{max}=100`$ keV , and $`\mathrm{\Lambda }=10`$, corresponding to $`h\nu _{min}4.5`$ eV .

The physical processes included in the code are: photoionization, Auger ionization, collisional ionization, radiative recombination, dielectronic recombination, bremsstrahlung, Compton scattering and collisionally-excited O$`\lambda `$1035 resonance line cooling. Given these processes, the ionization and thermal structure of the plasma is evolved from an initial state assumed to have a temperature of $`T=10^5`$ K , and equal ionization fractions in each state.

The ionization structure of the oxygen is governed by 9 equations of the form:

$$\frac{dn_i}{dt}=\underset{prod.}{}\delta _{prod.}\underset{dest.}{}\delta _{dest.}$$
(2)
where $`n_i`$ is the number density of state-i, $`\delta `$ stands for the ionization/recombination rate per unit volume, the first summation on the right hand side (RHS) is over all ionization and recombination mechanisms that produce state-i and the second summation is over all such mechanisms that destroy state-i.

The thermal evolution of the plasma is determined by the local heating/cooling rate and the macroscopic constraint. The net heating rate per unit volume is given by:

$$\frac{dQ}{dt}=\underset{heat}{}\mathrm{\Delta }_{heat}\underset{cool}{}\mathrm{\Delta }_{cool}$$
(3)
where $`\mathrm{\Delta }`$ stands for the heating/cooling rate per unit volume, the first summation on the RHS is over all heating processes and the second term is over all cooling processes. For a constant density plasma, we have $`du=dQ`$ where $`u`$ is the internal energy per unit volume and is given by:

$$u=\frac{3}{2}n_ek_BT.$$
(4)
Thus, the rate of change of temperature is related to the heating/cooling rate via:

$$\frac{dT}{dt}=\frac{2}{3n_ek_B}\frac{dQ}{dt}.$$
(5)
The state of the plasma is evolved using the ionization equations represented by (2) and the energy equation (3). A simple step-by-step integration of these differential equation is performed. The time step for the integration process is allowed to dynamically change and is set to be $`0.1\times t_{sh}`$ where $`t_{sh}`$ is the shortest relevant timescale444At each step the ionization, collisional and recombination timescales of each oxygen ionization state is calculated. The shortest of these timescales is taken to be $`t_{sh}`$.. In the case of a constant ionizing luminosity, it is found that the system always evolves to an equilibrium state, as expected.

The model has been compared against the photoionization code cloudy (version 9004, ) in the case of a plasma in photoionization equilibrium. Figures 2 and 3 report this comparison.There is a qualitative agreement between the model (solid line) and cloudy (dashed line) in figure 2 for the case of a pure hydrogen/oxygen plasma. The main goal of this paper is to reproduce the time variability of the oxygen edges, and, this figure shows how our model mirrors quite well the behavior of this element. cloudy includes many more physical processes and realistic elemental abundances and, as seen in figure 3, the agreement is not so good when more elements are considered555Realistic elemental abundances can be important both for the ionization structure and for the heating/cooling rate of the plasma, and since features of the stability curve result from peaks in the heating and cooling functions (which are in turn associated with particular ionic species), the comparison presented in figure 2 is expected to present a much better agreement. The consideration of realistic abundances is beyond the scope of this paper and is currently being studied..

We have also investigated the reaction of a constant density plasma to an inverted top hat function for the light-curve (i.e. initially $`T=10^5`$ K , equal ionization fractions in each state, the luminosity is $`L=3\times 10^{43}`$ $`\mathrm{erg}\mathrm{s}^1`$and held constant. The system is allowed to achieve an equilibrium at an ionization parameter $`\xi =\frac{L}{nR^2}=75`$ $`\mathrm{erg}\mathrm{cm}\mathrm{s}^1`$, where $`n`$ is the density of the warm plasma and $`R`$ is the distance of the slab of the plasma from the ionizing source of radiation with isotropic ionizing luminosity $`L`$. The luminosity is then halved for a short period, $`\mathrm{\Delta }t`$, and then back to its initial value) . During these changes, the state of the plasma is recorded as a function of time. This choice of $`\xi `$ places the inner warm absorber a region where $`f_{O8}L^1`$ (i.e. in this regime $`f_{O8}`$ is expected to be inversely proportional to $`\xi `$ or equivalently to $`L`$). The investigation has been carried out for different values of $`n`$ and $`R`$ keeping $`n>2\times 10^7`$ $`\mathrm{cm}^3`$and $`R<1.4\times 10^{17}`$ cm (i.e. the constraints given by Otani et al. for the inner warm absorber). Fig. 4 shows an example of the results obtained. As expected, $`f_{O8}`$ increases when $`L`$ is halved, showing a more noticeable increment when $`L`$ is halved over a longer period.

## 3 APPLICATION TO THE MCG$``$6-30-15 WARM ABSORBER

To reproduce the observational data using our model, it is required to obtain the value of the optical depth for OVIII once the ionization fraction $`f_{O8}`$ has been calculated by the code. Considering the cross section for OVIII, $`\sigma _8`$, as constant along the line-of-sight and using the fact that the material contains only ions of oxygen in an otherwise completely ionized plasma, the optical depth for OVIII, $`\tau (OVIII)`$, can be written as:

$$\tau (OVIII)=\frac{\sigma _8f_{O8}\mathrm{\Delta }Rn_e}{\frac{1}{Abun}+\underset{i=1}{\overset{9}{}}f_{Oi}\times (i1)}$$
(6)
where $`\sigma _8=10^{19}cm^2`$ , $`\mathrm{\Delta }R`$ is the line-of-sight distance through the ionized plasma, $`n_e`$ is the electron density, $`Abun`$ is the oxygen abundance relative to hydrogen and $`f_{Oi}`$ is the ionization fraction for i ionization species for oxygen.

The light-curve used to reproduced the ASCA observation is given by Otani et al. . For the BeppoSAX observation, the light-curve has been defined as a step function with a constant value for the luminosity $`L`$ over the different time periods given by Orr et al. . The constant value for $`L`$ over each period is chosen to be that given in Orr et al. (1997, Fig. 3) plus a period of $`2\times 10^4`$ s previous to epoch 1 with $`L=1.8\times 10^{43}`$ $`\mathrm{erg}\mathrm{s}^1`$(a standard value for the luminosity during the observation). Different parameters for the warm absorber have been investigated using the constraints given by Otani et al. . With the only exception of point 1 (see figure 1), a general good agreement for all other experimental points is found.

Following a suggestion by Orr et al. , we have modified the light-curve for the BeppoSAX observation including a previous epoch to it for which the luminosity is much lower (the range of values used is $`10^{40}10^{43}`$ $`\mathrm{erg}\mathrm{s}^1`$)666An example of a Seyfert 1 galaxy exhibiting an unusual low flux state (a decrease of more than one order of magnitude in luminosity) is presented in Uttley et al. .. The reason for using this very low value for the luminosity is that in a regime for which $`f_{O8}L^1`$, the highest values $`f_{O8}`$ are expected for low values of $`L`$. However, the suggested explanation does not reproduce the high value of $`\tau (OVIII)`$ for epoch 1 in the BeppoSAX observation. Even for those values of $`\xi `$ that give a maximum for $`f_{O8}`$ (i.e. $`\xi 50`$ $`\mathrm{erg}\mathrm{cm}\mathrm{s}^1`$), the ionization fraction for O(VIII) is still too low to account for the high $`\tau (OVIII)`$ value at epoch 1.

## 4 A THIRD WARM ABSORBER COMPONENT

The model we propose to explain both ASCA and BeppoSAX observations incorporates a new zone for the warm absorber. Let warm absorber 1 $``$ WA1 be the inner warm absorber which parameters are: $`R<1.4\times 10^{17}`$ cm , $`n>2\times 10^7`$ $`\mathrm{cm}^3`$and $`\mathrm{\Delta }R10^{14}`$ cm . The outer warm absorber will be warm absorber 2 $``$ WA2 with parameters $`R>3\times 10^{18}`$ cm , $`n<2\times 10^5`$ $`\mathrm{cm}^3`$and $`\mathrm{\Delta }R10^{14}`$ cm . In our model a warm absorber 3 $``$ WA3 is situated between WA1 and WA2. The WA3 radius and density will have values between those of WA1 and WA2. Therefore, while WA1 and WA2 respond on timescales of hours and weeks respectively, WA3 is expected to respond to variations in the ionizing flux on timescales of days. WA3 would have an ionization parameter $`\xi `$ of the order of 500 $`\mathrm{erg}\mathrm{cm}\mathrm{s}^1`$for $`L=3\times 10^{43}`$ $`\mathrm{erg}\mathrm{s}^1`$. Hence the ionization fraction for this value of $`L`$ is too low to be detected. Only when the luminosity of the source is sufficiently low (i.e. $`\xi 50`$ $`\mathrm{erg}\mathrm{cm}\mathrm{s}^1`$), WA3 reveals its presence by contributing to the total optical depth for OVIII. The medium between WA1 and WA2 would be constituted by a continuum of warm absorbers: clouds with different densities situated at different radii. Only some of them happen to be at the radius and have densities that efficiently absorbs the OVIII edge energy (i.e. most of them are undetectable). In Baldwin et al.  a model for the BLR is presented, in which individual BLR clouds can be thought of as machines for reprocessing radiation. As long as there are enough clouds at the correct radius and with the correct gas density to efficiently form a given line, the line will be formed with a relative strength which turns out to be very similar to the one observed. Similarly in our model, only WA1 and WA2 are detectable for the ordinary values of the ionizing flux. Only for the case of a state of low luminosity, WA3 will be unmasked. Other zones, as yet unseen, may be present.

Assuming then the presence of WA3 and also an epoch of low luminosity previous to the BeppoSAX observation, the expected WA3 behaviour would be:

i) when $`L0.4\times 10^{43}`$ $`\mathrm{erg}\mathrm{s}^1`$(previous to epoch 1), the ionization parameter $`\xi 50`$ $`\mathrm{erg}\mathrm{cm}\mathrm{s}^1`$, giving a high value for $`f_{O8}`$. This period lasts for approximately $`10^5`$ s , so the plasma has time to recombine.

ii) when $`L(1,5)\times 10^{43}`$ $`\mathrm{erg}\mathrm{s}^1`$(i.e. during the observation), then $`\xi 150,750`$ $`\mathrm{erg}\mathrm{cm}\mathrm{s}^1`$. For these high values of $`\xi `$, oxygen is practically fully stripped and therefore there is a very small contribution to the optical depth for OVIII.

The range of parameters investigated for WA3 is $`R=(2,8)\times 10^{17}`$ cm, $`n=5\times 10^5,10^7`$ $`\mathrm{cm}^3`$and $`\mathrm{\Delta }R`$ in the interval that gives a column density for WA3 approximately equal to $`3\times 10^{22}`$ $`\mathrm{cm}^3`$. After taking into account the soft X-ray absorption due to WA1, the response of WA3 has been calculated and an example of the results obtained is presented in figures 5 and 6, where the general good agreement is also extended to point 1.

WA3 has also been modelled using cloudy and we found a drop in the transmitted portion of the incident continuum at $`7.8`$ keV , of approximately $`1\%`$ (i.e. undetectable for current instruments). The coronal lines strength has also been checked using cloudy777See Sect. 3.3  for a discussion of the state of the coronal lines atomic data used in cloudy.. The ratio of the modelled to observed fluxes for WA3 are all below $`0.1\times f_c`$, where $`f_c`$ is the covering fraction888Porquet et al.  give some restrictions on the density of the WA in order to avoid producing coronal line equivalent widths larger than observed. Although WA3 presents no problems at all for any $`f_c`$, WA2 does, unless a low value of $`f_c`$ is considered. This possibility is currently being investigated..

Finally, the contribution to $`\tau `$(OVII) from each component for our model has been calculated and, as expected, we found that WA2 is the main responsible for the OVII edge (its contribution makes 96% of the total value of $`\tau `$(OVII)).

## 5 DISCUSSION

Warm absorbers have been the subject of extensive studies during the last decade. Such regions are not spatially resolved, and all the available information about their geometry is obtained from analysis of the variability of the oxygen edges. The explanation we offer for the time variability of the MCG$``$6-30-15 warm absorber during both ASCA and BeppoSAX observations does not invoke complex processes, but a very simple photoionization model together with the presence of a multi-zone warm absorber. This would be constituted by a continuum of clouds at different radii and different densities, such that only some of them contribute to the total optical depth for OVIII depending on the value of the luminosity.

As a final remark note how in our model WA3 is much more volume filling than WA1 ($`\mathrm{\Delta }R/R0.1`$ for WA3 while $`\mathrm{\Delta }R/R10^3`$ for WA1). This suggests that WA3 could be part of the inter-cloud medium of WA1.

## 6 ACKNOWLEDGMENTS

This work has been supported by PPARC and Trinity College (R.M.) and by the Royal Society (A.C.F.). C.S.R. thanks support from Hubble Fellowship grant HF-01113.01-98A. This grant was awarded by the Space Telescope Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., for NASA under contract NAS 5-26555. C.S.R. also thanks support from NASA under LTSA grant NAG5-6337.