text:
# Non Equilibrium Green’s Functions for Dummies: Introduction to the One Particle NEGF equations

## 1 Introduction

Non equilibrium Green’s function methods are regularly used to calculate current and charge densities in nanoscale (both molecular and semiconductor) conductors under bias. An overview of the theory of molecular electronics can be found in Ref.  and for semiconductor nanoscale devices see Ref. .

The aim of this text is to provide some intuitive explanations of one particle Green’s functions in a compact form together with derivations of the expressions for the current and the density matrix. It is not intended as a complete stand-alone tutorial, but rather as a complement to Ref. .

## 2 Green’s functions

Discrete Schrödinger equation:

$$H|n=E|n$$
(1)
We divide the Hamiltonian and wavefunction of the system into contact ($`H_{1,2}`$, $`|\psi _{1,2}`$) and device ($`H_d`$, $`|\psi _d`$) subspaces:

$$\left(\begin{array}{ccc}H_1& \tau _1& 0\\ \tau _1^{}& H_d& \tau _2^{}\\ 0& \tau _2& H_2\end{array}\right)\left(\begin{array}{c}|\psi _1\\ |\psi _d\\ |\psi _2\end{array}\right)=E\left(\begin{array}{c}|\psi _1\\ |\psi _d\\ |\psi _2\end{array}\right)$$
(2)
where $`\tau _{1,2}`$ describes the interaction between device and contacts. In general we have $`N`$ contacts ($`H_{1,\mathrm{},N}`$) connecting ($`\tau _{1,\mathrm{},N}`$) the device $`H_d`$ to the reservoirs. Here we will assume that the contacts are independent, i.e., there are no cross terms ($`\tau `$) between the different contacts.

We define the Green’s function111Others may (and do) use the opposite sign.:

$$\left(EH\right)G(E)=I$$
(3)
### 2.1 Why do we want to calculate the Green’s function?

The Green’s function gives the response of a system to a constant perturbation $`|v`$ in the Schrödinger equation:

$$H|\psi =E|\psi +|v$$
(4)
The response to this perturbation is:

$`\left(EH\right)|\psi `$ $`=`$ $`|v`$ (5)
$`|\psi `$ $`=`$ $`G(E)|v`$ (6)
Why do we need the response to this type of perturbation? Well, it turns out that it’s usually easier (see next section) to calculate the Green’s function than solve the whole eigenvalue problem222Especially for infinite systems. and most (all for the one-particle system) properties of the system can be calculated from the Green’s function. E.g., the wavefunction of the contact ($`|\psi _2`$) can be calculated if we know the wavefunction on the device ($`|\psi _d`$). From third row of Eq. 2:

$`H_2|\psi _2+\tau _2|\psi _d`$ $`=`$ $`E|\psi _2`$ (7)
$`\left(EH_2\right)|\psi _2`$ $`=`$ $`\tau _2|\psi _d`$ (8)
$`|\psi _2`$ $`=`$ $`g_2(E)\tau _2|\psi _d`$ (9)
where $`g_2`$ is the Green’s function of the isolated contact 2 ($`(EH_2)g_2=I`$).

It is important to note that since we have an infinite system, we obtain two types of solutions for the Green’s functions333 When the energy coincides with energy band of the contacts there are two solutions corresponding to outgoing or incoming waves in the contacts. , the retarded and the advanced444In practice these two solutions are usually obtained by adding an imaginary part to the energy. By taking the limit to zero of the imaginary part one of the two solutions is obtained. If the limit $`0^+`$ is taken the retarded solution is found, $`0^{}`$ gives the advanced. This can be seen from the Fourier transform of the time dependent Green’s function. solutions corresponding to outgoing and incoming waves in the contacts.

*Notation:* We will denote the retarded Green’s function with $`G`$ and the advanced with $`G^{}`$ (and maybe $`G^R`$ and $`G^A`$ occasionally). Here, CAPITAL $`G`$ denotes the full Green’s function and its sub-matrices $`G_1`$, $`G_d`$, $`G_{1d}`$ etc. Lowercase is used for the Green’s functions of the isolated subsystems, e.g., $`\left(EH_2\right)g_2=I`$.

Note that by using the retarded Green’s function of the isolated contact ($`g_2`$) in Eq. 9 we obtain the solution corresponding to a outgoing wave in the contact. Using the advanced Green’s function ($`g_2^{}`$) would give the solution corresponding to an incoming wave.

### 2.2 Self-Energy

The reason for calculating the Green’s function is that it is easier that solving the Schrödinger equation. Also, the Green’s function of the device ($`G_d`$) can be calculated separately without calculating the whole Green’s function ($`G`$). From the definition of the Green’s function we obtain:

$$\left(\begin{array}{ccc}EH_1& \tau _1& 0\\ \tau _1^{}& EH_d& \tau _2^{}\\ 0& \tau _2& EH_2\end{array}\right)\left(\begin{array}{ccc}G_1& G_{1d}& G_{12}\\ G_{d1}& G_d& G_{d2}\\ G_{21}& G_{2d}& G_2\end{array}\right)=\left(\begin{array}{ccc}I& 0& 0\\ 0& I& 0\\ 0& 0& I\end{array}\right)$$
(10)
Selecting the three equations in the second column:

$`\left(EH_1\right)G_{1d}\tau _1G_d`$ $`=`$ $`0`$ (11)
$`\tau _1^{}G_{1d}+\left(EH_d\right)G_d\tau _2^{}G_{2d}`$ $`=`$ $`I`$ (12)
$`\left(EH_2\right)G_{2d}\tau _2G_d`$ $`=`$ $`0`$ (13)
We can solve Eqs. 11 and 13 for $`G_{1d}`$ and $`G_{2d}`$:

$`G_{1d}`$ $`=`$ $`g_1\tau _1G_d`$ (14)
$`G_{2d}`$ $`=`$ $`g_2\tau _2G_d`$ (15)
substitution into Eq. 12 gives:

$$\tau _1^{}g_1\tau _1G_d+\left(EH_d\right)G_d\tau _2^{}g_2\tau _2G_d=I$$
(16)
from which $`G_d`$ is simple to find:

$$G_d=\left(EH_d\mathrm{\Sigma }_1\mathrm{\Sigma }_2\right)^1$$
(17)
where $`\mathrm{\Sigma }_1=\tau _1^{}g_1\tau _1`$ and $`\mathrm{\Sigma }_2=\tau _2^{}g_2\tau _2`$ are the so called self-energies.

Loosely one can say that the effect of the contacts on the device is to add the self-energies to the device Hamiltonian since when we calculate the Green’s function on the device we just calculate the Green’s function for the effective Hamiltonian $`H_{\text{effective}}=H_d+\mathrm{\Sigma }_1+\mathrm{\Sigma }_2`$. However, we should keep in mind that we can only do this when we calculate the Green’s function. The eigen-values and -vectors of this effective Hamiltonian are not quantities we can interpret easily.

For “normal” contacts, the surface Green’s functions $`g_1`$ and $`g_2`$ used to calculate the self-energies are usually calculated using the periodicity of the contacts, this method is described in detail in appendix B of Ref.  and in section 3 of Ref. .

### 2.3 The spectral function

Another important use of the Green’s function is the spectral function:

$$A=i\left(GG^{}\right)$$
(18)
which gives the DOS and *all* solutions to the Schrödinger equation!

To see this we first note that for any perturbation $`|v`$ we get two solutions ($`|\psi ^R`$ and $`|\psi ^A`$) to the perturbed Schrödinger equation:

$$(EH)|\psi =|v$$
(19)
from the advanced and retarded Green’s functions:

$`|\psi ^R`$ $`=`$ $`G|v`$ (20)
$`|\psi ^A`$ $`=`$ $`G^{}|v`$ (21)
The difference of these solutions ($`|\psi ^R|\psi ^A=(GG^{})|v=iA|v`$) is a solution to the Schrödinger equation:

$$(EH)(|\psi ^R|\psi ^A)=(EH)(GG^{})|v=(II)|v=0$$
(22)
which means that $`|\psi =A|v`$ is a solution to the Schrödinger equation for any vector $`|v`$!

To show that the spectral function actually gives *all* solutions to the Schrödinger equation is a little bit more complicated and we need the expansion of the Green’s function in the eigenbasis:

$$G=\frac{1}{E+i\delta H}=\underset{k}{}\frac{|kk|}{E+i\delta ϵ_k}$$
(23)
where the $`\delta `$ is the small imaginary part (see footnote 3), $`|k`$’s are all eigenvectors555Normalized! to $`H`$ with the corresponding eigenvalues $`ϵ_k`$. Expanding the spectral function in the eigenbasis gives:

$`A`$ $`=`$ $`i\left({\displaystyle \frac{1}{E+i\delta H}}{\displaystyle \frac{1}{Ei\delta H}}\right)`$ (24)
$`=`$ $`i{\displaystyle \underset{k}{}}|kk|\left({\displaystyle \frac{1}{E+i\delta ϵ_k}}{\displaystyle \frac{1}{Ei\delta ϵ_k}}\right)`$ (25)
$`=`$ $`{\displaystyle \underset{k}{}}|kk|{\displaystyle \frac{2\delta }{\left(Eϵ_k\right)^2+\delta ^2}}`$ (26)
where $`\delta `$ is our infinitesimal imaginary part of the energy. Letting $`\delta `$ go to zero gives:

$`A`$ $`=`$ $`2\pi {\displaystyle \underset{k}{}}\delta \left(Eϵ_k\right)|kk|`$ (27)
(here $`\delta (Eϵ_k)`$ is the delta function) which can be seen since $`\frac{2\delta }{\left(Eϵ_k\right)^2+\delta ^2}`$ goes to zero everywhere but at $`E=ϵ_k`$, integrating over $`E`$ (with a test function) gives the $`2\pi \delta \left(Eϵ_k\right)`$ factor. Eq. 27 shows that the spectral function gives us all solutions to the Schrödinger equation.

## 3 Response to an incoming wave

In the non-equilibrium case, reservoirs with different chemical potentials will inject electrons and occupy the states corresponding to incoming waves in the contacts. Therefore, we want to find the solutions corresponding to these incoming waves.

Consider contact $`1`$ isolated from the other contacts and the device. At a certain energy we have solutions corresponding to an incoming wave that is totally reflected at the end of the contact. We will denote these solutions with $`|\psi _{1,n}`$ where $`1`$ is the contact number and $`n`$ is a quantum number (we may have several modes in the contacts). We can find all these solutions from the spectral function $`a_1`$ of the isolated contact (as described above).

Connecting the contacts to the device we can calculate the wavefunction on the whole system caused by the incoming wave in contact $`1`$. To do this we note that a wavefunction should be of the form $`|\psi _{1,n}+|\psi ^R`$ where $`|\psi _{1,n}`$ is the totally reflected wave and $`|\psi ^R`$ is the retarded response of the whole system. Putting in the anzats $`|\psi _{1,n}+|\psi ^R`$ into the Schrödinger equation gives:

$`\left(\begin{array}{c}H_1+\tau _1+\\ H_d+\tau _1^{}+\tau _2^{}+\\ H_2+\tau _2\end{array}\right)\left(|\psi _{1,n}+|\psi ^R\right)`$ $`=`$ $`E\left(|\psi _{1,n}+|\psi ^R\right)`$ (31)
$`\left(\begin{array}{c}E|\psi _{1,n}+\\ \tau _1^{}|\psi _{1,n}+\\ 0\end{array}\right)+\left(\begin{array}{c}H_1+\tau _1+\\ H_d+\tau _1^{}+\tau _2^{}+\\ H_2+\tau _2\end{array}\right)|\psi ^R`$ $`=`$ $`E\left(|\psi _{1,n}+|\psi ^R\right)`$ (38)
$`H|\psi ^R`$ $`=`$ $`E|\psi ^R\tau _1^{}|\psi _{1,n}`$ (39)
(note the slight change in notation) and we see that $`|\psi ^R`$ is noting else the response of the whole system to a perturbation of $`\tau _1^{}|\psi _{1,n}`$, c.f., Eq. 9:

$$|\psi ^R=G\tau _1^{}|\psi _{1,n}$$
(40)
It is important to realize that the scattering states generated from Eq. 40, using all possible incoming waves from each contact, form a complete666Except for localized states in the device region. ON set of solutions to the full Schrödinger equation. Note that we have chosen the retarded response which means that the only part of the wave that is traveling towards the device is the incoming wave (part of $`|\psi _{1,n}`$). We will make full use of this fact below.

It will be useful to have the expressions for the device wavefunction $`|\psi _d`$ and contact wavefunction ($`|\psi _{1,2}`$). The device part is straightforward:

$$|\psi _d=G_d\tau _1^{}|\psi _{1,n}$$
(41)
and from Eq. 9 or Eq. 15:

$$|\psi _2=g_2\tau _2|\psi _d=g_2\tau _2G_d\tau _1^{}|\psi _{1,n}$$
(42)
Note, to calculate the wavefunction in the contact containing the incoming wave (contact $`1`$) we need to add the incoming wave, giving a slightly more complicated expression:

$$|\psi _1=\left(1+g_1\tau _1G_d\tau _1^{}\right)|\psi _{1,n}$$
(43)
Knowing the wavefunctions corresponding to incoming waves in different contacts enables us to fill up the different solutions according to the electron reservoirs filling the contacts.

## 4 Charge density matrix

In the non equilibrium case we are often interested in two quantities: the current and the charge density matrix. Lets start with the charge density (which allows us to use a self-consistent scheme to describe charging).

The charge density matrix is defined as:

$$\rho =\underset{k}{}f(k,\mu )|\psi _k\psi _k|$$
(44)
where the sum runs over all states with the occupation number $`f(E_k,\mu )`$ (pure density matrix) (note the similarity with the spectral function $`A`$, in equilibrium you find the density matrix from $`A`$ and not as described below). In our case, the occupation number is determined by the reservoirs filling the incoming waves in the contacts such that:

$$f(E_k,\mu _1)=\frac{1}{1+e^{(E_k\mu _1)/k_BT}}$$
(45)
is the Fermi-Dirac function with the chemical potential ($`\mu _1`$) and temperature ($`T`$) of the reservoir responsible for injecting the electrons into the contacts.

The wavefunction on the device given by an incoming wave in contact $`1`$ (see Eq. 41) is:

$$|\psi _{d,k}=G_d\tau _1^{}|\psi _{1,k}$$
(46)
Adding up all states from contact $`1`$ gives:

$`\rho _d[\text{contact }1]`$ $`=`$ $`{\displaystyle \underset{E=\mathrm{}}{\overset{\mathrm{}}{}}}\text{d}E{\displaystyle \underset{k}{}}f(E,\mu _1)\delta (EE_k)|\psi _{d,k}\psi _{d,k}|`$ (47)
$`=`$ $`{\displaystyle \underset{E=\mathrm{}}{\overset{\mathrm{}}{}}}\text{d}Ef(E,\mu _1){\displaystyle \underset{k}{}}\delta (EE_k)G_d\tau _1^{}|\psi _{1,k}\psi _{1,k}|\tau _1G_d^{}`$ (48)
$`=`$ $`{\displaystyle \underset{E=\mathrm{}}{\overset{\mathrm{}}{}}}\text{d}Ef(E,\mu _1)G_d\tau _1^{}\left({\displaystyle \underset{k}{}}\delta (EE_k)|\psi _{1,k}\psi _{1,k}|\right)\tau _1G_d^{}`$ (49)
$`=`$ $`\left[\text{Eq. }\text{27}\right]={\displaystyle \underset{E=\mathrm{}}{\overset{\mathrm{}}{}}}\text{d}Ef(E,\mu _1)G_d\tau _1^{}{\displaystyle \frac{a_1}{2\pi }}\tau _1G_d^{}`$ (50)
introducing the new quantity $`\mathrm{\Gamma }_1=\tau _1^{}a_1\tau _1=i\left(\mathrm{\Sigma }_1\mathrm{\Sigma }_1^{}\right)`$ we obtain the simple formula:

$$\rho [\text{from contact }1]=\frac{1}{2\pi }\underset{E=\mathrm{}}{\overset{\mathrm{}}{}}\text{d}Ef(E,\mu _1)G_d\mathrm{\Gamma }_1G_d^{}$$
(51)
The total charge density thus becomes a sum over all contacts:

$$\rho =\frac{2\text{(for spin)}}{2\pi }\underset{E=\mathrm{}}{\overset{\mathrm{}}{}}\text{d}E\underset{i}{}f(E,\mu _i)G_d\mathrm{\Gamma }_iG_d^{}$$
(52)
## 5 Probability Current

Having different chemical potentials in the reservoirs filling the contacts gives rise to a current. In the next section we will calculate this current in a similar way as the charge density was calculated. But to do this we need an expression for the current from the wavefunction.

In the continuum case we can calculate the current from the velocity operator. However, for a discrete Hamiltonian it is not so clear what the velocity operator is. Therefore, we derive an expression for the current from the continuity equation (using two contacts). In steady-state, the probability to find an electron on the device ($`\underset{i}{}|\psi _i|^2`$ where the sum runs over the device subspace) is conserved:

$`0`$ $`=`$ $`{\displaystyle \frac{\underset{i}{}|\psi _i|^2}{t}}={\displaystyle \underset{i}{}}{\displaystyle \frac{\psi |ii|\psi }{t}}={\displaystyle \underset{i}{}}\left({\displaystyle \frac{\psi |i}{t}}i|\psi +\psi |i{\displaystyle \frac{i|\psi }{t}}\right)`$ (53)
$`=`$ $`{\displaystyle \frac{i}{\mathrm{}}}{\displaystyle \underset{i}{}}\left(\psi |H|ii|\psi \psi |ii|H|\psi \right)={\displaystyle \frac{i}{\mathrm{}}}\left(\psi |H|\psi _d\psi _d|H|\psi \right)`$ (54)
$`=`$ $`{\displaystyle \frac{i}{\mathrm{}}}\left(\psi |H_d+\tau _1+\tau _2|\psi _d\psi _d|H_d+\tau _1^{}+\tau _2^{}|\psi \right)`$ (55)
$`=`$ $`{\displaystyle \frac{i}{\mathrm{}}}\left(\left[\psi _1|\tau _1|\psi _d\psi _d|\tau _1^{}|\psi _1\right]+\left[\psi _2|\tau _2|\psi _d\psi _d|\tau _2^{}|\psi _2\right]\right)`$ (56)
We interpret the term in the first (square) bracket as the incoming probability current into the device from contact $`1`$ and the second bracket from contact $`2`$. Generalizing to an arbitrary contact $`j`$ gives us the electric current (at one energy) as the charge (-e) times the probability current:

$$i_j=\frac{ie}{\mathrm{}}\left(\psi _j|\tau _j|\psi _d\psi _d|\tau _j^{}|\psi _j\right)$$
(57)
where $`i_j`$ is defined as positive for a current from the contacts into the device. We can now put in the expressions for the wavefunctions in the same way as for the density matrix.

## 6 Electrical Current

To calculate the total current through the device we only need to put in the wavefunction of the device and the contacts ($`|\psi _d`$, $`|\psi _1`$, $`|\psi _2`$) from Eqs. 41, 43 and 9 and add all the contributions together. Thus the current into the device from a incoming wave of one energy ($`E`$) in contact $`1`$ ($`|\psi _{1,n}`$) through the coupling defined by $`\tau _2`$ is:

$`i_{\text{2 from 1}}`$ $`=`$ $`{\displaystyle \frac{ie}{\mathrm{}}}\left(\psi _2|\tau _2|\psi _d\psi _d|\tau _2^{}|\psi _2\right)`$ (58)
$`=`$ $`{\displaystyle \frac{ie}{\mathrm{}}}(\psi _{1,n}|\tau _1G_d^{}\tau _2^{}g_2^{}\tau _2G_d\tau _1^{}|\psi _{1,n}\psi _{1,n}|\tau _1G_d^{}\tau _2^{}g_2\tau _2G_d\tau _1^{}|\psi _{1,n})`$ (59)
$`=`$ $`{\displaystyle \frac{ie}{\mathrm{}}}\psi _{1,n}|\tau _1G_d^{}\tau _2^{}(g_2^{}g_2)\tau _2G_d\tau _1^{}|\psi _{1,n})`$ (60)
$`=`$ $`{\displaystyle \frac{e}{\mathrm{}}}\psi _{1,n}|\tau _1G_d^{}\mathrm{\Gamma }_2G_d\tau _1^{}|\psi _{1,n}`$ (61)
Adding over the modes $`n`$ and noting that the levels are filled from the reservoir connected to contact $`1`$ gives (2 for spin):

$`I_{\text{2 from 1}}`$ $`=`$ $`2{\displaystyle \frac{e}{\mathrm{}}}{\displaystyle \underset{E=\mathrm{}}{\overset{\mathrm{}}{}}}\text{d}Ef(E,\mu _1){\displaystyle \underset{n}{}}\delta (EE_n)\psi _{1,n}|\tau _1G_d^{}\mathrm{\Gamma }_2G_d\tau _1^{}|\psi _{1,n}`$ (62)
$`=`$ $`{\displaystyle \frac{2e}{\mathrm{}}}{\displaystyle \underset{E=\mathrm{}}{\overset{\mathrm{}}{}}}\text{d}Ef(E,\mu _1){\displaystyle \underset{m,n}{}}\delta (EE_n)\psi _{1,n}|\tau _1|mm|G_d^{}\mathrm{\Gamma }_2G_d\tau _1^{}|\psi _{1,n}`$ (63)
$`=`$ $`{\displaystyle \frac{2e}{\mathrm{}}}{\displaystyle \underset{E=\mathrm{}}{\overset{\mathrm{}}{}}}f(E,\mu _1){\displaystyle \underset{m}{}}m|G_d^{}\mathrm{\Gamma }_2G_d\tau _1^{}\left({\displaystyle \underset{n}{}}\delta (EE_n)|\psi _{1,n}\psi _{1,n}|\right)\tau _1|m`$ (64)
$`=`$ $`{\displaystyle \frac{2e}{\mathrm{}}}{\displaystyle \underset{E=\mathrm{}}{\overset{\mathrm{}}{}}}\text{d}Ef(E,\mu _1){\displaystyle \underset{m}{}}m|G_d^{}\mathrm{\Gamma }_2G_d\tau _1^{}{\displaystyle \frac{a_1}{2\pi }}\tau _1|m`$ (65)
$`=`$ $`{\displaystyle \frac{e}{\pi \mathrm{}}}{\displaystyle \underset{E=\mathrm{}}{\overset{\mathrm{}}{}}}\text{d}Ef(E,\mu _1)\text{Tr}\left(G_d^{}\mathrm{\Gamma }_2G_d\mathrm{\Gamma }_1\right)`$ (66)
To get the total current through the device the current from contact two have to be subtracted away:

$$I=\frac{e}{\pi \mathrm{}}\underset{E=\mathrm{}}{\overset{\mathrm{}}{}}\text{d}E\left(f(E,\mu _1)f(E,\mu _2)\right)\text{Tr}\left(G_d^{}\mathrm{\Gamma }_2G_d\mathrm{\Gamma }_1\right)$$
(67)
which is exactly the Landauer formula for the current.

## Acknowledgments

The author is grateful for the inspiration from Supriyo Datta and Mads Brandbyge. The revised version was initiated from discussions with Casper Krag and Carsten Rostgård.


format: text

text:
# Measurement of refractive index and equation of state in dense He, H2, H2O, Ne under high pressure in a diamond-anvil cell

## I Introduction

There exist considerable evidences that pressure can dramatically change the electronic properties of a system. On the atomic scale, the variation of the electronic kinetic and potential energies with volume tends to make a contraction of the atomic size and a delocalisation of valence electrons. Consequently, spectacular transformations have been observed, such as rehybridization of atomic orbitals, like in carbon (sp2) to (sp3) across the graphite to diamond transition ; insulator-metal phase transitions, like in O2 Desgreniers90 , Xe Eggert89  or CsI Eremets98 ; the transformation of ice from an hydrogen bonded solid to an ionic solid in symmetric ice Goncharov96 . Most changes of the properties of a system under pressure are related to the microscopic electronic reorganization, and the determination of the refractive index over a large frequency domain is probably the most direct method to probe the electronic structure. The refractive index expresses the response of the distribution of charges to a perturbing electric field. In order to extract a detailed change of the charge distribution of a system with pressure through the measurement of the refractive index, a large frequency domain has to be investigated. This is unfortunately very difficult at high pressure. However, already interesting electronic evolutions have been obtained from measurements of the refractive index performed in the visible, for instance in hydrogen with the red shift of an electronic transition that parallel the closure of the electronic gap (see references below).

Also, an accurate determination of the refractive index under pressure has a clear metrological importance. In Brillouin spectroscopy, where the refractive index of the sample and of the pressure transmitting medium have to be known to extract the sound velocity from the frequency shift Grimsditch88 . In the newly developed laser shock measurements on pre-compressed targets to extract the shock front velocity from the phase shift of the VISAR (velocity interferometric system for any reflector) diagnostics Lee02 . Finally, to measure easily the thickness of the sample chamber in the diamond-anvil cell (DAC) with a single spectrometer, when the sample or the pressure transmitting medium is transparent. This can be directly used to measure the Equation of State (EoS) of the sample (see below), or help to control the experimental conditions, for instance by evaluating the strain of the studied sample.

In the DAC, the tips of the diamond anvils that enclose the sample form a Fabry-Perot interferometer. This geometry have already been exploited to measure the refractive index of transparent samples at high pressure. Most of the experimental methods developed so far have been based on the measurement of fringes contrast (VanStraaten82 , VanStraaten88 , Hemley90b , Evans98 ). Other were based on a controlled sample thickness (Eggert92 Balzaretti93 ). In one work LeToullec89 , the combination of the measurement of parallel white beam interference spectrum and convergent monochromatic beam inferference rings allowed the determination of both refractive index $`n`$ and thickness $`t`$ of the sample. We present here an improvement of this method.

With this setup, we have also been able to measure the Equation of State of helium, by estimating the sample chamber volume changes upon compression. Such a method has been previously described for H2 VanStraaten88 , Evans98 . We show here that with numerical photograph technology and the accurate determination of the refractive index, its accuracy becomes comparable to other techniques. Therefore, it becomes an interesting alternative to X-Ray methods (X-Ray diffraction and X-Ray absorption Shen02 ) for lights elements, particularly light liquids, or to adiabatic compressibility measurement of Brillouin spectroscopy Grimsditch96 .

Three goals have mainly motivated the present accurate refractive index measurements on four systems, He, H2, Ne and H2O. First, H, He and H2O are the main constituents of planetary interiors and the determination of their properties up to planetary core conditions is motivating a lot of effort. Second, He and Ne are the most commonly used quasi-hydrostatic pressure transmitting media in DAC. In both cases, the determination of the refractive index should be useful to analyze some spectroscopic measurements. Finally, they have simple electronic structures with different electronic organisations and the present data set could be interesting to test the confidence of ab-initio calculations in predicting the optical properties.

The experimental method and a discussion of the uncertainties is presented in section II. The refractive index measurements are presented in section III. The analysis of the data in terms of the changes in the polarizability with density is presented in section IV. Applications of the method to the determination of the equation of state and to the study of phase transitions are described in section V.

## II Experimental procedure

### II.1 High pressure techniques

In these experiments we used a membrane diamond anvil cell described elsewhere (Letoullec88 ) with a large optical aperture ($`\pm `$ 30) and equipped with 400 or 500 $`\mu `$m table diamonds. Samples have been loaded in a rhenium gasket at room temperature, by the use of high pressure loading techniques. In the case of H2O, ultra-pure water was inserted in the gasket hole using a syringe. As explained below, thick samples are needed for a better accuracy of the refractive index measurements. For this reason, double rhenium gaskets have been used in the highest pressure runs. These gaskets are prepared in the following manner: two rhenium gaskets are pre-indented, then polished on one side and glued together. With this method, the thickness of an helium sample was of the order of 50 $`\mu `$m at 200 kbar, 30 % more than with a simple gasket. The pressure gauge was ruby, calibrated with the quasi-hydrostatic pressure scale Mao86 , that allows pressure measurement with an accuracy of the order of 0.5 kbar. The absolute accuracy is the one of the ruby scale, i.e at maximum 2% over the pressure range investigated here.

### II.2 White light interferences

The refractive index setup is modified from LeToullec et al.’s one (see Fig. 1). The white light interference pattern in transmission is recorded over the range 500nm to 700nm with a double substractive spectrometer (Dilor). A small area, varying between 20 x 20 and 10 x 10 $`\mu `$m2, was selected to perform the measurement. Under high pressures ($`P>250`$ kbar), the contrast drastically decreased out of a central flat region, because of the strain of the diamonds under pressure. The wavelength difference between two intensity maxima :

$$\mathrm{\Delta }\lambda =\frac{\lambda ^2}{2n(\lambda )t},$$
(1)
where $`t`$ is the thickness of the sample, $`n(\lambda )`$ its refractive index, was determined. In that wavelength range, dispersion could be observed only in H2O samples; in that case, a linear dependency of $`1/n`$ on $`\lambda `$ was assumed in order to fit Eq. 1.

### II.3 Monochromatic interferences

The setup presented in Fig. 1 allows the measurement of the reflective Fabry-P rot rings produced by the cavity between the diamonds, when illuminated by a convergent monochromatic beam. We used an Argon laser ($`\lambda =`$ 488.0 nm) for helium samples 1 to 3, and a He-Ne laser ($`\lambda =`$ 632.8 nm) for H2O, H2, Ne and He samples.

The interference order $`k`$ at angle $`i`$ is defined by:

$$k=\frac{1}{\lambda }2t(n^2\mathrm{sin}^2i)^{1/2};$$
(2)
if $`k`$ is an integer, a minimum in intensity is observed.

At the center of the pattern, the interference order $`k_0=2n(\lambda )t/\lambda ,`$ can be deduced from white light interferences measurement. The angular position $`i_p`$ of the $`p^{th}`$ intensity minimum has an interference order $`k_p=[k_0]p+1`$. Refractive index and thickness can be deduced from two fringes ($`p`$ and $`p^{}`$) angular positions:

$$n^2=\frac{k_p^{}^2\mathrm{sin}^2i_pk_p^2\mathrm{sin}^2i_p^{}}{k_p^{}^2k_p^2}$$
(3)
and

$$t^2=\frac{\lambda ^2}{4}\frac{k_p^{}^2k_p^2}{\mathrm{sin}^2i_p\mathrm{sin}^2i_p^{}}$$
(4)
Interference rings were projected and recorded by a CCD camera (C1, see Fig. 1). The radius $`R_p`$ of the rings on C1 is thus measured. Because of optical aberrations, mainly caused by the objective lenses L2, the relationship between $`R_p`$ and $`i_p`$ was non trivial. To correct these abberations, a reference spectrum was recorded. This spectrum consisted in fringes formed by an empty indented and drilled gasket, placed in the DAC. The adjustment of the optical setup was performed as it is described herebelow. First, the cell was centered and perpendicularly oriented on the laser beam. Then M1, C1, L3, L2, L1 and D1 were positioned and the calibration spectrum recorded. All theses elements were mounted on ($`x`$, $`y`$, $`z`$, $`\theta `$) linear and rotation stages. After calibration, all the elements, except the cell, remained fixed. In the first experiments, the cell was regularly removed to record the white light interferences pattern. The best position of the cell (in the $`x`$ direction) was determined, within $`\pm `$ 10 $`\mu `$m, by optimization of the rings pattern stability with respect to Y and X translations. This allowed a reproducibility of the index measurement of $`\mathrm{\Delta }n=\pm 5.10^4`$. Each digital Fabry-Perot pattern was automatically processed by a program that evaluates thickness and index for all available pair of rings. This procedure allows the analysis of numerous spectra, and again, optimization of the $`x`$ position of the sample. Tentative error bars have been estimated for each measurement, taking into account two factors: first, the uncertainty on rings radius $`R_p`$; second, the scatter of the refractive index given by different couples of rings (Eq. 3).

The precision of refractive index measurements is limited by the following effects:

\- the precision of $`k_0`$ measurement by white light interferences. The relative precision on $`k_0`$ estimate is better than 10-3 ($`\mathrm{\Delta }k<0.2`$) for weakly dispersive media (He, Ne and H2) and of approximately 2.10-3 for dispersive media (H2O). This measurement must also be carried out in the same sample region than the region illuminated by the laser beam for monochromatic interferences, in order to prevent any error caused by non-parallelism of the diamonds.

\- the number of Fabry-P rot rings. They must be as numerous as possible to statistically increase the precision of the measurement. In these experiments, a minimum of four rings has been observed. For this purpose, the angular aperture of the optical setup was increased to 25. The number of rings is also proportional to the sample thickness; therefore, thick samples have been used.

\- the elastic strain (cupping) of the diamonds. This is the main source of uncertainty of our measurements, that limits the pressure range of this method. This cupping can be first observed by translating the cell in $`y`$ or $`z`$ directions and looking at the modification of the rings pattern. It can be noted that these measurements constitute an accurate mean of studying the strain of the anvils. At the beginning of the compression, inverse cupping, caused by the gasket strength, is observed. Upon pressure increase, the diamonds become flat and then cupping occurs. It has been observed approximately at the same thickness of a simple rhenium gasket (40 to 45 $`\mu `$m for H2O samples), and for pressures varying between 100 and 200 kbar. When pressure increases, it perturbs the rings pattern out of a central flat zone, up to a point that makes measurements impossible. To reduce this intrinsic limitation, thick gaskets can be used because they increase the pressure at which cupping appears. This elastic strain also makes the focusing diagnostic less precise.

\- the pressure gradients in the solid phase. The ruby chip can not be directly located where the sample was illuminated, because it perturbs the Fabry-P rot interferometric rings. It was located at $``$ 60 $`\mu `$m from the refractive index measurement point. In separate runs we have estimated the radial pressure gradient in ice VII, in the same sample geometry; this gradient is negligible up to 200 kbar and reaches 0.33 kbar/$`\mu `$m at 350 kbar. In ice, the pressure was corrected by this factor. Actually, this pressure gradient is also an intrinsic limitation of the method, because the extend of the source of interference rings is 100 $`\mu `$m, if $`i<25^{}`$; at high pressure, rings with $`i<15^{}`$ are analyzed in order to reduce this error.

### II.4 EoS measurement

An analysis of numerical images of the sample taken in transmitting light allows the measurement of the surface $`S`$ of the sample, with a precision of 1 %. The volume $`V_S`$ is then obtained by multiplying $`S`$ with the thickness of the sample, obtained at the same time as its refractive index. The uncertainty on volume determination is mainly caused by the uncertainty on $`S`$. No corrections taking into account the cupping of the diamond was necessary for the measurements presented in Sec. V, but these corrections are in principle possible. We made the assumption that the amount of matter trapped in the sample chamber remained constant upon pressure increase. In that case, the EoS $`V(P)`$ can be simply deduced from $`V_S(P)`$ by choosing a reference point $`V_0(P_0)`$. For a given run, the molar volume of the sample can be expressed as:

$$V_{\mathrm{mol}}(P)=V(P)\frac{V_{\mathrm{mol}}(P_{\mathrm{ref}})}{V(P_{\mathrm{ref}})},$$

$`V(P)`$ being the volume of the sample at the pressure $`P`$ and $`V(P_{\mathrm{ref}})`$ the volume of the sample at the reference pressure $`P_{\mathrm{ref}}`$. We noticed that the precision of the volume measurement increased when the thickness of the sample decreases. A similar observation has already been made by Evans and Silvera, 1998, Evans98 . This is probably because the geometrical imperfections of the gasket hole, extensively described VanStraaten88 , are smaller for thin samples.

## III Refractive index data

### III.1 Helium

The refractive index data obtained in 4 runs up to 210 kbar are shown in Fig. 2. During the same experimental run, the orders of interference of the rings of the fringe pattern have been followed between each pressure steps. The good agreement between the data of the different runs data indicates that the uncertainty associated to the white light determination of the order of interference is negligible. For all data shown in Fig. 2, the error bars are smaller than the size of the symbols. The difference between the fitted index law and LeToullec measurements is smaller than 2.10-3, except for the highest pressure points where LeToullec measurements were less precise.

### III.2 Neon

One run has been performed with a neon sample, up to 27 GPa. The evolution of the refractive index of neon with pressure is plotted in Fig. 3, The refractive index of Ne had been determined previously to a maximum pressure of 1 GPa from capacitance measurements Lallemand77 . The present data are in very good agreement with this previous determination, in the pressure range of overlap.

### III.3 Hydrogen

The evolution of the refractive index of hydrogen with pressure is plotted in Fig. 4.

There is a systematic difference between our data and previous determinations, either the one by brillouin scattering Shimizu81  (1.5 %) or the one by white light interference Evans98  (around 1 %, which is within their estimated error bars). The matching of refractive index of ruby and hydrogen has been observed (annulation of the contrast) at 35 $`\pm 10`$ GPa. This provides an anchorage point to the refractive index curve. At this pressure, the refractive index of ruby is expected to be 1.77 Balzaretti93 ; this value is in good agreement with the present measurements and that gives a good confidence in our highest pressure runs.

### III.4 Water and ice

For H2O, 5 experimental runs have been made, up to 35 GPa. Due to the dispersion of the refractive index of H2O, of the order of 1 % over the 500 - 800 nm range, the determination of the order of interference from white light measurements is less precise and consequently small deviations between different experimental runs can be detected in the fluid phase. However, our measurements agree very well with precise interferometric measurements in the fluid phase Vedam78 . At higher pressures, our data are within the scatter of a previous Brillouin scattering determination Shimizu96 . At pressures higher than 20 GPa, the strength of ice causes pressure gradients in the sample, that can perturb the pressure estimation. To correct it, the pressure measured by the ruby chip $`P_R`$ was then corrected by a factor $`\mathrm{\Delta }P=d\times \mathrm{d}P/\mathrm{d}r`$, $`d`$ being the distance between the ruby chip and the center of the sample, and the last factor has been estimated in a separate run by mapping the sample chamber with ruby (in kbar/$`\mu `$m): $`\mathrm{d}P/\mathrm{d}r=(P_R\times 0.425)/t.`$ The estimated error bars plotted in Fig. 5 increase drastically as pressure increases, showing the decrease of the quality of the fringe pattern associated to the strain in the diamonds. Refractive index matching between ice VII and ruby has been observed at $`P=40.5\pm 10`$ GPa. At this pressure, the refractive index of ruby is approximately 1.765 Balzaretti93 , 1.2 % lower than the extrapolation of our measurements. That is within the uncertainty of our highest pressure measurements.

The functions obtained by a fitting procedure to our experimental data are gathered in table 1. The use of power laws is justified, to a certain extent, by the fact that the Lorentz-Lorenz factor defined in Eq. 5 is approximately constant (see below). We used as least parameters as possible to fit the data. The quality of the various fits is shown in the various figures that present the refractive index of He, H2, Ne, H2O versus pressure.

## IV Discussion

### IV.1 Polarizability

The refractive index is generally related to the polarizability of the molecular entities of the system through the the Lorentz-Lorenz formulation:

$$LL=\frac{1}{\rho }\frac{ϵ1}{ϵ+2}=\frac{4\pi N_A}{3}\alpha .$$
(5)
In Eq. 5, $`N_A`$ is the Avogadro number, $`ϵ`$ the relative dielectric constant, and $`\alpha `$ the molecular polarizability of the compound. In the static case, this equation is called the Clausius-Mossoti relation. It arises from the expression of the local electrical field for dense materials: $`𝐄_{\mathrm{𝐥𝐨𝐜}}(𝐫)=𝐄_{\mathrm{𝐦𝐚𝐜𝐫𝐨}}(𝐫)+4\pi 𝐏(𝐫)/3+𝐄^{}(𝐫)`$. $`𝐏(𝐫)`$ is the polarization and $`𝐄^{}(𝐫)`$ is the field created by the close neighbors of the electronic cloud centered on $`𝐫`$. Eq. 5, corresponds to $`𝐄^{}(𝐫)=0`$. This is the case on symmetry center of a a cubic crystal. For hydrogen, helium and neon, Eq. 5 is commonly used Evans98 LeToullec89  because they are constituted by spherical molecular entities (in the case of hydrogen the molecules are freely rotating), organized in highly symmetric structures. Also, for these systems in the pressure range covered here, the electronic transitions are far from the visible. The imaginary part of the refractive index should then be very small and $`ϵ=n^2`$ is a very good approximation.

In the case of water, the polarizability is generally written as the sum of three components: the polarizability due to the orientation of the permanent dipole of the water molecule, the polarizability due to the electronic changes associated to the inter and intra- molecular modes and the polarizability associated to the electronic transitions within the water molecules. The contribution of the permanent dipoles operates at frequencies below 105 Hz. The contribution due to the intermolecular modes is important below 1011 Hz in the liquid and 105 Hz in the solid Petrenko99 . The contributions due to the intra-molecular modes should be important below 1014 Hz (the highest frequency vibration mode of H2O is a stretching mode at 3656 cm-1 in the gas phase). Above 1016 Hz electronic transitions should be considered. The optical window covered here is just in between these dispersion phenomena, which implies a small dispersion of the refractive index and a small absorption. We did not observe any pressure effect on the dispersion. In Fig. 6, the Lorentz-Lorenz factor of water and ice calculated using the present data and published equations of state of liquid water, Ice VI and Ice VII (Saul89 , Bridgman42 , Hemley87 , Loubeyre99 ) is plotted, as a function of molar density $`\rho `$. It is interesting to note that the variation of $`LL`$ plotted on Fig. 6 does not depend on the microscopic arrangement of the molecules but only of the density. This tends to prove that the hypothesis $`𝐄^{}(𝐫)=0`$ in the local field is also reasonably valid in water. Consequently, Eq. 5 should be valid here to analyze the data. It is seen in figure 6 that the $`LL`$ factor decreases by approximately 20 % up to the maximum pressure of investigation. The approximation $`LL`$ constant, that is used to analyse spectroscopic data of water under pressure, such as the Brillouin scattering sound velocity measurements Grimsditch96 , is not correct because it leads to an underestimation of the refractive index of the order of 3 % at 35 GPa.

It has been proposed from experimental data below 1 GPa, that the evolution of the Lorentz-Lorenz factor of various molecular systems should follow an universal curve when plotted in reduced units Lallemand77 , after small statistical dipole-induced dipole corrections. In our case, their relative effect on $`LL`$ is smaller than 10-3 and can be neglected. In figure 7, the change of the Lorentz-Lorenz factors of He, Ne, H2, H20 vs. density is plotted in reduced units, i.e. as $`LL/LL_{gas}`$ (because $`LL_{gas}`$, taken in Achtermann93 Lallemand77 Schiebener90 , contains the isolated atom polarizability) vs. $`N_A\rho \sigma ^3`$, (where $`\sigma `$ is the Lennard-Jones parameter of the considered compound). The following Lennard-Jones parameters have been used: He, $`\sigma =2.556`$ Å Klein76 ; Ne, $`\sigma =2.83`$ Å Vos91 ; H2, $`\sigma =3.06`$ Å Ross83 ; H2O, $`\sigma =2.85`$ Å Belonoshko91 . Calculations have been performed using the following published EoS: He LeToullec89 Loubeyre93 , Ne Vos91 Lallemand77 , Finger81 Kortbeek88 , and H2Mills77 Hemley90 Pratesi97 . The scatter of the points give an idea of the experimental errors. The large scatter at low pressure observed for He and H2, the most compressible materials, is mainly caused by the uncertainty on density calculation; moreover, uncertainty on $`LL`$ is larger when $`n`$ is close to one. At higher compression, relative uncertainty decreases.

In Fig. 7, a common trend is observed, namely a decrease of the polarizability with density for the four systems. In fact, the polarizability is related to the extent of the electronic cloud. In a dense environment, the atoms and molecules themselves compress to lower their interactions with their neighbors, and a decrease of the polarizability is expected. Over the compression range covered here, the scaled $`LL`$ factors of H2 and of He, two isoelectronic species, exhibit the same density evolution. However, in that pressure range, the fundamental difference between the density effect on electronic properties of these two media would be extracted from dispersion data. In He, it has been shown that the main contribution of the change of polarizability with density is taken into account by the blue shift of the 1S-2P transition energy LeToullec89 . In contrast, the change of polarizability in hydrogen has been related to the red shift of the excitonic level $`X^1\mathrm{\Sigma }_g^+`$ to $`B^1\mathrm{\Sigma }_u^+`$ that should lead to the metal hydrogen VanStraaten88b Garcia92 . But different evolutions are observed in Ne and H2O; the proposed universal behavior Lallemand77  is nor valid nor theoretically justified for a large perturbation by density. In the past 20 years, numerous computations and measurements have been performed for rare gas fluids (see Dacre81 Achtermann93 Hattig99  and references therein), in order to reproduce, with success for rare-gas fluids, the two and three order refractivity virial coefficients. Yet, they have not been extended in the dense fluid or the solid phases. We are not aware of any calculation on the dielectric properties of water under pressure. On the other hand, the calculation of the dielectric properties of hydrogen is a subject of great current interest (see Souza00  and references therein). We present here an interesting database of the optical dielectric properties of low Z compounds under pressure. We hope that it could be helpful to validate ab-initio approaches.

### IV.2 Gladstone-Dale relation

The Gladstone-Dale relation expresses as Setchell02 :

$$\frac{\mathrm{d}n}{n1}=\frac{\mathrm{d}\rho }{\rho },$$
(6)
which leads to $`n=1+a\rho `$. It can be derived by differentiating Eq. 5 and expanding it as a function of $`u=n1`$, if the polarizability is assumed to remain constant :

$$\frac{\mathrm{d}u}{u}\left(1\frac{1}{6}u\frac{5}{6}u^2+\mathrm{}\right)=\frac{\mathrm{d}\rho }{\rho }.$$

Thus, this relation corresponds to the zero-order truncated series expansion of $`\mathrm{d}u`$ vs. $`\mathrm{d}\rho `$. For $`n<1.8`$, the first order term remains smaller than 13 % of the zero-order term. Moreover, the volume dependency of $`\alpha `$ will tend to compensate the effect of the first and second order terms.

The Gladstone-Dale relation has been widely used to analyze the shock-wave VISAR data. In fact, a transparent window is often attached to the rear face of the sample, together whith an optically reflecting interface, to prevent the formation of a reflected rarefaction wave. The high pressure refractive index of this window is needed to determine the rear face velocity of the sample by VISAR technique Setchell02 . It can also been useful if the sample itself is transparent, and if the shock front becomes reflective. This is the case for pre-compressed samples Lee02 . An advantage of the Gladstone-Dale relation is that refractive index expresses as a function of the density of the material alone (intrinsic temperature effect is not taken into account), and very simply. In this peculiar form, it allows to express straightforwardly the actual velocity of the reflecting interface as a function of its apparent velocity, for a planar single shock wave.

For this reason, we think that a presentation of our results in a form of Glastone-Dale fit could be useful for the high pressure community. They are presented in Fig. 8. The data are close to straight lines, even if a small curvature is obvious in all curves. The Gladstone-Dale relation is often presented under the following form: $`n=a+b\rho `$ Setchell02 . We used this function to fit our data. The results are presented on Table 2.

## V Applications

### V.1 Equation of State

The EoS of helium in the fluid and solid phase has been measured up to 13 GPa, as presented in Fig. 9, during runs 1 to 3. The volume of the sample was obtained by the multiplication of its thickness by its surface. Sample thickness was varying between 37 $`\mu `$m and 23 $`\mu `$m (simple gaskets).

The evolution of the sample thickness during run 4 (double rhenium gasket) is also plotted in Fig. 10. It is interesting to note that the thickness of the sample decreased non-linearly upon pressure increase but remained approximately constant upon pressure decrease: from 47 $`\mu `$m to 48.5 $`\mu `$m. This can be used experimentally, when the sample thickness knowledge is needed Shen02 : one can perform measurements upon pressure decrease.

Our EoS is compared with Brillouin spectroscopic measurements LeToullec89  in the fluid phase and x-ray data in the solid phase Mao88 Loubeyre93 . At moderate pressures ($`P<3`$ GPa), it seems that the compressibility of the sample is overestimated; this effect can be attributed to either leakage of the sample or geometrical effects. Between 3 and 14 GPa, agreement with other data is excellent. No evident leakage has been observed in that pressure range. This shows that this method for EoS measurements is very interesting for liquids, where x-ray diffraction measurements are still difficult to perform and where there is no radial pressure gradient in the sample. In particular, these measurements could be carried out for high temperature and high pressure liquids.

### V.2 Detection of phase changes

The volume discontinuities at each phase changes can be deduced from measurement. In fact, we observed no evidence the the $`LL`$ factor was discontinuous at any phase change studied. Thus, we assumed that this function depends on density alone. By differentiating Eq. 5, $`\mathrm{\Delta }\rho /\rho `$ can then be deduced from the observed $`\mathrm{\Delta }n`$. The results are summarized table 3. Our data compare very well with those obtained by other techniques. Thus, careful refractive index measurement can detect phase changes with volume changes as small as 2 %.

## VI Conclusion

We have performed accurate optical refractive index measurements of H2O, He, H2 and Ne up to 35 GPa and at ambient temperature in a diamond anvil cell. These data should b useful to analyze some spectroscopic measurements on these systems at high pressure. The change of the refractive index has been expressed in terms of the change of the polarizability of the molecular entities in their dense surroundings. It is still challenging to compute such electronic properties. The present set of data could be useful in validating the various theoretical approaches. Our method also allows the direct measurements of the Equations of State of transparent media: this has been carried out for fluid and solid He, with less than 1 % error. This method should then be interesting to determine the equation of state of dense fluids or dense amorphous media.

###### Acknowledgements.

We thank M. Berhanu for his help during Neon experiments.


text:
# Cauchy horizon stability in self-similar collapse: scalar radiation.

## I Introduction

Perhaps the richest source of examples of space-times admitting naked singularities is the class of spherically symmetric self-similar space-times. There is an extensive literature on the topic; the recent review of self-similarity in general relativity by Carr and Coley carr (1) provides a suitable bibliography. Of particular note in this class are the perfect fluid solutions studied by Ori and Piran OP (2), the massless scalar field solutions studied by Christodoulou christo1 (3) and by Brady brady (4) and the $`SU(2)`$ sigma model solutions studied by Bizon and Wasserman bizon (5). We mention these because (i) the matter model has particular interest, for either physical or mathematical reasons and (ii) these self-similar solutions are of interest in studies of critical phenomenon gundlach (6). More generally, self-similar solutions admitting naked singularities are of interest because of what they may tell us about cosmic censorship. Intriguingly, the evidence is not all in one direction. Recent work has indicated the stability of perfect fluids admitting naked singularities in the class of perfect fluid space-times harada (7), while for the case of the massless scalar field, generic spherical perturbations of self-similar initial data which correspond to naked singularities will lead to censored singularities christo2 (8). Also, within the class of self-similar spherically symmetric space-times, the sectors corresponding to censored and to naked singularities are both topologically stable nolan (9).

With these results in mind, the aim here is to begin a comprehensive study of the stability of Cauchy horizons in self-similar collapse. In the case of charged rotating black holes, the instability of the Cauchy (or inner) horizon has been firmly established (see brady2 (10) for a review). This instability is in one way easily understood; an observer crossing the inner horizon views the entire history of the external universe in a finite amount of proper time, and so time-dependent perturbations of the exterior suffer an infinite blue-shift on crossing the horizon. This instability mechanism which can be “read off” the conformal diagram, does not have a counterpart in self-similar collapse which leads to globally naked singularities (see Figures 1 and 2). At best, one can speculate that the curvature at the regular center which diverges in the limit as the scaling origin is approached makes itself felt by perturbations approaching the Cauchy horizon. This is by no means convincing, and so a rigorous analysis is required. We begin this analysis here by studying the propagation of scalar radiation in a fixed background (spherically symmetric, self-similar) space-time which admits a Cauchy horizon.

In the following section, we define the class of space-times of interest and obtain some useful relations for the metric functions thereof. We consider spherically symmetric space-times admitting a homothetic Killing vector field whose energy-momentum tensor obeys the dominant energy condition. (A complete account of energy conditions in spherical symmetry is given in the appendix.) For generality, no further restrictions are imposed at this stage, although some differentiability conditions at the past null cone of the scaling origin and at the Cauchy horizon will be imposed. Using co-ordinates adapted to the homothety and to the past null cones of the central world-line, simple conditions can be given on the metric which determine the visibility or otherwise of the singularity at the scaling origin $`𝒪`$. This allows a simple way of identifying both the past null cone $`𝒩`$ of $`𝒪`$ and the Cauchy horizon $``$. In Section 3, we determine the behaviour of completely general time-like geodesics (i) crossing $`𝒩`$ and (ii) crossing the Cauchy horizon. These are used to calculate fluxes of the scalar field at the respective surfaces. The minimally coupled scalar wave equation is studied in the next section. A mode decomposition relying on the Mellin transform is used, and the asymptotic behaviour of the general solution at $`𝒩`$ is determined. This is used to impose the boundary condition that an arbitrary observer with unit time-like tangent $`v^a`$ measures a finite flux $`v^a_a\mathrm{\Phi }|_𝒩`$. We also demand that the influx at $`𝒥^{}`$ be finite. The modes not ruled out by these boundary conditions are then allowed to evolve up to the Cauchy horizon and the flux $`v^a_a\mathrm{\Phi }|_{}`$ is calculated. Our principal result is that this flux is finite for all the cases we consider.

## II Self-similar spherically symmetric space-times admitting a naked singularity.

We will consider the class of space-times which have the following properties. Space-time $`(,g)`$ is spherically symmetric and admits a homothetic Killing vector field. These symmetries pick out a scaling origin $`𝒪`$ on the central world-line $`r=0`$ (which we will refer to as the axis), where $`r`$ is the radius function of the space-time. We assume regularity of the axis to the past of $`𝒪`$ and of the past null cone $`𝒩`$ of $`𝒪`$. We will use advanced Bondi co-ordinates $`(v,r)`$ where $`v`$ labels the past null cones of $`r=0`$ and is taken to increase into the future. Translation freedom in $`v`$ allows us to situate the scaling origin at $`(v=0,r=0)`$ and identifies $`v=0`$ with $`𝒩`$. The homothetic Killing field is

$$\stackrel{}{\xi }=v\frac{}{v}+r\frac{}{r}.$$

The line element may be written

$$ds^2=2Fe^{2\psi }dv^2+2e^\psi dvdr+r^2d\mathrm{\Omega }^2,$$
(1)
where $`d\mathrm{\Omega }^2`$ is the line element of the unit 2-sphere. The homothetic symmetry implies that $`F(v,r)=F(x),\psi (v,r)=\psi (x)`$ where $`x=v/r`$. The only co-ordinate freedom remaining in (1) is $`vV(v)`$; this is removed by taking $`v`$ to measure proper time along the regular center $`r=0`$.

We will not specify the energy-momentum tensor of $`(,g)`$, but will demand that it satisfies the dominant energy condition. A complete description of energy conditions in spherical symmetry is given in Appendix A. Of these, the following will be used. (These are equations (41), (42) and (46) respectively.)

$`x\psi ^{}0,`$ (2)
$`e^\psi (F^{}+xF^2e^\psi \psi ^{})0,`$ (3)
$`12F+2x(F^{}+F\psi ^{})0.`$ (4)
We impose the following regularity conditions at the axis. As previously mentioned, we take $`v`$ to be proper time along the axis for $`v<0`$. Noting that $`x\mathrm{}`$ on this portion of the axis, (1) then gives

$$\underset{x\mathrm{}}{lim}2Fe^{2\psi }=1.$$
(5)
The other regularity condition that we use is that all curvature invariants are finite on $`r=0,v<0`$. In the present case, the (invariant) Misner-Sharp mass is given by

$$E=\frac{r}{2}(12F).$$

Then $`E/r^3`$ is a curvature invariant; this term has the same units as e.g. $`R`$ and $`\mathrm{\Psi }_2`$. Demanding that $`E/r^3`$ be finite on the axis yields

$$\underset{x\mathrm{}}{lim}F=\frac{1}{2}.$$
(6)
Combining (5) and (6) gives these regularity conditions:

$$F(\mathrm{})=\frac{1}{2},\psi (\mathrm{})=0.$$
(7)
We define the interior region $`_{int}`$ of space-time to be the interior of $`𝒩`$, i.e. the interior of the causal past of $`𝒪`$. The exterior region $`_{ext}`$ is defined to by $`=_{int}𝒩_{ext}`$. (These definitions are in line with those of christo1 (3).) We assume that the metric is regular throughout $`_{int}𝒩`$ \- this set does not include $`𝒪`$ \- by which we mean $`F,\psi C^2(\mathrm{},0]`$. As a Cauchy horizon can only form in $`_{ext}`$, we assume further that $`F,\psi C^2(\mathrm{},x_{})`$ for some $`x_{}>0`$. As we will see, if a Cauchy horizon develops, it must be of the form $`x=x_c`$ for some $`x_c>0`$. Our assumption is that the metric is regular at least up until the Cauchy horizon.

Since we are studying collapse, our assumptions must include some statement of regularity - in the sense of the absence of trapped surfaces - of an initial configuration. The 2-sphere $`(v,r)`$ is trapped if and only if

$$\chi (v,r):=g^{ab}_ar_br<0.$$

In the present case, this is equivalent to $`F<0`$, and implies that the condition for an apparent horizon is $`F=0`$. So in order to express the notion that the matter is initially in some non-extreme state, we rule out trapped or marginally trapped surfaces in the interior region $`_{int}`$. We will also demand that $`𝒩`$ is not foliated by marginally trapped surfaces, and so we take

$$F(x)>0\text{for all}x(\mathrm{},0].$$

Next, we point out the inevitability of there being a curvature singularity at $`𝒪`$. Any curvature invariant which has units $`L^2`$ is of the form $`c(x)r^2`$. For example,

$$\frac{E}{r^3}=\frac{12F}{2r^2}.$$

This term diverges as we approach $`𝒪`$ along the null line $`x=0`$ unless $`F(0)=\frac{1}{2}`$. But subject to the assumption that $`F>0`$ for $`x<0`$, we see that the surfaces $`x=x_c<0`$ are time-like. So we may also approach $`𝒪`$ along $`x=x_c<0`$, and we then see that $`E/r^3`$ diverges unless $`F\frac{1}{2}`$ on $`(\mathrm{},0]`$. Applying the same reasoning to the invariant

$$\frac{E}{r^3}+\mathrm{\Psi }_2+\frac{R}{12}=\frac{1}{2r^2}(12F+2x(F^{}+F\mathrm{\Psi }^{})),$$

regularity at $`𝒪`$ would require $`\mathrm{\Psi }0`$ on $`(\mathrm{},0]`$ (we have used the boundary condition (7) here). Hence $`_{int}`$ is a portion of flat space-time. So avoiding the trivial case implies the existence of a curvature singularity at $`𝒪`$.

Let us now prove the assertion above regarding the Cauchy horizon.

###### Proposition 1

Let $`x_c`$ be the first positive root of $`G(x):=Fe^\psi =\frac{1}{x}`$, if such exists and $`x_c=+\mathrm{}`$ otherwise. Then there are no future pointing outgoing radial null curves emanating from $`𝒪`$ in the region $`x(0,x_c)`$.

Proof The outgoing radial null curves of (1) satisfy

$$\frac{dr}{dv}=F(x)e^\psi (x)=G(x).$$
(8)
Let $`(v_i,r_i)`$ be a point on a solution curve $`\gamma _i`$ of (8) in the region $`0<x<x_c`$. Then $`x_i=v_i/r_i<x_c`$, and so $`G(x_i)<\frac{1}{x_i}`$. If $`x_c`$ is finite, we note that $`x=x_c`$ is a solution of (8),and so by uniqueness, $`\gamma _i`$ cannot cross $`x=x_c`$ away from $`𝒪`$ i.e. for $`v>0`$. Thus

$$\frac{dx}{dv}|_{\gamma _i}=\frac{1}{r}(1xG)>0$$

for $`v(0,v_i]`$. Note that this inequality is immediate when $`x_c=+\mathrm{}`$. So the inequality applies generally and says that as $`v0`$, $`x`$ decreases and is bounded below by $`0`$. Hence the limit

$$x_l=\underset{v0}{lim}x(v)|_{\gamma _i}<x_c$$

exists and is non-negative. Thus either $`rr_{}>0`$ as $`v0`$ \- in which case the singularity is avoided - or $`r0`$ in the limit. In this case,

$`x_l`$ $`=`$ $`\underset{v0}{lim}{\displaystyle \frac{v}{r}}`$
$`=`$ $`\underset{v0}{lim}{\displaystyle \frac{1}{r^{}(v)}}`$
$`=`$ $`\underset{v0}{lim}{\displaystyle \frac{1}{G(x)}}={\displaystyle \frac{1}{G(x_l)}},`$
where all limits are taken along $`\gamma _i`$ and l’Hopital’s rule is used in the second line. The conclusion that $`x_l<x_c`$ is a root of $`xG=1`$ contradicts minimality of $`x_c`$ and completes the proof. $`\mathrm{}`$

###### Corollary 1

If $`G(x)<\frac{1}{x}`$ for all $`x>0`$, then the singularity $`𝒪`$ is censored. $`\mathrm{}`$

###### Corollary 2

If $`G(x)=\frac{1}{x}`$ for some values of $`x>0`$, then $`x=x_c`$ is the Cauchy horizon $``$ of the space-time, where $`x_c`$ is the smallest positive root of $`xG=1`$. $`\mathrm{}`$

These results show an advantage of describing self-similar collapse in the co-ordinates $`v`$ and $`r`$: the visibility of the singularity at $`𝒪`$ (and indeed the presence of an apparent horizon $`F=0`$) can be read-off from the metric. More accurately, the presence of a naked singularity can be determined by tracking the evolution of metric functions, and without having to integrate geodesic equations.

An apparent horizon may form either before or after the Cauchy horizon. This horizon must be space-like, and the region lying to its future is trapped:

###### Proposition 2

If $`F(x_a)=0`$ for some $`x_a>0`$, then $`x=x_a`$ is space-like and the region $`xx_a`$ is trapped or marginally trapped.

Proof Restricting to $`\mathrm{\Sigma }_a:\{v=x_ar\}`$ in (1) gives

$$ds^2|_{\mathrm{\Sigma }_a}=2x_ae^{\psi (x_a)}(1x_aG(x_a))dr^2+r^2d\mathrm{\Omega }^2,$$

which has spatial signature at $`G=F=0`$ when $`x_a>0`$. From (3), we see that $`F^{}0`$ at an apparent horizon. Hence $`F(x)0`$ for $`xx_a`$. $`\mathrm{}`$

We conclude this section with a lemma which will play a central role in determining the stability of $``$ with respect to scalar radiation.

###### Lemma 1

$`G^{}<0`$ prior to the formation of a Cauchy horizon.

Proof We note first that the results of Propositions 1 and 2 show that $`0<xG<1`$ for $`x(0,x_c)`$. Then (2) gives

$$xGF\psi ^{}=xF^2e^\psi \psi ^{}>F\psi ^{},$$

and using (3) we get

$$F^{}+F\psi ^{}<F^{}+xF^2e^\psi \psi ^{}0,$$

i.e. $`G^{}(x)<0`$ for $`x(0,x_c)`$. $`\mathrm{}`$

###### Corollary 3

$`G^{}(x_c)0`$. $`\mathrm{}`$

We note that if $`G^{}(x_c)=0`$, then $`R_{ab}k^ak^b|_{}=0`$, where $`k^a`$ is tangent to the outgoing radial null direction. This implies that there is no ingoing radiative flux of energy-momentum crossing the Cauchy horizon. We rule out this situation as being physically unrealistic and so we will assume that $`G^{}(x_c)<0`$.

## III Time-like geodesics crossing $`𝒩`$ and $``$.

The stability of the Cauchy horizon will be studied from the point of view of the behaviour of the flux of scalar radiation measured by an observer crossing the horizon. This flux is $`=v^a_a\mathrm{\Phi }`$, where $`\mathrm{\Phi }`$ is the scalar field and $`v^a`$ is the unit tangent to an arbitrary time-like geodesic. Thus we will need to determine the behaviour of the tangent $`v^a`$ for such arbitrary geodesics at the Cauchy horizon. Since we will impose boundary conditions on $`\mathrm{\Phi }`$ in terms of the fluxes at $`𝒩`$, we will need to do the same at this surface. The full set of equations governing time-like geodesics may be written in the form

$`\ddot{v}{\displaystyle \frac{1}{r}}(x(G^{}+G\psi ^{})\psi ^{})\dot{v}^2{\displaystyle \frac{e^\psi }{r^3}}L^2,`$ $`=`$ $`0`$ (9)
$`2Ge^\psi \dot{v}^2+2e^\psi \dot{v}\dot{r}+{\displaystyle \frac{L^2}{r^2}}`$ $`=`$ $`1,`$ (10)
$`\dot{\mathrm{\Omega }}`$ $`=`$ $`{\displaystyle \frac{L}{r^2}},`$ (11)
where the overdot represents differentiation with respect to proper time $`\tau `$, $`L`$ is the conserved angular momentum and $`\mathrm{\Omega }`$ is an azimuthal angular variable. (11) plays no further role below, but is given for completeness. It is convenient to rewrite (9) and (10) as a first order system. Defining $`X=(r,v,u)^T`$ where $`u:=\dot{v}`$, these equations may be written as

$$\dot{X}=H(X)=\left(\begin{array}{c}\frac{e^\psi }{2u}(2Ge^\psi u^2\frac{L^2}{r^2}1)\\ u\\ \frac{1}{r}(x(G^{}+G\psi ^{})\psi ^{})u^2+\frac{L^2}{r^3}e^\psi \end{array}\right).$$
(12)
A future-pointing time-like geodesic crossing $`𝒩`$ corresponds to a solution of (12) with initial values $`r_0>0`$, $`v_0=0`$, $`u_0>0`$. The assumptions of the previous section indicate that $`H`$ is $`C^1`$ in a neighbourhood of $`(r_0,v_0,u_0)^3`$, and so standard theorems imply the existence of a $`C^1`$ solution for $`X`$ which exists for (at least) finite duration. Note that this implies that both $`v`$ and $`r`$ (via (10)) are $`C^2`$ functions of proper time $`\tau `$ in a neighbourhood of $`𝒩`$. Thus we can apply Taylor’s theorem and write def1 (11)

$`v(\tau )`$ $`=`$ $`u_o\tau +v_2\tau ^2+O(\tau ^3),`$
$`r(\tau )`$ $`=`$ $`r_0+r_1\tau +O(\tau ^2),`$
where the coefficients $`v_2,r_1`$ can be given in terms of the initial data and metric functions and we have set $`\tau =0`$ at $`𝒩`$. From this we may write down the following result which will be required below.

###### Proposition 3

For any future-pointing time-like geodesic crossing $`𝒩`$, we have

$`v`$ $``$ $`u_0\tau ,\dot{v}u_0,`$ (13)
$`x`$ $``$ $`{\displaystyle \frac{u_0}{r_0}}\tau ,\dot{x}{\displaystyle \frac{u_0}{r_0}},`$ (14)
as $`\tau 0`$ where $`\tau =0`$ on the geodesic at $`𝒩`$. $`\mathrm{}`$

Obtaining equivalent results at the Cauchy horizon is more difficult, as this corresponds to a singular point of the geodesic equations. Two things must be established: the existence of time-like geodesics crossing the horizon and the limiting values of the components of the tangent vector at the horizon. The proof below requires an assumption on the level of differentiability at the horizon which it would be desirable to remove.

###### Proposition 4

Suppose that $`G`$ and $`\psi `$ are differentiable at $`x=x_c`$. Then all radial time-like geodesics whose initial points are sufficiently close to the Cauchy horizon will cross the horizon in finite time. For any time-like geodesic crossing the horizon, the components of the tangent $`\dot{x}`$ and $`\dot{v}`$ have finite non-zero values at the horizon which, denoting them by $`\dot{x}_c`$ and $`\dot{v}_c`$ respectively, satisfy the relation

$$\dot{v}_c=\frac{1}{2}\frac{x_c^2}{v_c\dot{x}_c}e^{\psi _c},$$
(15)
where the subscript refers to the value of a quantity at $`x=x_c`$.

Proof (i) First, we establish a first order non-autonomous system for the geodesics. If $`\xi ^a`$ is the homothetic Killing vector field and $`u^a`$ is tangent to a time-like geodesic, then

$$\frac{d}{d\tau }(\xi _au^a)=1,$$

where $`\tau `$ is proper time along the geodesic (see e.g. Appendix C of wald (12)). Integrating yields

$$g_{ab}\xi ^au^b=2\frac{v}{x}e^\psi (1xG)\dot{v}\frac{v^2}{x^2}e^\psi \dot{x}=k\tau ,$$

for some $`k`$ which is constant along the geodesic. Combining with

$$\frac{2}{x}e^\psi (1xG)\dot{v}^2\frac{2}{x^2}e^\psi v\dot{v}\dot{x}=(1+\frac{L^2}{r^2})$$
(16)
(which is (10) written in terms of $`v`$ and $`x`$) we obtain the first order system

$`\dot{x}`$ $`=`$ $`\pm {\displaystyle \frac{x^2}{v^2}}e^\psi Y^{1/2},`$ (17)
$`\dot{v}`$ $`=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle \frac{x}{v}}e^\psi (1xG)^1((\tau k)\pm Y^{1/2}),`$ (18)
where

$$Y=(\tau k)^2+2\frac{v^2}{x}(1xG)e^\psi (1+\frac{L^2}{r^2}).$$

We choose the upper sign, which corresponds to future-pointing geodesics.

(ii) For radial $`(L=0)`$ time-like geodesics, (9) becomes

$$\ddot{v}=(x^2G^{}+x\psi ^{}(xG1))\frac{\dot{v}^2}{v}.$$

Since $`G^{}(x_c)<0`$, the coefficient on the right hand side is negative for values of $`x<x_c`$ sufficiently close to $`x_c`$. Hence a geodesic with initial value $`x_0=x(\tau _0)`$ sufficiently close to $`x_c`$ in this sense satisfies $`\ddot{v}<0`$ for $`\tau \tau _0`$, and so $`v`$ cannot diverge to infinity in finite time.

(iii) Next, we establish that if $`v\mathrm{}`$ as $`\tau \mathrm{}`$ along a geodesic which does not cross the Cauchy horizon, then $`xx_c`$ as $`\tau \mathrm{}`$. From (16), we can write

$$\frac{\dot{v}}{v}<\frac{\dot{x}}{x(1xG)}.$$

Integrating both sides yields and taking $`v=v_0`$ at some $`0<x_0<x_c`$, we get

$$v<v_0\mathrm{exp}\left(_{x_0}^x\frac{dy}{y(1yG)}\right)$$

for $`x_0<x`$. Thus if $`v`$ diverges to $`+\mathrm{}`$, then so too must the integral. This can only occur if the integrand diverges, i.e. if $`xx_c`$. Now we show that provided a geodesic has initial point sufficiently close to $`x=x_c`$, it cannot behave in this way.

(iv) Consider a radial time-like geodesic for which $`v\mathrm{}`$ and $`xx_c`$ as $`\tau \mathrm{}`$. We have from (9)

$$\frac{v\ddot{v}}{\dot{v}^2}x_c^2G^{}(x_c)=k^2<0$$

as $`\tau \mathrm{}`$. Integrating and reusing this relation yields the asymptotic relations

$`\dot{v}`$ $``$ $`c_1v^{k^2},`$ (19)
$`\ddot{v}`$ $``$ $`k^2c_1v^{2k^21}`$ (20)
as $`\tau \mathrm{}`$ for some $`c_1>0`$. Using these and (16) we obtain

$`\dot{x}`$ $``$ $`c_2v^{k^21}`$ (21)
as $`\tau \mathrm{}`$ for some $`c_2>0`$. We must have $`\dot{x}0`$ as $`\tau \mathrm{}`$, for otherwise $`\dot{x}`$ is positive and bounded away from zero for an infinite amount of time and so $`x`$ reaches $`x_c`$ in finite time. Our present assumption is that this does not happen, so we must have $`k^2<1`$.

The geodesic equations yield

$$2A\ddot{v}+vB\ddot{x}+2A^{}\dot{x}\dot{v}+vB^{}\dot{x}^2=0,$$
(22)
where

$$A(x)=\frac{2}{x}e^\psi (1xG),B=\frac{2}{x^2}e^\psi .$$

Using the assumption that these terms are differentiable at the Cauchy horizon, we have from (19)-(21), def1 (11)

$`2A\ddot{v}`$ $``$ $`2Ak^2c_1^2v^{2k^21}=o(v^{2k^21}),`$
$`2A^{}\dot{x}\dot{v}`$ $``$ $`2A^{}c_1c_2v^1=O(v^1),`$
$`vB^{}\dot{x}^2`$ $``$ $`B^{}v^{2k^21}.`$
Comparing these with (22), we see that we must have

$$\underset{v\mathrm{}}{lim}v\ddot{x}=\underset{v\mathrm{}}{lim}v\frac{B^{}}{B}\dot{x}^2.$$

We have

$$\frac{B^{}}{B}=2x^1+\psi ^{}.$$

Using the energy condition (2), we see that this term is strictly negative at the Cauchy horizon, and so this implies that $`\ddot{x}`$ is positive for sufficiently large values of $`\tau `$. However this contradicts the fact that $`\dot{x}(\tau )>0`$ with $`\dot{x}0`$ as $`\tau 0`$. Hence the geodesic cannot extend to arbitrarily large values of $`v`$ without first crossing the Cauchy horizon.

(v) To conclude, it has been established that (at least some) radial geodesics cross the horizon in finite time and so with a finite value $`v_c`$ of $`v`$. For any such geodesic, including non-radials, we can read off from (17) the non-zero and finite limiting value of $`\dot{x}`$ and from (16) we obtain (15). $`\mathrm{}`$

## IV The scalar field on the Cauchy horizon.

Now we are in a position to examine the stability of the Cauchy horizon by measuring the flux of the scalar field in different regions of the space-time.
In order to measure the flux $`=v^a_a\mathrm{\Phi }`$ we need first the solution of the scalar wave equation,

$$\mathrm{}\mathrm{\Phi }=(g)^{\frac{1}{2}}_a\left[(g)^{\frac{1}{2}}g^{ab}_b\mathrm{\Phi }\right]=0.$$

We exploit the spherical symmetry of the space time and split the scalar field,

$$\mathrm{\Phi }(v,x,\theta ,\varphi )=T(v,x)A(\theta ,\varphi ),$$

where we use the advanced null coordinate $`v`$, the homothetic coordinate $`x`$, and the standard angular coordinates $`\theta ,\varphi `$. Then the line element in these co-ordinates reads

$$ds^2=2e^\psi \left(\frac{1}{x}G\right)dv^2\frac{2e^\psi v}{x^2}dvdx+\frac{v^2}{x^2}d\mathrm{\Omega }^2.$$

By using separation of variables we arrive at a p.d.e. in $`v,x`$

$`2x^2({\displaystyle \frac{1}{x}}G)T,_{xx}+2vT,_{xv}2x^2G^{}T,_x`$
$`2vT,_v\rho e^\psi T=0`$ (23)
where $`\rho =l(l+1)`$ is the separation constant, $`l=0,1,2\mathrm{}`$ is the multipole mode number, and  denotes differentiation w.r.t.$`x`$. The complementary p.d.e. in $`\theta ,\varphi `$ reduces to a form of Legendre’s equation and is solved by the spherical harmonic functions, $`P_l^m(\theta ,\varphi )`$.
We can perform a Mellin transformation on the field, defined by

$$M\{T\}(x,n)=H_n(x)=_0^{\mathrm{}}T(v,x)v^{n1}𝑑v$$

which amounts to replacing $`T(v,x)`$ with $`v^nH_n(x)`$, where $`n`$ is an as yet unconstrained complex parameter. Equation (23) thus reduces to an o.d.e. in $`H(x)`$,

$`2x^2\left({\displaystyle \frac{1}{x}}G\right)H^{\prime \prime }+(2n2x^2G^{})H^{}`$
$`\left({\displaystyle \frac{2n}{x}}+\rho e^\psi \right)H=0`$ (24)
where we have suppressed the subscript $`n`$. Performing the inverse Mellin transform on the solution of this o.d.e. over a contour in the viable range of $`n`$ will return the solution to (23).
This o.d.e. has a number of singular points, namely $`x=0`$ and the roots of $`xG=1`$, the lowest of which we have defined to be $`x_c`$. The canonical form of a second order linear o.d.e. in a neighborhood of $`x=x_0`$ is

$$(xx_0)^2H^{\prime \prime }+(xx_0)q(x)H^{}+p(x)H=0,$$

and when we write equation (24) in its canonical form in the neighborhood of $`x=0`$, we find

$`q(x)={\displaystyle \frac{nx^2G^{}}{1xG}},p(x)={\displaystyle \frac{2n+\rho e^\psi x}{2(1xG)}}.`$
Since $`q(x)`$ and $`p(x)`$ are both $`C^1`$ in a neighborhood of $`x=0`$ we can use the method of Frobenius to solve (24) on $`𝒩`$ foot (13) (see e.g. Chapter 3 of BendOr (14)). The indicial exponents are $`1,n`$.
As it stands we cannot make any assumptions about $`n`$, however later analysis shows if $`\text{Re}(n)1`$ the flux of the scalar field will be always infinite on $`𝒩`$, thus we only consider $`\text{Re}(n)<1`$.
It is possible for $`1`$ and $`n`$ to differ by an integer and so the method of Frobenius yields the following expression for the solution to (24) in a neighborhood of $`x=0`$,

$`H(x)=c_1{\displaystyle \underset{m=0}{\overset{\mathrm{}}{}}}a_mx^{m+1}`$
$`+c_2\left\{k\mathrm{ln}x{\displaystyle \underset{m=0}{\overset{\mathrm{}}{}}}a_mx^{m_+1}+{\displaystyle \underset{m=0}{\overset{\mathrm{}}{}}}b_mx^{mn}\right\}.`$ (25)
In this expression, $`c_1`$ and $`c_2`$ are arbitrary constants, $`a_0=b_0=1`$ with $`k=0`$ if $`1`$ and $`n`$ do not differ by an integer, $`a_0=1,b_0=0`$ with $`k=1`$ if $`1`$ and $`n`$ are equal, and $`a_0=b_0=1`$ with $`k0`$ if $`1+n=m`$ for some positive integer $`m`$.
After some rearranging and some cancellations, the expression for the flux on $`𝒩`$ is

$`_1(v,r)=\dot{v}{\displaystyle \underset{m=0}{\overset{\mathrm{}}{}}}a_m(m+n+1){\displaystyle \frac{v^{n+m}}{r^{m+1}}}`$
$`\dot{r}{\displaystyle \underset{m=0}{\overset{\mathrm{}}{}}}a_m(m+1){\displaystyle \frac{v^{m+n+1}}{r^{m+2}}}`$ (26)
$`_2(v,r)=\dot{v}{\displaystyle \underset{m=0}{\overset{\mathrm{}}{}}}b_{m+1}(m+1){\displaystyle \frac{v^m}{r^{mn+1}}}`$
$`\dot{r}{\displaystyle \underset{m=0}{\overset{\mathrm{}}{}}}b_m(mn){\displaystyle \frac{v^m}{r^{mn+1}}}`$
$`\dot{r}k{\displaystyle \underset{m=0}{\overset{\mathrm{}}{}}}\left[1+(m+1)\mathrm{ln}\left({\displaystyle \frac{v}{r}}\right)\right]a_m{\displaystyle \frac{v^{n+m+1}}{r^{m+2}}}`$
$`+\dot{v}k{\displaystyle \underset{m=0}{\overset{\mathrm{}}{}}}\left[1+(m+n+1)\mathrm{ln}\left({\displaystyle \frac{v}{r}}\right)\right]a_m{\displaystyle \frac{v^{m+n}}{r^{m+1}}}.`$ (27)
where the $`1`$ subscript denotes the $`c_1`$ part, and likewise the $`2`$ subscript.
The components of the velocity have been shown to be finite on $`𝒩`$ in Proposition 3, and we see that for the flux to have a finite measure on $`𝒩`$, that is when $`v=0`$, we require

$$\text{Re}(n)>0.$$

Under this condition we let the scalar field evolve towards $``$, and examine its flux there.
Note: A scalar field coming from past null infinity will have a finite flux thereon if $`\text{Re}(n)1`$. While this physically desirable condition should be imposed, it does not play any role in later analysis.
When we write (24) in its canonical form around $`x=x_c`$, we find the coefficients are

$`q(x)=\left({\displaystyle \frac{xx_c}{x}}\right)\left({\displaystyle \frac{nx^2G^{}}{1xG}}\right)`$
$`p(x)=\left({\displaystyle \frac{xx_c}{x}}\right)^2\left({\displaystyle \frac{n+\frac{\rho e^\psi x}{2}}{1xG}}\right).`$
Now we reach an important distinction, whether $`G(x)`$ has a unique lowest root or multiple lowest roots. We distinguish the two cases so:

###### Lemma 2

When $`xG=1`$ has a unique lowest root,

$`x_c^2G^{}(x_c)+1>0.`$
When $`xG=1`$ has a multiple lowest root,

$`x_c^2G^{}(x_c)+1=0.`$
$`\mathrm{}`$

The two cases will lead to very different analyses, thus we treat them separately.
(i) The first case leads to $`q(x),p(x)`$ being $`C^1`$ on $`x=x_c`$, thus $`x_c`$ is a regular singular point and hence we can use the method of Frobenius. The indicial exponents are $`0,1q_0`$ where

$`q_0={\displaystyle \frac{x_c^2G^{}(x_c)n}{x_c^2G^{}(x_c)+1}}.`$
Since $`n>0`$, Lemma $`1`$ and $`2`$ tell us $`q_0<0`$, hence $`1q_0>0`$, which gives us

$`H(x)=C_1{\displaystyle \underset{m=0}{\overset{\mathrm{}}{}}}A_m\zeta ^{m+1q_0}`$
$`+C_2\left\{k\mathrm{ln}\zeta {\displaystyle \underset{m=0}{\overset{\mathrm{}}{}}}A_m\zeta ^{m+1q_0}+{\displaystyle \underset{m=0}{\overset{\mathrm{}}{}}}B_m\zeta ^m\right\}`$ (28)
where $`\zeta =xx_c`$, and the coefficients have the same structure as (25). From this we calculate the flux,

$`_1=\dot{x}v^nC_1{\displaystyle \underset{m=0}{\overset{\mathrm{}}{}}}(m+1q_0)A_m\zeta ^{mq_0}`$
$`+\dot{v}nv^{n1}C_1{\displaystyle \underset{m=0}{\overset{\mathrm{}}{}}}A_m\zeta ^{m+1q_0}`$ (29)
$`_2=\dot{x}v^nC_2[k{\displaystyle \underset{m=0}{\overset{\mathrm{}}{}}}A_m[\mathrm{ln}\zeta (m+1q_0)+1]\zeta ^{mq_0}`$
$`+{\displaystyle \underset{m=0}{\overset{\mathrm{}}{}}}B_mm\zeta ^{m1}]`$
$`+\dot{v}nv^{n1}C_2\left[{\displaystyle \underset{m=0}{\overset{\mathrm{}}{}}}B_m\zeta ^m+k\mathrm{ln}\zeta {\displaystyle \underset{m=0}{\overset{\mathrm{}}{}}}A_m\zeta ^{m+1q_0}\right].`$ (30)
Using the finiteness of $`\dot{v},\dot{x}`$ given in Proposition 4, we see that if $`q_0<0`$, that is if $`n>0`$, this expression is finite on $``$, i.e. when $`xx_c=\zeta =0`$.
Thus in the case of $`xG=1`$ having a unique lowest root, a scalar field measuring a finite flux entering the region will measure a finite flux on the Cauchy horizon.
(ii) If $`x_c^2G^{}(x_c)+1=0`$, $`x_c`$ is an irregular singular point of (24). Note that this is a special case which one would expect to correspond to a set of measure zero in the class of space-times under consideration. We label $`\eta =x_cx`$ and examine solutions to the o.d.e. in the asymptotic limit $`\eta 0`$ (see e.g. Chapter 3 of BendOr (14)).
We assume the solution to (24) can be written in the form

$$H(\eta )=e^{S(\eta )},$$

reducing (24) to an o.d.e. in $`S`$. Now we assume the common property near irregular singular points,

$$\ddot{S}=o(\dot{S}^2),\eta 0$$

where the overdot denotes differentiation with respect to $`\eta `$. (24) becomes a quadratic in $`\dot{S}`$,

$`\dot{S}^2\left\{(x_c\eta )(x_c\eta )^2G\right\}\left(n+(x_c\eta )^2\dot{G}\right)\dot{S}`$
$`{\displaystyle \frac{n}{x_c\eta }}+{\displaystyle \frac{\rho e^\psi }{2}},\eta 0`$ (31)
If we consider $`xG=1`$ to have a lowest root of multiplicity $`k`$, then we can write its Taylor series around $`\eta =0`$ as

$$1(x_c\eta )G(\eta )=P(\eta )=\eta ^k\frac{P^{(k)}(0)}{k!}+O(\eta ^{k+1})$$

This means if the lowest root is of multiplicity $`k`$, we need the metric functions to be $`C^k`$. This is not too much of a restriction however, since the class of functions with roots of multiplicity $`k`$ becomes very small as $`k`$ increases, meaning we are dealing with a very special case in this analysis.
We can make the approximation

$`n+(x_c\eta )^2\dot{G}`$ $``$ $`n+1,\eta 0,`$
and since we assume the metric coefficients are at least $`C^2`$, we can approximate $`e^\psi `$ by the first term in its expansion, $`c_0`$, in the limit $`\eta 0`$. Thus we arrive at a quadratic in $`S`$,

$`\eta ^k(\dot{S})^2\alpha \dot{S}\beta ,\eta 0,`$ (32)
$`\alpha ={\displaystyle \frac{k!(n+1)}{x_cP^{(k)}(0)}},\beta ={\displaystyle \frac{k!}{x_cP^{(k)}(0)}}\left({\displaystyle \frac{n}{x_c}}+{\displaystyle \frac{\rho c_0}{2}}\right)`$
where $`\alpha ,\beta >0`$ (if $`Re(n)>0`$) and constant in the limit $`\eta 0`$, and $`k>1`$. This quadratic has two solutions corresponding to two linearly independent solutions of (24), which are

$`S_1`$ $``$ $`{\displaystyle \frac{\alpha }{(k1)}}\eta ^{1k}+O(\eta )`$
$`S_2`$ $``$ $`{\displaystyle \frac{\beta }{\alpha }}\eta +{\displaystyle \frac{\beta ^2}{\alpha ^3}}{\displaystyle \frac{\eta ^{k+1}}{(k+1)}}+O(\eta ^{2k+1}),\eta 0.`$
At this point we verify our earlier assumption, namely

$`\ddot{S}=o(\dot{S}^2),\eta 0.`$
Thus we have constructed two solutions to (24),

$`H_1(\eta )`$ $`=`$ $`\eta ^k\mathrm{exp}\left\{{\displaystyle \frac{\alpha }{(k1)}}\eta ^{1k}+O(\eta )\right\}`$ (33)
$`H_2(\eta )`$ $`=`$ $`\mathrm{exp}\{O(\eta )\}`$ (34)
Both of these functions and their derivatives are finite in the limit $`\eta 0`$, $`xx_c`$ if $`Re(n)>0`$, and thus the resulting expressions for the flux are finite, where again we use Proposition 4.
We summarize thus:

###### Proposition 5

Let space-time ($``$,g) satisfy the requirements of Section II and admit a Cauchy horizon $`x=x_c`$. Assume also that $`g_{ab}C^2`$ at $`x=x_c`$. Then a scalar field which has a finite flux on $`𝒩`$, the past null cone of $`𝒪`$, will also have a finite flux on the Cauchy horizon, $``$.

## V Conclusions

We have shown that the Cauchy horizon formed by collapse in a self similar, spherically symmetric space-time is stable with respect to scalar radiation. This space-time is very general, the only other constraints being that the field satisfies the dominant energy condition, and, other than the special case discussed in Section IV(ii), we require the metric functions to be $`C^2`$ on $`𝒩`$ and $``$. These differentiability conditions are stronger than one would like to assume (cf. the $`C^0`$ Cauchy horizons appearing in the collapse of wave maps bizon (5)), but are as low as one can go without having to resort to a generalised solution concept for the wave equation.
The next step is to examine whether linear perturbations of the metric functions will lead to an unstable Cauchy horizon, as is seen, for example, in the Reissner-Nordström solution. Such an examination would be more significant in considering cosmic censorship. Is it difficult to anticipate the general outcome of such an examination. One expects to observe instability for the case of a massless scalar field christo2 (8), but stability for (some sectors of) perfect fluid collapse harada (7). The present results and the Cauchy horizon stability conjecture would lead one to expect stability in general debbie (15).

## Acknowledgement

This research is supported by Enterprise Ireland grant SC/2001/199.

## Appendix A Energy conditions in spherical symmetry

### A.1 Spherical symmetry

We write the line element in double null coordinates;

$$ds^2=2e^{2f}dudv+r^2d\mathrm{\Omega }^2,$$

where $`f=f(u,v)`$, $`r=r(u,v)`$ and $`d\mathrm{\Omega }^2`$ is the line element on the unit 2-sphere. The non-vanishing Ricci tensor terms are

$`R_{uu}`$ $`=`$ $`2r^1(r_{uu}+2r_uf_u),`$
$`R_{vv}`$ $`=`$ $`2r^1(r_{vv}+2r_vf_v),`$
$`R_{uv}`$ $`=`$ $`2r^1(r_{uv}rf_{uv}),`$
$`R_{\theta \theta }`$ $`=`$ $`\mathrm{csc}^2\theta R_{\varphi \varphi }=2{\displaystyle \frac{E}{r}}+2e^{2f}rr_{uv},`$
where $`E`$ is the Misner-Sharp mass,

$$E=\frac{r}{2}(1+2e^{2f}r_ur_v).$$

Subscripts on $`f,r`$ denote partial derivatives. The only non-vanishing Weyl tensor term is

$`\mathrm{\Psi }_2`$ $`=`$ $`{\displaystyle \frac{1}{3}}{\displaystyle \frac{E}{r^3}}+{\displaystyle \frac{1}{3}}e^{2f}(f_{uv}+r^1r_{uv})`$
$`=`$ $`{\displaystyle \frac{E}{r^3}}{\displaystyle \frac{1}{12}}g^{AB}R_{AB}+{\displaystyle \frac{R_{\theta \theta }}{3r^2}},`$
where $`x^A`$ are co-ordinates in the Lorentzian 2-space.

### A.2 The strong energy condition

Our aim is to write down a set of conditions on the curvature terms listed above which are equivalent to the strong energy condition:

$$R_{ab}v^av^b0$$

for all (future-pointing) causal vectors $`\stackrel{}{v}`$.

#### A.2.1 Null vectors

The radial null directions are $`\delta _u^a,\delta _v^a`$. These give

$$R_{uu}0,R_{vv}0.$$

At any point, the general non-radial null direction may be written as

$$v^a=\alpha \delta _u^a+\beta \delta _v^a+\gamma \delta _\varphi ^a.$$
(35)
The null condition is

$$\alpha \beta =2r^2e^{2f}\gamma ^2.$$

We find

$$R_{ab}v^av^b=\alpha ^2R_{uu}+\beta ^2R_{vv}+2\alpha \beta (2f_{uv}+2\frac{E}{r^3}e^{2f}).$$

This is non-negative for all non-radial null vectors if and only if it is non-negative for all values of $`\alpha ,\beta `$ with $`\alpha \beta >0`$. In turn, this is true if

$$\underset{\mu >0}{\mathrm{min}}Q(\mu )0,$$

where $`\mu =\alpha /\beta >0`$ and

$$Q(\mu )=\mu ^2R_{uu}+2\mu (2f_{uv}+2\frac{E}{r^3}e^{2f})+R_{vv}.$$

If $`R_{uu}=R_{vv}=0`$, this is simply

$$f_{uv}+\frac{E}{r^3}e^{2f}0.$$

If $`R_{uu}=0`$ and $`R_{vv}0`$, the condition is equivalent to

$$2\mu (2f_{uv}+2\frac{E}{r^3}e^{2f})+R_{vv}0$$

for all $`\mu >0`$. This is satisfied if $`f_{uv}+\frac{E}{r^3}e^{2f}0`$. If $`f_{uv}+\frac{E}{r^3}e^{2f}<0`$, then the condition will be violated for sufficiently large values of $`\mu `$ (which can always be chosen). The same holds for $`R_{uu}0,R_{vv}=0`$. Thus if $`R_{uu}R_{vv}=0`$, the strong energy condition holds for null directions if and only if

$$f_{uv}+\frac{E}{r^3}e^{2f}0.$$

So now assume that $`R_{uu}>0`$, $`R_{vv}>0`$. The quadratic $`Q(\mu )`$ has a global minimum at

$$\mu _{}=2R_{uu}^1(f_{uv}+\frac{E}{r^3}e^{2f}),$$

while $`Q(0)=R_{vv}>0`$. Thus $`Q(\mu )>0`$ for $`\mu >0`$ if and only if either $`\mu _{}0`$ or $`\mu _{}>0`$ and $`Q(\mu _{})0`$.

$`\mu _{}0`$ if and only if $`f_{uv}+\frac{E}{r^3}e^{2f}0`$.

$`\mu _{}>0`$ if and only if $`f_{uv}+\frac{E}{r^3}e^{2f}<0`$. In this case,

$`Q(\mu _{})`$ $`=`$ $`4R_{uu}^1(f_{uv}+{\displaystyle \frac{E}{r^3}}e^{2f})^2+R_{vv}0`$
$``$ $`|f_{uv}+{\displaystyle \frac{E}{r^3}}e^{2f}|{\displaystyle \frac{1}{2}}(R_{uu}R_{vv})^{1/2}`$
$``$ $`f_{uv}+{\displaystyle \frac{E}{r^3}}e^{2f}{\displaystyle \frac{1}{2}}(R_{uu}R_{vv})^{1/2}.`$
Combining these results, we can say:

$`R_{ab}v^av^b0`$ for all null $`v^a`$ if and only if

$`R_{uu}`$ $``$ $`0`$ (36)
$`R_{vv}`$ $``$ $`0`$ (37)
$`{\displaystyle \frac{1}{2}}(R_{uu}R_{vv})^{1/2}+f_{uv}+{\displaystyle \frac{E}{r^3}}e^{2f}`$ $``$ $`0.`$ (38)
#### A.2.2 Time-like vectors

Again we write

$$v^a=\alpha \delta _u^a+\beta \delta _v^a+\gamma \delta _\varphi ^a,$$

and we can use the time-like condition $`g_{ab}v^av^b=1`$, so that

$$\gamma ^2=r^2(2e^{2f}\alpha \beta 1)0.$$

So in this case we are minimizing over the set $`\alpha \beta \frac{1}{2}e^{2f}`$. we do this by minimizing over the hyperbola $`\alpha \beta =c`$ and then minimizing over all hyperbolas $`c\frac{1}{2}e^{2f}`$. This yields the conditions above and the extra condition

$$\frac{1}{2}(R_{uu}R_{vv})^{1/2}+f_{uv}r^1r_{uv}0.$$
(39)
### A.3 The weak energy condition

The weak energy condition $`T_{ab}v^av^b0`$ for all causal $`v^a`$ can be written, using Einstein’s equation, as $`R_{ab}v^av^bϵR/2`$, where $`ϵ=g_{ab}v^av^b`$. ($`R=`$ Ricci scalar.) Thus the only extra work to do is for time-like vectors. The algebra involved in the previous section only needs minute changes, and we can show that the weak energy condition is equivalent to (36) - (38) and

$$\frac{1}{2}(R_{uu}R_{vv})^{1/2}+r^1r_{uv}+\frac{E}{r^3}e^{2f}0.$$
(40)
### A.4 The dominant energy condition

This states that for every future-pointing timelike $`v^a`$, the vector $`T^{ab}v_b`$ is non-spacelike and future-pointing, and $`T_{ab}v^av^b0`$. Using the usual general form for $`v^a`$, we again have quadratic inequalities for the parameters $`\alpha `$ and $`\beta `$ which can be treated in the usual way. (The non-spacelike condition is $`g_{ac}T^{ab}T^{cd}v_bv_d0`$; the left hand side is homogeneous of degree 2 in $`(\alpha ,\beta )`$ and so quadratic in $`\mu `$.) The future-pointing condition is simple to examine by assuming that $`u,v`$ increase into the future. The resulting inequalities are

$`R_{uu}`$ $``$ $`0,`$
$`R_{vv}`$ $``$ $`0,`$
$`{\displaystyle \frac{E}{r^3}}e^{2f}+r^1r_{uv}`$ $``$ $`0,`$
$`{\displaystyle \frac{1}{2}}(R_{uu}R_{vv})^{1/2}+{\displaystyle \frac{E}{r^3}}e^{2f}+r^1r_{uv}`$ $``$ $`|f_{uv}r^1r_{uv}|.`$
Using the first three of these, we see that the left hand side of the fourth is non-negative, and so the fourth is equivalent to the two inequalities

$`{\displaystyle \frac{1}{2}}(R_{uu}R_{vv})^{1/2}+{\displaystyle \frac{E}{r^3}}e^{2f}+2r^1r_{uv}f_{uv}`$ $``$ $`0,`$
$`{\displaystyle \frac{1}{2}}(R_{uu}R_{vv})^{1/2}+{\displaystyle \frac{E}{r^3}}e^{2f}+f_{uv}`$ $``$ $`0.`$
Note how (as expected) some of these are the same as some of the strong and weak energy conditions.

### A.5 Summary: covariant form of the energy conditions

The energy conditions are given here in terms that use $`r_{uv}`$ and $`f_{uv}`$. A more transparently covariant form is obtained by using $`R`$ and $`\mathrm{\Psi }_2`$ instead of these two. Then the results are as follows (we note that the signs of $`R_{uu}`$, $`R_{vv}`$ and the term $`e^{4f}R_{uu}R_{vv}`$ are invariants - the last of these is defined in terms of contractions of Ricci with the two invariantly defined radial null directions):

The strong energy condition is equivalent to

$`R_{uu}`$ $``$ $`0,`$ (41)
$`R_{vv}`$ $``$ $`0,`$ (42)
$`{\displaystyle \frac{1}{2}}e^{2f}|R_{uu}R_{vv}|^{1/2}+2{\displaystyle \frac{E}{r^3}}+2\mathrm{\Psi }_2{\displaystyle \frac{R}{12}}`$ $``$ $`0,`$ (43)
$`{\displaystyle \frac{1}{2}}e^{2f}|R_{uu}R_{vv}|^{1/2}+{\displaystyle \frac{E}{r^3}}+\mathrm{\Psi }_2{\displaystyle \frac{R}{6}}`$ $``$ $`0.`$ (44)
The weak energy condition is equivalent to (41), (42), (43) and

$`{\displaystyle \frac{1}{2}}e^{2f}|R_{uu}R_{vv}|^{1/2}+{\displaystyle \frac{E}{r^3}}+\mathrm{\Psi }_2+{\displaystyle \frac{R}{12}}`$ $``$ $`0.`$ (45)
The dominant energy condition is equivalent to (41), (42), (43) and

$`{\displaystyle \frac{E}{r^3}}+\mathrm{\Psi }_2+{\displaystyle \frac{R}{12}}`$ $``$ $`0,`$ (46)
$`{\displaystyle \frac{1}{2}}e^{2f}|R_{uu}R_{vv}|^{1/2}+{\displaystyle \frac{R}{4}}`$ $``$ $`0.`$ (47)