text:
# Kohn-Sham equations for nanowires with direct current

## I Introduction

The goal of building real electronic devices from individual molecules joachim2000  has spurred experimental and theoretical studies of molecular conductivities. Measurements of current-voltage (I-V) characteristics have been reported for single molecules reed97 ; cui2001  but mechanisms for electron transmissions through molecular interfaces remain largely unexplored and unexplained nitzan2001 ; grossmann2002 . It appears that the miniaturization of microelectronic devices shrinks the size of the system to a level when the simulations of transport characteristics of molecular electronic devices can be performed by ab initio electronic structure methods (with single benzene-1,4-dithiolate molecule as the molecular wire, the contact wire-contact “supermolecule” consists of just 32 atoms). However, the main principle hurdle for the application of the first principle electronic structure methods is not the size of the system but conceptually unsolved problems of modeling of the electronic current and interpretation of the applied voltage for correlated electronic systems. In traditional methods, which are based on the Landauer theory landauer70 , the molecule and the edge part of the metal contacts can be in principle treated ab initio, e.g. in the density functional theory (DFT). But then to induce the electronic current it is assumed that the system is attached to two equilibrium electronic reservoirs with the different chemical potentials datta95 . Being computationally not very efficient, since it requires to perform the calculations over the whole coordinate space, these traditional approaches are also pivotally relied on the electronic reservoir concept and often use ad hoc reservoir-molecule interaction as an essential input parameter.

We have previously shown that the reservoir model can be conceptually circumvented in electronic transport calculations if the electronic current is used as an input information instead the applied voltage bias kosov2001 ; kosov2002 . In our approach the electronic current enters a quantum variational problem via the corresponding Lagrange multiplier. In this paper, we continue development of our method and discuss detailed derivation of the Kohn-Sham equations for a nanowire with steady current. We begin by describing the variational principle for an inhomogeneous electron gas with steady current. We then derive the Kohn-Sham equations with steady current.

## II Variational principle

We start our derivation by defining density and current density of an inhomogeneous, interacting $`N`$-electron gas in an external scalar potential $`v(𝐫)`$ as (a.u. are used throughout the paper):

$`\rho (𝐫)`$ $`=`$ $`N{\displaystyle \mathrm{\Psi }(𝐫,𝐫_2,\mathrm{};𝐫_N)\mathrm{\Psi }^{}(𝐫,𝐫_2,\mathrm{},𝐫_N)\underset{i=2}{\overset{N}{}}d𝐫_i},`$ (1)
$`𝐣(𝐫)`$ $`=`$ $`N{\displaystyle \text{Im}\left\{\mathrm{\Psi }^{}(𝐫,𝐫_2,\mathrm{},𝐫_N)\mathrm{\Psi }(𝐫,𝐫_2,\mathrm{},𝐫_N)\right\}\underset{i=2}{\overset{N}{}}d𝐫_i},`$ (2)
where $`\mathrm{\Psi }`$ is the many-electron wave-function. For simplicity of notation we do not write explicitly the spin indeces but our derivation can be readily extended to spin-polarized systems.

Next, we introduce a subsidiary condition for the current density distribution by specifying the quantum wire constraint on the current kosov2001 ; kosov2002 :

$$𝑑y𝑑zj_x(𝐫)=I(x),$$
(3)
where $`I(x)`$ is a given input value of the steady current through the electronic device. With this constraint the current is aligned along the x-axis and the net current flow across a cross section $`𝑑y𝑑zj_x(𝐫)`$ is constrained, and this quantity is available experimentally.

For an inhomogeneous, interacting electron gas in scalar external potential $`v(𝐫)`$, the ground state energy is given by the Hohenberg-Kohn energy functional hohenberg64 :

$$E_{HK}[\rho ]=T_o[\rho ]+𝑑𝐫v(𝐫)\rho (𝐫)+\frac{1}{2}𝑑𝐫𝑑𝐫^{}\frac{\rho (𝐫)\rho (𝐫^{})}{|𝐫𝐫^{}|}+E_{xc}[\rho ].$$
(4)
The first term, $`T_0`$, is the kinetic energy functional of noninteracting electrons with given density and current density distribution. The second and third terms, are the interaction energy with external potential and the electrostatic interactions. The last term, $`E_{xc}[\rho ]`$, is the exchange and correlation energy functional. The current dependent correction to the standard exchange-correlation energy functional is $`O(j^2)`$ vignale87 . A dependence of the exchange-correlation energy functional upon $`𝐣(𝐫)`$ can be readily included in our approach, yielding the additional exchange-correlation vector potential vignale87 :

$$𝐀_{xc}[\rho ,𝐣](𝐫)=\frac{\delta E_{xc}[\rho ,𝐣]}{\delta 𝐣(𝐫)}$$
(5)
which modifies the kinetic energy operator in the final Kohn-Sham equation in the same way as an external vector potential does. The current dependent correction to the standard exchange-correlation energy functional is $`O(j^2)`$ vignale87 . Molecular and atomic wires can sustain a current of about $`10^4a.u.`$ for typical experiments on molecular conduction reed97 . Therefore the direct current-dependent corrections to the total electronic energy are only of the order of $`10^8a.u.`$ and can be ignored in most molecular transport calculations.

Likewise our derivation of the Schrödinger equation for current carrying states kosov2002 , we use the constrained variational principle in our derivations. The problem is to minimize Hohenberg-Kohn energy functional (4) subject to the constraint on the current density (3). The constraint is explicitly achieved by introduction of an auxiliary functional

$$\mathrm{\Omega }[\mathrm{\Psi }]=E_{HK}[\rho ]𝑑xa(x)\left(𝑑y𝑑zj_x(𝐫)I(x)\right).$$
(6)
The first term giving the total energy is standard in the variational formulation of the DFT. The second term with the Lagrange multipliers $`a(x)`$ has been introduced to impose the subsidiary constraint for the current density (3).

## III Kohn-Sham equations with direct current

Following the Kohn and Sham approach kohn65 , we introduce a reference fermion system with orthonormal single-particle orbitals $`\psi _i(𝐫)`$ and occupation numbers $`n_i`$ to reproduce the charge and current densities:

$$\rho (𝐫)=\underset{i}{}n_i\psi _i^{}(𝐫)\psi _i(𝐫),$$
(7)
$$𝐣(𝐫)=\frac{1}{2i}\underset{i}{}n_i\left(\psi _i^{}(𝐫)\psi _i(𝐫)\psi _i(𝐫)\psi _i^{}(𝐫)\right).$$
(8)
According to the variational principle and assuming a fixed set of the occupation numbers $`n_i`$, we carry out the minimization of the Hohenberg-Kohn energy functional with respect to the single-electron orbitals $`\psi _i(𝐫)`$ subject to the constraint for the current density (3), i.e. the variation of auxiliary functional $`\mathrm{\Omega }[\mathrm{\Psi }]`$. The minimization yields the following single-electron self-consistent equations:

$$\left(\frac{1}{2}+v(𝐫)+𝑑𝐫\frac{\rho (𝐫)}{|𝐫𝐫^{}|}+v_{xc}[\rho ](𝐫)+\frac{1}{2i}[_x,a(x)]_+\right)\psi _i(𝐫)=E_i\psi _i(𝐫),$$
(9)
where the anti-commutator term,

$$[_x,a(x)]_+=_x(a(x))+2a(x)_x$$
(10)
is the additional imaginary potential arising directly from the constraint on electronic current (Appendix A). This equation is a particular case of the more general equation derived recently by Kosov and Greer kosov2001 . Relating a position gradient with the corresponding momentum operator ($`\widehat{p}_x=i_x`$) we can rewrite the imaginary potential as the momentum-dependent operator:

$$\frac{1}{2i}[_x,a(x)]_+=\frac{1}{2}[\widehat{p}_x,a(x)]_+.$$
(11)
Since the operator (11) is hermitian the Konh-Sham equation (9) is a hermitian eigenproblem.

We have not specified yet the occupation numbers $`n_i`$ which enter the Kohn-Sham equations (9) via the density and the current density. To fix the set of occupation numbers, we need to know the number of electrons in the nanowire. Although we deal with the open quantum system where electrons can leave and enter the wire, we still can put a constraint on the total number of electrons because of time and space independence of the electronic current (i.e. steady state direct current). Then, the continuity equation tells us that the total number of electrons does not depend upon time. If we assume that the charge neutrality is maintained after the establishment of the current currying state, the average number of electrons in the nanowire $`N`$ is the same as the number of electrons in the zero current wire. Therefore to obtain the occupation numbers we can minimize $`\mathrm{\Omega }`$ (6) with respect to $`n_i`$ subject to the particle conserving constraint imposed with the additional Lagrange multiplier $`\mu `$ (chemical potential):

$$\delta \left[\mathrm{\Omega }\mu (\underset{i}{}n_iN)\right]0.$$
(12)
The direct differentiation of the Hohenberg-Kohn energy is straitforward and results into the following expression (Janak theorem) janak78 :

$`{\displaystyle \frac{E_{HK}}{n_i}}={\displaystyle 𝑑𝐫\psi _i^{}(𝐫)\left(\frac{1}{2}+v(𝐫)+𝑑𝐫\frac{\rho (𝐫)}{|𝐫𝐫^{}|}+v_{xc}[\rho ](𝐫)\right)\psi _i(𝐫)}+`$
$`[{\displaystyle \underset{k}{}}n_k{\displaystyle }d𝐫{\displaystyle \frac{\psi _k^{}(𝐫)}{n_i}}({\displaystyle \frac{1}{2}}+v(𝐫)+{\displaystyle }d𝐫{\displaystyle \frac{\rho (𝐫)}{|𝐫𝐫^{}|}}+v_{xc}[\rho ](𝐫))\psi _k(𝐫)+c.c.]`$
The details of the differentiation of the constrained functional with respect to an occupation number $`n_i`$ are given in Appendix B. The differentiation leads to the following expression for the derivative:

$`{\displaystyle \frac{\mathrm{\Lambda }}{n_i}}`$ $`={\displaystyle \frac{1}{2i}}{\displaystyle 𝑑𝐫\psi _i^{}(𝐫)[_x,a(x)]_+\psi _i(𝐫)}+`$ (13)
$`[{\displaystyle \frac{1}{2i}}{\displaystyle \underset{k}{}}n_k{\displaystyle }d𝐫{\displaystyle \frac{\psi _k^{}(𝐫)}{n_i}}[_x,a(x)]_+\psi _k(𝐫)+c.c.].`$
Combining these two derivatives we arrive to the following expression for the derivative of the auxiliary functional:

$`{\displaystyle \frac{\mathrm{\Omega }}{n_i}}={\displaystyle 𝑑𝐫\psi _i^{}(𝐫)\left(\frac{1}{2}+v(𝐫)+𝑑𝐫\frac{\rho (𝐫)}{|𝐫𝐫^{}|}+v_{xc}[\rho ](𝐫)+\frac{1}{2i}[_x,a(x)]_+\right)\psi _i(𝐫)}+`$
$`[{\displaystyle \underset{k}{}}n_k{\displaystyle }d𝐫{\displaystyle \frac{\psi _k^{}(𝐫)}{n_i}}({\displaystyle \frac{1}{2}}+v(𝐫)+{\displaystyle }d𝐫{\displaystyle \frac{\rho (𝐫)}{|𝐫𝐫^{}|}}+v_{xc}[\rho ](𝐫)+{\displaystyle \frac{1}{2i}}[_x,a(x)]_+)\psi _k(𝐫)`$
$`+c.c.].`$ (14)
Given that $`\psi _i`$ satisfy eq.(9) it is easy to see that:

$$\frac{\mathrm{\Omega }}{n_i}=E_i𝑑𝐫\psi _i^{}(𝐫)\psi _i(𝐫)+\underset{k}{}n_kE_k\frac{}{n_i}𝑑𝐫\psi _k^{}(𝐫)\psi _k(𝐫)=E_i.$$
(15)
From eq.(12) with the use of eq.(15) we arrive to the inequality

$$\underset{i}{}(E_i\mu )\delta n_i0$$
(16)
which gives rise to the standard “build up principle”( orbitals with $`E_i<\mu `$ have $`n_i=1`$ and those with $`E_i>\mu `$ have $`n_i=0`$) which is still hold within our approach: a maximum two electrons are put into current carrying orbitals in the order of increasing orbital energy. The approach can be readily extended to the finite temperature case replacing of the integer $`n_i=0\text{ or }1`$ occupation number by the Fermi-Dirac expression. The formal aspects of the DFT with direct current and non-zero temperature are discussed elsewhere.

The Kohn-Sham equations with current (9) are not yet in a form allowing for a solution to be found as the Lagrange multiplier $`a(x)`$ is not known yet. The additional equation for the Lagrangian multiplier $`a(x)`$ can be obtained if we require that the Kohn-Sham orbitals from the eq.(9) yield the required current $`I(x)`$:

$$\frac{1}{2i}\underset{i}{}n_i𝑑y𝑑z\left(\psi _i^{}(𝐫)_x\psi _i(𝐫)\psi _i(𝐫)_x\psi _i^{}(𝐫)\right)=I(x).$$
(17)
The system of nonlinear equations (9, 17) are to be solved simultaneously and being combined together composes a set of Kohn-Sham equations for nanowire with direct current.

## IV Self-consistent absorbing-emitting boundary conditions

Multiplying the eq. (9) from the left on the $`n_i\psi _i(𝐫)^{}`$, summing up over the Kohn-Sham orbitals and subtracting from the obtained equation its own complex conjugate yields the continuity equation:

$$𝐣(𝐫)=_x(a(x)\rho (𝐫)).$$
(18)
After the integration of eq.(18) over the $`\{y,z\}`$ plane and substitution of the constraint for the current density eq.(3) into the equation the following first order differential equation for the Lagrange multiplier $`a(x)`$ can be deduced:

$$_x(a(x)\rho _{yz}(x)I(x))=0,$$
(19)
where we have introduced the quantity

$$\rho _{yz}(x)=𝑑y𝑑z\rho (𝐫).$$
(20)
The solution of the differential equation (19) depends upon the asymptotic behavior of the specified function for the current $`I(x)`$. It is computationally convenient to assume the space-localized boundary conditions ($`I(x),\rho _{yz}(x)`$ $``$ 0, when $`x\pm \mathrm{}`$). This space-localized boundary conditions yields zero integration constant for the differential equation (19) and we obtain the following expression for the Lagrange multiplier:

$$a(x)=\frac{I_x}{\rho _{yz}(x)}.$$
(21)
With this choice of the Lagrange multiplier (21) we complete the Kohn-Sham equations with the current $`I`$ in the self-consistent and closed form:

$$\left(\frac{1}{2}+v(𝐫)+𝑑𝐫\frac{\rho (𝐫)}{|𝐫𝐫^{}|}+v_{xc}[\rho ](𝐫)\frac{1}{2i}[\frac{I(x)}{\rho _{yz}(x)},\frac{}{x}]_+\right)\psi _i(𝐫)=E_i\psi _i(𝐫).$$
(22)
Some aspects of eq.(22) deserve special discussion. Suppose that the constrained current starts at the left system boundary L and completely absorbed at the right system boundary R:

$$I(x)=I\theta (xL)\theta (Rx),$$
(23)
where $`I`$ is the steady current through the system and $`\theta `$ is the step function. Then one of the terms in the anticommutator imaginary potential becomes proportional to the derivative of the step function and yields the singular imaginary $`\delta `$-function potential on the system boundaries:

$$\frac{Ii}{2\rho _{yz}(x)}(\delta (xL)\delta (xR)).$$

This absorbing-emitting imaginary boundary potential is optimized due to its self-consistency. It depends upon the density and value of the input current and emits the current $`I`$ at the left boundary and completely absorbs the same current at the right boundary.

## V Conclusions

In this paper, we have derived the set of Kohn-Sham equations for a nanowire with direct, steady current. Not restricted to the linear response, our approach uses a constrained minimization of the total energy with a subsidiary condition for the current density to formulate the Kohn-Sham equation with direct current. The subsidiary condition for the current density is maintained during the variation by a pointwise Lagrange multiplier. Being entirely based on the DFT, our method incudes electronic correlations into the transport calculations at the same level as they are taken into account into the corresponding exchange-correlation functionals.

## Appendix A Functional derivatives of the constraint functional with respect to orbitals

In this Appendix we compute the functional derivative of the constraint functional. The constraint functional has the following form

$$\mathrm{\Lambda }[\psi _i,\psi _i^{}]=𝑑xa(x)\left(𝑑y𝑑zj_x(𝐫)I_x\right),$$
(24)
with

$$j_x(𝐫)=\frac{1}{2i}\underset{i}{}n_i\left(\psi _i^{}(𝐫)_x\psi _i(𝐫)\psi _i(𝐫)_x\psi _i^{}(𝐫)\right)$$
(25)
The direct variation of the constraint functional with respect to the $`\psi ^{}(𝐫)`$ yields:

$`{\displaystyle \frac{\delta \mathrm{\Lambda }[\psi ]}{\delta \psi _i^{}(𝐫)}}`$ $`=`$ $`{\displaystyle \frac{1}{2i}}{\displaystyle 𝑑𝐫^{}a(x^{})\frac{\delta }{\delta \psi ^{}(𝐫)}\underset{k}{}n_k\left\{\psi _k^{}(𝐫^{})_x\psi _k(𝐫^{})\psi _k(𝐫^{})_x\psi _k^{}(𝐫^{})\right\}}`$
$`=`$ $`{\displaystyle \frac{1}{2i}}a(x)_x\psi _i(𝐫){\displaystyle \frac{1}{2i}}{\displaystyle \frac{\delta }{\delta \psi _i^{}(𝐫)}}{\displaystyle 𝑑𝐫^{}a(x^{})\underset{k}{}n_k\psi _k(𝐫^{})\psi _k^{}(𝐫^{})}`$
$`=`$ $`{\displaystyle \frac{1}{2i}}a(x)_x\psi _i(𝐫){\displaystyle \frac{1}{2i}}{\displaystyle \frac{\delta }{\delta \psi _i^{}(𝐫)}}{\displaystyle 𝑑𝐫^{}_x\left\{\underset{k}{}n_ka(x^{})\psi _k(𝐫^{})\psi _k^{}(𝐫^{})\right\}}`$
$`+`$ $`{\displaystyle \frac{1}{2i}}{\displaystyle \frac{\delta }{\delta \psi _i^{}(𝐫)}}{\displaystyle 𝑑𝐫^{}\underset{k}{}n_k\psi _k^{}(𝐫^{})_x\left\{a(x^{})\psi _k(𝐫^{})\right\}}`$
$`=`$ $`{\displaystyle \frac{1}{2i}}a(x)_x\psi _i(𝐫)+{\displaystyle \frac{1}{2i}}_x\{a(x)\psi _i(𝐫)\}={\displaystyle \frac{1}{2i}}[a(x),_x]_+\psi _x(𝐫)`$
## Appendix B Differentiation of the constrained functional with respect to an occupation number $`n_i`$

In this appendix the technical details of the direct differentiation of the constraint functional (24) is presented.

$`{\displaystyle \frac{}{n_i}}\mathrm{\Lambda }[\psi _i,\psi _i^{}]={\displaystyle \frac{}{n_i}}{\displaystyle 𝑑xa(x)𝑑y𝑑zj_x(𝐫)}=`$
$`{\displaystyle 𝑑xa(x)j_{ix}(x)}+{\displaystyle \frac{1}{2i}}{\displaystyle 𝑑xa(x)𝑑y𝑑z\underset{k}{}n_k\frac{}{n_i}\left(\psi _k^{}(𝐫)_x\psi _k(𝐫)\psi _k(𝐫)_x\psi _k^{}(𝐫)\right)}`$ (26)
Here we have introduced the new quantity which is the current through the orbital $`i`$ integrated over the plane $`\{y,z\}`$:

$$j_{ix}(x)=\frac{1}{2i}𝑑y𝑑z\left(\psi _i^{}(𝐫)_x\psi _i(𝐫)\psi _i(𝐫)_x\psi _i^{}(𝐫)\right)$$

Then

$`{\displaystyle \frac{}{n_i}}\mathrm{\Lambda }[\psi _i,\psi _i^{}]={\displaystyle 𝑑xa(x)j_{ix}(x)}+{\displaystyle \frac{1}{2i}}{\displaystyle \underset{k}{}}n_k{\displaystyle 𝑑𝐫a(x)\left(\frac{\psi _k^{}(𝐫)}{n_i}_x\psi _k(𝐫)\frac{\psi _k(𝐫)}{n_i}_x\psi _k^{}(𝐫)\right)}+`$
$`{\displaystyle \frac{1}{2i}}{\displaystyle \underset{k}{}}n_k{\displaystyle 𝑑𝐫a(x)\left(\psi _k^{}(𝐫)_x\frac{\psi _k(𝐫)}{n_i}\psi _k(𝐫)_x\frac{\psi _k^{}(𝐫)}{n_i}\right)}`$
$`={\displaystyle 𝑑xa(x)j_{ix}(x)}+{\displaystyle \frac{1}{2i}}{\displaystyle \underset{k}{}}n_k{\displaystyle 𝑑𝐫a(x)\left(\frac{\psi _k^{}(𝐫)}{n_i}_x\psi _k(𝐫)\frac{\psi _k(𝐫)}{n_i}_x\psi _k^{}(𝐫)\right)}+`$
$`{\displaystyle \frac{1}{2i}}{\displaystyle \underset{k}{}}n_k{\displaystyle 𝑑𝐫\left(\frac{\psi _k^{}(𝐫)}{n_i}_x(a(x)\psi _k(𝐫))\frac{\psi _k(𝐫)}{n_i}_x(a(x)\psi _k^{}(𝐫))\right)}=`$
$`{\displaystyle }dxa(x)j_{ix}(x)+{\displaystyle \frac{1}{2i}}[{\displaystyle \underset{k}{}}n_k{\displaystyle }d𝐫{\displaystyle \frac{\psi _k^{}(𝐫)}{n_i}}[_x,a(x)]_+\psi _k(𝐫)c.c.]=`$
$`{\displaystyle \frac{1}{2i}}{\displaystyle }d𝐫a(x)(\psi _i^{}(𝐫)_x\psi _i(𝐫)\psi _i(𝐫)_x\psi _i^{}(𝐫))+{\displaystyle \frac{1}{2i}}[{\displaystyle \underset{k}{}}n_k{\displaystyle }d𝐫{\displaystyle \frac{\psi _k^{}(𝐫)}{n_i}}[_x,a(x)]_+\psi _k(𝐫)c.c.]=`$
$`{\displaystyle \frac{1}{2i}}{\displaystyle }d𝐫\psi _i^{}(𝐫)[_x,a(x)]_+\psi _i(𝐫)+{\displaystyle \frac{1}{2i}}[{\displaystyle \underset{k}{}}n_k{\displaystyle }d𝐫{\displaystyle \frac{\psi _k^{}(𝐫)}{n_i}}[_x,a(x)]_+\psi _k(𝐫)c.c.]`$ (27)


format: text

text:
# 1 Abstract

## 1 Abstract

We propose that the overlapping shock fronts from young supernova remnants produce a locally unsteady, but globally steady large scale spiral shock front in spiral galaxies, where star formation and therefore massive star explosions correlate geometrically with spiral structure. This global shock front with its steep gradients in temperature, pressure and associated electric fields will produce drifts, which in turn give rise to a strong sheet-like electric current, we propose. This sheet current then produces a large scale magnetic field, which is regular, and connected to the overall spiral structure. This rejuvenates the overall magnetic field continuously, and also allows to understand that there is a regular field at all in disk galaxies. This proposal connects the existence of magnetic fields to accretion in disks. We not yet address all the symmetries of the magnetic field here; the picture proposed here is not complete. X-ray observations may be able to test it already.

## 2 Introduction

Magnetic fields are found almost everywhere in the Universe, and usually they are strong, but rarely so strong as to dominate locally, as, e.g., in the environment of pulsars. After magnetic fields had been predicted to exist in the disk of our Galaxy inside the interstellar medium on the basis of cosmic ray arguments, , and the presence of optical polarization from dust scattering, , with a strength predicted to be around 5 microGauss, to within 20 percent, they were indeed detected, with a strength as expected, see, e.g., .

Now the best measurement for the Solar neighborhood is 6 - 7 microGauss, with about 1/3 to 1/2 regular, and 2/3 to 1/2 irregular components, .

### 2.1 Critical observations

The magnetic field is a special challenge to understand in physics terms. Other critical measurements are:

The order of the magnetic field in the Galaxy is destroyed by the interstellar medium stirring on time scales of order 30 million years, as can be deduced from Cosmic Ray transport arguments. This is less than the rotation period of the Galaxy, and so any diffusive process taking many rotation periods is ruled out. In contrast, inside stars there is enough time.

Starburst galaxies such as M82 have a magnetic field, which is ordered and also in near-equipartition. And yet the starburst is just a few tens of millions of years old. The conclusion is once more that the resurrection time scale of the magnetic field is of order 30 million years.

It appears that the magnetic field has an overall orientation, with its direction along the spiral arms inwards always, . This seems to be a general pattern for spiral galaxies, as long as they are not interacting. There is also a different point of view as regards our Galaxy, , but a worry is that local effects overshadow the global picture for us; and so many expect that also in our Galaxy the magnetic field is pointed globally along the spiral arms inwards. However, we have to keep in mind that the observational statistics are still small.

Galaxies which have a strong magnetic field, but also solid body rotation over that radial range of the disk where the magnetic field has been measured, are another important test.

Clusters of galaxies have a magnetic field which normally appears to be quite strong, with a lower limit of several percent of equipartition in energy density, . There are two examples of clusters which seem to have fields close to equipartition, e.g. . Even large scale sheets in the distribution of galaxies have been successfuly measured to have magnetic fields, ; in that case the strength is not very certain, between 10 and 100 nanoGauss range the best estimates, but the fields could also be stronger; just the existence is certain, since that argument is based on the observed synchrotron emission.

### 2.2 History and reviews

The history of this topic is large: The dynamo process was invented twice, by Steenbeck and his associates, in ,  and by Parker, in  and . Recent reviews are in . Large scale magnetic fields are discussed in .

### 2.3 The challenge

On the basis of the data the challenge is to find a mechanism to generate an ordered magnetic field which is fast. We outline a proposal here that may be able to do this. First attempts along these lines were made in . This work reported here is another step along the way, but we are nowhere near finished.

## 3 The equations for mass and current

We follow the arguments of , in the notation of .

Ions have mass $`m_i`$, charge $`Z`$, and velocity $`v_i`$ and electrons mass $`m_e`$ and velocity $`v_e`$, $`𝐄`$ and $`𝐁`$ are the electric and magnetic field. Ions have a density $`n_i`$ of with an average charge of $`Z`$, and electrons have a density $`n_e`$. Then the two equations of motion for ions and electrons are

$`n_im_i({\displaystyle \frac{}{t}}𝐯_𝐢+(𝐯_𝐢)𝐯_𝐢)=`$ (1)
$`+n_im_i\nu _i\mathrm{\Delta }𝐯_𝐢+n_im_i(\nu _i/3+\zeta _i)\mathrm{grad}\mathrm{div}𝐯_𝐢+`$ (2)
$`2n_im_i\omega _𝐏\times 𝐯_𝐢+n_iZe(𝐄+{\displaystyle \frac{𝐯_𝐢}{c}}\times 𝐁)`$ (3)
$`P_in_im_i\mathrm{\Phi }+𝐏_{\mathrm{𝐢𝐞}}`$ (4)
and

$`n_em_e({\displaystyle \frac{}{t}}𝐯_𝐞+(𝐯_𝐞)𝐯_𝐞)={\displaystyle \frac{\sigma _T𝐅n_e}{c}}+`$ (5)
$`+n_em_e\nu _e\mathrm{\Delta }𝐯_𝐞+n_em_e(\nu _e/3+\zeta _e)\mathrm{grad}\mathrm{div}𝐯_𝐞+`$ (6)
$`2n_em_e\omega _𝐏\times 𝐯_𝐞n_ee(𝐄+{\displaystyle \frac{𝐯_𝐞}{c}}\times 𝐁)`$ (7)
$`P_en_em_e\mathrm{\Phi }+𝐏_{\mathrm{𝐞𝐢}}`$ (8)
where $`\omega _𝐏`$ is the pattern velocity of any preferred rotating system of reference, of interest in the case we consider a spiral stationary pattern. We need to point out that these equations describe the forces and so the changes in the orbit of any ion or electron (positron), but they do not contain the steady flow of ions/electrons, which are carried around in a circle perhaps; this can happen in the neighborhood of the intense radiation field of an accretion disk around a black hole: The radiation field pushes on the skin of interstellar material close to it, and so a charged layer is formed; this charged layer is carried around in its orbit with the material, and so a current is built up. This current then produces a specific unique symmetry, as discussed below.

The stress-tensors $`𝚿_𝐞`$ and $`𝚿_𝐢`$ are replaced with viscosity terms and the pressure gradient, $`\mathrm{\Phi }`$ is the gravitational potential, $`𝐅`$ the radiative flux, acting only on electrons via Thomson scattering (cross section $`\sigma _T`$), and $`𝐏_{\mathrm{𝐞𝐢}}=𝐏_{\mathrm{𝐢𝐞}}`$ is the frictional exchange of momentum between ions and electrons, while the term involving the viscosity $`\nu _i`$ and $`\nu _e`$ allows for accretion disks.

We assume, that we can ignore the dependence of the density in the viscosity terms, and that the viscosity is the same for electrons and ions to first order, since it clearly depends on large scale flows, and at the level at which we are initially interested in, not on the microphysics. Therefore we set $`\nu _i=\nu _e=\nu `$ for the kinematic shearflow viscosity, and also $`\zeta _i=\zeta _e=\zeta `$ for the specific volume viscosity. These last assumptions are made here only for initial simplicity.

We define the density $`\rho `$, flow velocity $`𝐯_𝐢`$, $`\sigma _{i,e}`$ the net charge density, and electric current $`𝐣`$ in an obvious way,

$`\rho =`$ $`n_im_i+n_em_e`$ (9)
$`\rho 𝐯=`$ $`n_im_i𝐯_𝐢+n_em_e𝐯_𝐞`$ (10)
$`𝐣=`$ $`e(n_iZ𝐯_𝐢n_e𝐯_𝐞)`$ (11)
$`\sigma _{i,e}=`$ $`n_iZn_e`$ (12)
$`n_i𝐯_𝐢m_i(1+Z{\displaystyle \frac{m_e}{m_i}})=`$ $`\rho 𝐯+{\displaystyle \frac{m_e}{e}}𝐣`$ (13)
$`n_e𝐯_𝐞{\displaystyle \frac{m_i}{Z}}(1+Z{\displaystyle \frac{m_e}{m_i}})=`$ $`\rho 𝐯{\displaystyle \frac{m_i}{Ze}}𝐣`$ (14)
Cosmic Ray particles can be part of these equations, but due to the two-stream instability  they cannot flow faster than the local Alfvén speed. There is a possibility to create secondaries, electron-positron pairs from charged pion decay, arising from p-p collisions between cosmic ray particles and the interstellar medium; another way to consider this to allow for charge exchanges, recombination and ionization: Let then $`S_e/n_e`$ be the rate of such processes. $`S_i`$ and $`S_e`$ is the rate of density increase of ions and electrons, respectively.

The Maxwell equations are

$`\mathrm{div}𝐄=`$ $`4\pi e\sigma _{i,e}`$ (15)
$`\mathrm{div}𝐁=`$ $`0`$ (16)
$`\mathrm{curl}𝐄=`$ $`{\displaystyle \frac{1}{c}}{\displaystyle \frac{}{t}}𝐁`$ (17)
$`\mathrm{curl}𝐁=`$ $`{\displaystyle \frac{4\pi }{c}}𝐣+{\displaystyle \frac{1}{c}}{\displaystyle \frac{}{t}}𝐄`$ (18)
The Maxwell equations allow to relate observed magnetic fields to the underlying symmetries, and electric currents. In the following we ignore the electric field inside a shockwave at first. We add the two main equations, to obtain the standard equation of motion, including magnetic fields and currents, and also subtract these equations to obtain the equation of motion for electric currents.

### 3.1 Mass and momentum flow

We consider the galactic disk, in the observer’s inertial frame system, i.e. keep the rotational flow $`𝐯_\varphi `$ with all spatial gradients. The equations allow for a wind from the disk through the term $`v_z`$, a Galactic wind as an analogy to the Solar wind, driven by Supernova explosions, individual massive star winds, and HII regions.

$`\rho {\displaystyle \frac{}{t}}𝐯={\displaystyle \frac{\sigma _T𝐅}{c}}{\displaystyle \frac{Z}{m_i}}{\displaystyle \frac{\rho }{1+Zm_e/m_i}}+`$ (19)
$`+\rho \nu \mathrm{\Delta }𝐯+\rho (\nu /3+\zeta )\mathrm{grad}\mathrm{div}𝐯+`$ (20)
$`2\rho \omega _𝐏\times 𝐯+𝐓_𝐢+\sigma _{i,e}e𝐄+{\displaystyle \frac{1}{c}}𝐣\times 𝐁P\rho \mathrm{\Phi }`$ (21)
plus non-linear terms

$`𝐓_𝐢=𝐯_𝐢m_i(\mathrm{div}(n_i𝐯_𝐢)S_i)`$ (22)
$`𝐯_𝐞m_e(\mathrm{div}(n_e𝐯_𝐞)S_e)+`$ (23)
$`+𝐯\mathrm{div}(\rho 𝐯)𝐯(m_iS_i+m_eS_e)`$ (24)
$`m_in_i(𝐯_𝐢)𝐯_𝐢`$ (25)
$`m_en_e(𝐯_𝐞)𝐯_𝐞`$ (26)
Momentum exchanges $`𝐏_{\mathrm{𝐢𝐞}}+𝐏_{\mathrm{𝐞𝐢}}=0`$ balance out, and the total pressure sum is $`P=P_i+P_e`$.

### 3.2 Electric current and momentum flow

The momentum exchange term is $`𝐏_{\mathrm{𝐞𝐢}}=en_e\eta 𝐣`$ with $`\eta `$ the resistivity. Subtracting the two main initial equations for ions and electrons gives

$`𝐓_𝐞=Z𝐯_𝐢(\mathrm{div}(n_i𝐯_𝐢)S_i))+`$ (27)
$`+𝐯_𝐞(\mathrm{div}(n_e𝐯_𝐞)S_e)+`$ (28)
$`Zn_i(𝐯_𝐢)𝐯_𝐢+n_e(𝐯_𝐞)𝐯_𝐞`$ (29)
$`{\displaystyle \frac{1}{e}}{\displaystyle \frac{}{t}}𝐣=𝐓_𝐞+{\displaystyle \frac{1}{m_e}}P_e{\displaystyle \frac{Z}{m_i}}P_i`$ (30)
$`2{\displaystyle \frac{1}{e}}\omega _𝐏\times 𝐣{\displaystyle \frac{\sigma _T𝐅}{m_ec}}{\displaystyle \frac{Z}{m_i}}{\displaystyle \frac{\rho }{1+Zm_e/m_i}}`$ (31)
$`\sigma _{i,e}\mathrm{\Phi }+{\displaystyle \frac{1}{e}}\rho (\nu /3+\zeta )\mathrm{grad}\mathrm{div}{\displaystyle \frac{𝐣}{\rho }}+{\displaystyle \frac{1}{e}}\rho \nu \mathrm{\Delta }{\displaystyle \frac{𝐣}{\rho }}`$ (32)
$`{\displaystyle \frac{\sigma _{i,e}}{m_e}}(1{\displaystyle \frac{Zm_e}{m_i}})e𝐄+{\displaystyle \frac{eZ\rho }{m_em_i}}(𝐄\eta 𝐣+{\displaystyle \frac{𝐯}{c}}\times 𝐁+`$ (33)
$`+{\displaystyle \frac{m_i}{Z\rho ec}}({\displaystyle \frac{Zm_e}{m_i}}1)𝐣\times 𝐁)`$ (34)
and with a properly modified term $`𝐓_{𝐞,}`$:

$`{\displaystyle \frac{m_em_i}{e^2Z\rho }}{\displaystyle \frac{}{t}}𝐣={\displaystyle \frac{m_em_i}{eZ\rho }}𝐓_{𝐞,}+𝐣{\displaystyle \frac{S_e}{n_e}}{\displaystyle \frac{m_em_i}{e^2Z\rho }}+`$ (35)
$`+{\displaystyle \frac{m_i}{eZ\rho }}(P_e{\displaystyle \frac{Zm_e}{m_i}}P_i)`$ (36)
$`{\displaystyle \frac{\sigma _T𝐅}{ec}}{\displaystyle \frac{1}{1+Zm_e/m_i}}+`$ (37)
$`+{\displaystyle \frac{m_em_i}{e^2Z}}(\nu /3+\zeta )\mathrm{grad}\mathrm{div}{\displaystyle \frac{𝐣}{\rho }}+{\displaystyle \frac{m_em_i}{e^2Z}}\nu \mathrm{\Delta }{\displaystyle \frac{𝐣}{\rho }}`$ (38)
$`{\displaystyle \frac{m_i}{Z\rho }}\sigma _{i,e}(1{\displaystyle \frac{Zm_e}{m_i}})𝐄\eta 𝐣+𝐄+{\displaystyle \frac{𝐯}{c}}\times 𝐁+`$ (39)
$`{\displaystyle \frac{m_em_i}{eZ\rho }}\sigma _{i,e}\mathrm{\Phi }+{\displaystyle \frac{m_i}{Z\rho ec}}({\displaystyle \frac{Zm_e}{m_i}}1)𝐣\times 𝐁`$ (40)
In deriving these equations we have made various simplifying assumptions in order to keep them tractable, as we will discuss below; we have not insisted to keep all terms involving the lack of charge neutrality.

### 3.3 Number estimates

From the saturation argument of  as applied to Galaxy, and its interstellar medium (hot phase, detected by ROSAT, see , ASCA observations , RXTE measurements, ) we obtain estimates for the various currents and scales.

* The $`\mathrm{curl}`$-operator scale is $`10^{23.2}`$ cm, from $`2\pi r`$, where $`r`$ is the radial distance from us to the center of the Galaxy, about 8 kpc. This we use to infer from Maxwell’s equation the average electric current in the Galaxy, mentioned below.
* The z-gradient, the $``$-operator, here has a scale of 2 kpc = $`10^{21.8}`$ cm, from the scale of the hot gas, and the scaleheight of the magneto-ionic disk is 1.8 kpc , basically the same.
* The resistivity according to Spitzer, , is of order $`10^{17}`$ sec, but highly uncertain; the uncertainty is many powers of ten. This resistivity corresponds to the limit that binary collisions dominate the exchange of momentum, but if the Larmor frequency were to dominate, then the resistivity could be very much higher, possibly by more than ten orders of magnitude.
* The hot gas has density is $`\mathrm{3\hspace{0.17em}10}^3`$ cm-3, and its temperature is $`\mathrm{4\hspace{0.17em}10}^6`$ K , with a total gas pressure of $`10^{11.4}`$ dyn/cm2.
* The total magnetic field is approximately $`10^{5.2}`$ Gauss, and so the total magnetic field pressure is $`10^{11.8}`$ dyn/cm2; we have therefore approximately $`P_{CR}=P_B`$, $`P_{CR}+P_B=10^{11.5}`$ dyn/cm2.
* The electric current is $`j=\mathrm{\hspace{0.17em}10}^{18.9}\mathrm{cgs}\mathrm{units}`$. For the observed symmetry of the magnetic field there must be a current. The existence of this current is the key to our approach.
* These numbers are uncertain, probably to within somewhat less than a factor of 2 (see E. Berkhuijsen in Beck et al. ). A much higher uncertainty derives from the assumption, that everything is basically smooth.

We make the implicit assumption here, that all relevant time scales are of order $`\mathrm{3\hspace{0.17em}10}^7`$ years, or longer; we also assume that the topology of the medium, permeated by magnetic fields, hot and cold gas, and cosmic rays is basically smooth, with no critical different topologies, such as shocks. We basically assume at this stage, that shocks exist, but do not drive any relevant part of the system.

### 3.4 Symmetries

In the following we will use the symmetry with respect to the symmetry plane $`z=0`$ as a first guiding principle:

We introduce the following language:

* We will call a quantity, vector component or scalar “even$`(z,z^n)`$” in the case, that the quantity can be expanded into a power series in $`z`$, that contains only even terms, and starts with $`z^n`$, where $`n`$ is an even number. We will use the expression even$`(z,1)`$ in the case that the quantity is a constant in its first term, so $`n=0`$.
* We will call a quantity, vector component or scalar “odd$`(z,z^n)`$” in the case, that the quantity can be expanded into a power series in $`z`$, that contains only odd terms, and starts with $`z^n`$, where $`n`$ is an odd number.
* We will have to check whether any argument forces us to use a mixed expresssion - but we have not found one yet.
* We will use the language (even, even, odd), when the three vector components in cylindrical coordinates $`r,\varphi ,z`$ are in that sequence even, even, and odd, or, e.g., (odd, odd, even) in a different case.

Looking at Maxwell’s equations and the data, we find then that the main symmetries of different quantities are as follows:

* $`𝐁`$ is (even, even, odd).
* $`𝐣`$ is (odd, odd, even).
* In the observer frame $`𝐄`$ is also (odd, odd, even), as $`𝐣`$. $`𝐄`$ derives from a Lorentz transformation of a shock frame $`𝐁`$-field, and so is not an independent source term.
* The expression $`𝐯\times 𝐁`$ is (odd, odd, even).
* The expression $`𝐣\times 𝐁`$ is (even, even, odd).
* The expression $`P`$ is (even, even, odd).

There is an alternate pattern, in which the magnetic field shows (odd, odd, even), and the electric current shows (even, even, odd) symmetry. It cannot be excluded, that some galaxies have this symmetry. In the following we will pursue the discussion of the main symmetry.

What we find then in the case of the equations for mass flow and electric current flow is that the following is true:

* In the equation for mass flow we have only terms that are either independent of $`𝐣`$, or $`𝐁`$, or quadratic. All these terms have the symmetry (even, even, odd). There are some linear cross-terms, which should cancel each other out; they cannot act as source terms due to a failure to match symmetries.
* In the equation for the electric current the entire equation splits up into terms, that are either (odd, odd, even), or (even, even, odd). Those terms which are linear in either $`𝐣`$, $`𝐄`$ or $`𝐁`$ are all (odd, odd, even). Those terms which are either a source term, or quadratic, are all (even, even, odd).
* Interestingly, the terms involving the pressure gradient, and also the Hall term, which are of similar numerical magnitude, do have the same symmetry as well. They do “belong together”.

The first conclusion is then, with the equations for mass flow, and the electric current as written here, the equations are independent of the (+,-) symmetry with respect to the sign of $`𝐁`$ and $`𝐣`$. Replacing both $`𝐁`$ and also $`𝐣`$ with the same numbers, but with all signs changed to the opposite, does not change the equation. We conclude, that the physics which determines the sign is not yet adequately described by these equations.

### 3.5 Battery effect and dynamo action

Now two arguments can be made:

#### 3.5.1 The dynamo process

Considering the overall equation for driving a current, and taking the curl, and dropping all irrelevant terms, one obtains the basic equation for a discussion of the dynamo process, allowing for Coriolis forces and turbulence in a rotating system. The large number of original references were already given in the introduction.

#### 3.5.2 Battery effect

Assuming no magnetic field to exist a priori, one obtains the Battery effect, . Again taking the $`\mathrm{curl}`$ of the entire equation for the electric current, we have one important term

$$\mathrm{curl}\frac{m_i}{eZ\rho }P_e(\frac{m_i}{eZ\rho })\times P_e$$
(41)
In an accretion disk, where density and temperature are usually governed by very different microphysics, this term is certainly non-zero, and drives an electric current. This the original effect pointed out by the early authors . This is a slow process, although in its limit it gives the reasonable estimates for magnetic fields, see . The key argument here is that an electric current drives the generation of the magnetic field; this electric current is driven by the rotation of the system.

### 3.6 Scaling the hierarchy

Using the observations, and Maxwell’s equations we can discern the numerical relevance of all the terms in the two main equations, for mass flow and for the electric current; the numbers have been given above already. In this we assume at first that there are no discontinuities, to obtain a first guide line of what the observations are telling us.

Three levels of numerical strength exist in the system:

The level structure exists already in the equations connecting momentum flow and electric current to electron and ion velocities. There are two terms also at a lower level: the momentum flow of the electric current is lower than the real mass momentum flow by $`10^{14.4}`$ and $`10^{17.7}`$, respectively. As a consequence this structuring also exists in the fuller equation above, the equation describing the driving of the electric current.

* $`𝐄`$ and $`\frac{𝐯}{c}\times 𝐁`$ of order $`10^8`$ cgs units, using our Galaxy as the prime example.
* Pressure gradients and the Hall term are about $`10^{15}`$ to $`10^{18}`$ times smaller. The radiation driving term is about a factor of 100 weaker again.
* In $`𝐓_𝐞`$ the $`(𝐣)𝐯`$-terms, the time dependent term, the ionization term (assuming an ionization time scale of $`\mathrm{3\hspace{0.17em}10}^7`$ years as for the time dependence), the viscosity terms are all again about $`10^{15}`$ to $`10^{18}`$ times smaller. The resistivity term is nominally a factor of about $`10^4`$ larger, but very badly understood.

We introduce here a way to “scale the hierarchy”, by noting that there may be a large scale shock wave pattern in the Galaxy:

In a spiral galaxy there is a spiral pattern of young massive star formation, and an associated spiral pattern of supernova explosions of just those stars. These supernova explosions produce nearly spherical shockwaves, that propagate outwards, and can reach quite large distances in the tenuous medium. It is easy to show that these shockwaves can overlap in a way that there is always a coherent large scale connected net of shock wave surfaces, locally highly time-dependent, but globally stationary. Locally these shock waves are very time-dependent, and quite fast, obviously faster then the speed of sound in the hot medium, but the global pattern is slow, with a speed much less than the speed of sound in the hot medium. So a comparison with the waves on a beach may be appropiate, very strong waves, as in the Bretagne on the north-western shore of France, but limited by the beach itself and the tides, running up the beach as breakers with a speed that children match only with difficulty.

Let us check on the time scales: In a galaxy like ours we have about one supernova every 100 years, and maybe 1 supernova every $`10^3`$ years in the appropiate mass range, those corresponding to zero age main sequence mass above about 15 to 20 solar masses. This means that a supernova expands into the tenuous hot interstellar medium to a radius of about 50 pc within $`\mathrm{2.5\hspace{0.17em}10}^4`$ years, and the shock then still has a velocity of 500 km/s, a small multiple of the speed of sound in the hot medium. Taking then the notion, that the younger supernovae occur all along the line of the spiral arms, we have about $`10^6`$ yr-1 and pc-1. If each supernova reaches out to 100, pc, within $`\mathrm{2.5\hspace{0.17em}10}^4`$ years, then we have at this frequency about 2 supernovae occurring per 100 pc and per this time, so just a sufficient number.

We need to inquire about the “front” and the “back” of this geometrical set of exploding stars: There is a minimum age for stars to explode as supernovae, but there is a long tail of larger ages, and so one can expect that there is a front pattern, but no back-pattern. At the backside of the spiral pattern there is no corresponding shock; everything peters out.

This concept is similar to the concept for the young supernova shock itself, as in .

In such a shock an electric field exists, driven by the fact that ions and electrons have a different mass. Then the overall global spiral shock drives a current system, a sheet current, and the current is very strong inside the sheet, but weak elsewhere. We note that a sheet current produces a homogeneous magnetic field, with a different sign on opposing sides. The current may be driven by the $`𝐄\times 𝐁`$ drift inside the shock sheet; also, the gradient in pressure and/or temperature in this shock can drive a current. Once one sign is established, it may prevail due to the overall structure winning by emphasizing the common vector components, i.e. those geometrically aligned with the overall shock structure. Since an $`𝐄\times 𝐁`$ drift pushes positive and negative charges in the same direction, only inertia can produce a real electric current; inertia is only important if the velocities reach close to thermal velocities (or relativistic, unlikely here). The requirement to transport angular momentum transport implies that the sign permeates everything in the same way, since any region, where the sign changes would violate the strict direction of angular momentum transport.

Does this solve the problem?

It goes part of the way towards solving the problem by making all three levels equivalent, since the second level was small just by the ratio of the Larmor radius of thermal hot ions to the overall scale, and a shock wave pushes the strength of these terms up by exactly this factor. With a strong sheet-like electric current, stronger by about $`10^{14}`$, but geometrically limited to the shock sheet, we have then equivalence between the first and the second level. And again the third level was smaller than the second by again the same factor. The third level involves a spatial derivative and the current, and so is doubly pushed up in level. Therefore, only in the shock sheet and its immediate environment do we have equivalence of all three levels, previously disjunct. With a shock wave all three levels are becoming equivalent in numerical strength; however, the shock does not change the symmetries.

#### 3.6.1 The current sheet

This means that we have to consider whether such a current can actually be carried by the sheet: the thickness of the sheet is about a postshock thermal Larmor radius $`r_g`$, and so about $`10^9`$ cm, stronger than an average current by the ratio of length scales, $`\pi r`$ vs. $`r_g`$. So the electric current in the sheet is about $`10^5`$ cgs units. This means that the drift velocity is required to be about thermal, so all ions and electrons participate, with a velocity difference about equal to the thermal velocity. This thermal velocity here is the preshock velocity, so the drift can be slow compared to the postshock velocity; however, in order to get a current, we require that electron and ion velocities are quite different. We emphasize again, that this very strong electric current is only in the shock-sheet, and in order to provide the magnetic field observed, has to go vertically.

#### 3.6.2 Angular momentum transport

In this model of expanding supernova remnants angular momentum is mixed, and exchanged, which means that it is effectively transported. The transport is with the azimuthal velocity difference across a supernova remnant, and with the diameter of the supernova remnant. And as shown by Lynden-Bell & Pringle, , and by Duschl et al., , this leads to an angular momentum transport that scales with some fraction of the radius $`r`$, and the azimuthal velocity $`v_\varphi `$, with a factor of order 0.01. This picture developed here does just that. So this way of using large scale shock waves, locally highly non-stationary but globally stationary shock waves, does transport angular momentum. Whether this is already at the level implied by the observations is not clear; the numbers suggest that this mechanism is not sufficient. We need an additional enhancement.

### 3.7 What determines the sign?

Now we have found that even with all the details we do not get a sign prevailing in the detailed description of the flow; this is exactly the same conclusion that Krause & Beck came to in their original work on the sign of magnetic fields in disk galaxies.

#### 3.7.1 The role of the central black hole

Here we wish to explore one way to obtain a specific sign: In the case that a powerfully radiating accretion disk surrounds a central massive black hole, then we have a strongly ionizing ultraviolet photon field. This photon field will ionize and actually charge up a skin in the surrounding matter. Carrying this skin of net charge around produces a circular electric current, with the flow direction exactly the same as the rotation of the material, presumably with the same sense of rotation as the central black hole. This will then produce in turn a magnetic field with a specific sign, aligned with the rotation axis of the black hole, and pointing in the direction of the spin vector. Related questions have been explored in full relativistic treatment by .

This then may determine a definitive sign overall, because it should infect the surrounding matter; a wave of a specific sign will then propagate outwards from this environment. At this point we do not specify how the sign of this inner magnetic field can be turned into a sign of the electric current further out. But once a specific sign is established, it will sustain itself, within the frame work of the mechanism outlined above.

From this picture there are some predictions and speculations: a) Galaxies without a central black hole will not show a specific sign in their magnetic field configuration. b) Galaxies which did not have an accretion event in a very long time, may lose their memory of the original sign. Memory may also be lost, at least for a while, in a merger event with another galaxy, especially if the merger leads to a spin-flip of the central black hole. c) In galaxies at high redshift we may be able to observe the propagation of the wave of a specific sign.

However, this mechanism does not really prefer the observed symmetry over another one, in which the azimuthal magnetic field changes sign in midplane, so a symmetry of the magnetic field which is (odd, odd, even). This latter symmetry is the one that is implicit in any symmetry coming from the central black hole and a charge distribution carried around in corotation. The effect of the radiation of the black hole prefers the alternate symmetry introduced above.

#### 3.7.2 Transport processes?

So we have to look at another option, really a program of investigation:

Considering how the basic equation was derived from the Boltzmann equation, we realize that the stress tensor, the heat flux tensor, and the transport of angular momentum all are dumped into some convenient functions, and simple dependencies that hide the micro-physics. It appears clear from the derivation (see, e.g. Spitzer’s book on plasma physics, ) that we require a vector, and the angular momentum flux, that is probably non-linear in $`𝐣`$ and/or $`𝐁`$, and has the symmetry (even, even, odd) in the expression for the mass flux, while a corresponding vector is required, which has the symmetry (odd, odd, even) in the equation for the electric current. From the idea, that this is Nature’s message, it is clear that the resistivity has to approach the limit of using the gyrofrequency instead of the collision frequency, and that a non-linear behaviour is required, that leads from both the main symmetry and the alternate symmetry only to the main symmetry. Just geometrically, the vector expression $`𝐁\times \mathrm{curl}𝐣`$ would do just that; this expression leads both the main symmetry and the alternate symmetry into driving an electric current with the main symmetry - but the relevance for Nature is at this stage mere speculation. A physical picture for this microphysics has to be left to a later discussion. What may appear at the end, is that the requirement that angular momentum flow is outwards drives the sign of the symmetries, and that we have a highly non-linear expression, using shock structures in a disk to scale the hierarchies. Shocks may appear if the state of the gas is too extreme in its failure to transport angular momentum, and magnetic fields may correspondingly appear. It can be expected that the expressions and coefficients are not independent, just as in the theory of the thermodynamics of irreversible processes.

### 3.8 Where magnetic fields may come from

Therefore we arrive at following speculative picture of the origin of magnetic fields: Stars inject a fairly strong magnetic field, which is highly chaotic. This field is injected with extreme chaos, but from the ram pressure condition for the winds with an Alfvénic Machnumber of order 3, , and so probably with only a factor of order 10 below equipartition. It is actually open to question, whether more is actually required in an inhomogeneous picture of the interstellar medium. The sheet current in the global shock pattern keeps restablishing the strength and the symmetries of the overall regular field, while all the stellar activity produces all the small scale perturbances.

The sheet current involves a very large scale shock, locally unsteady, but globally steady. This is reminiscent of the early discussions in spiral shock waves, . We are reinventing a large scale shock, highly unsteady locally, with a very high maximum local shock velocity, but a rather small pattern velocity.

The global shock proposed here is intimately connected to spiral structure, and thus to the overall flow of angular momentum outwards in any accretion disk. One might speculate whether such a spiral shock could also occur in accretion disks without star formation, and whether in that case it could again drive a sheet current, and thus provide a power source for overall magnetic fields. Could it be that magnetic field generation is connected to the accretion phenomenon for ionized disk as soon as they show a spiral shock structure?

One may ask whether such a shock pattern should not be visible? And indeed this is an interesting question worthwile to pursue with X-ray spectroscopy and high spatial resolution data. X-rays may provide a crucial test for this speculative proposal.

## 4 Conclusions - Future

Magnetic fields are present in the universe, usually strong and relevant.

Where they come from is uncertain, and an attempt outlined here is still incomplete, but we hope, we will be able to figure out using just stars, and the physics of the ionized interstellar medium, how magnetic fields may be generated.

We have outlined a program of investigating the Boltzmann equation to resolve the connections between heat flow, angular momentum transport, and driving currents in an accretion disk. Such an investigation may finally yield the sign of the prevailing magnetic field, and if so, also elucidate the physics of angular momentum transport and magnetic field generation in disks at the same time.

What we have proposed here is to suggest that the overlapping shockwaves from supernova remnants produce a locally unsteady, but globally steady shock structure, which is rather similar in its pattern with the early attempts to find a global spiral shock; however, through its local unsteadiness it is really quite different in its physical properties. This global shock structure drives a sheet current through its electric field and strong gradients in pressure and temperature, which in turn produces by its overall symmetry a global magnetic field, and so rejuvenates continuously the overall topology of the Galactic magnetic magnetic field.

Cosmic magnetic fields are key to the transport of charged particles throughout the Universe, at any energy.

## 5 Acknowledgements

P.L. Biermann would like to acknowledge the hospitality during his prolonged stay at the University of Paris VII, in the spring 2002, offered by his colleagues Norma Sanchez and Hector de Vega, as well as their hospitality at many meetings in Erice, Paris and Palermo. Work with PLB is being supported through the AUGER theory and membership grant 05 CU1ERA/3 through DESY/BMBF (Germany); further support for the work with PLB comes from the DFG, DAAD, Humboldt Foundation (all Germany), grant 2000/06695-0 from FAPESP (Brasil) through G. Medina-Tanco, a grant from KOSEF (Korea) through H. Kang and D. Ryu, a grant from ARC (Australia) through R.J. Protheroe, and European INTAS/ Erasmus/ Sokrates/ Phare grants with partners V. Berezinsky, M.V. Rusu, and V. Ureche. The authors would like to thank R. Beck, E. Berkhuijsen, N. Ikhsanov, H. Kang, P.P. Kronberg, H. Lee, G. Medina-Tanco, D. Mitra, and D. Ryu for many discussions of these questions.


text:
# Gravitational Collapse of Null Radiation and a String fluid

## I Introduction

The end state of a collapsing star has attracted much attention in its own right and especially so since the statement of the Cosmic Censorship Conjecture (CCC) (See joshi  for a detailed discussion of the CCC.). While the end state of a sufficiently massive collapsing star is a singularity he , there is no guarantee that such a singularity should be hidden as is required by the CCC penrose . Indeed, a number of papers have been written describing various scenarios in which naked singularities form (See joshi2  for a recent review.). Such examples clearly violate the CCC. Future formulations of the CCC will have to be carefully re-considered in light of these examples.

Here we present another spherically symmetric spacetime which admits a naked singularity. This spacetime is not asymptotically flat and describes the collapse of null radiation with a string fluid first introduced by Glass and Krisch gk . A full physical description of the model is contained in gk1 . We merely quote the relevant information for our analysis of the singularities.

However, we do make one observation: In their analysis Glass and Krisch gk ; gk1  state that the matter content of their model is a radiation-string two fluid. The matter content is taken to be a string fluid together with Vaidya null (pure) radiation. We proceed in a similar manner for the main part of our paper. However, in §V we consider the special case of Ricci flat spacetimes. Here, we do not assume that the string satisfies a diffusion equation. The combination of these two characteristics (Ricci flatness and non-diffusing strings) now allows the matter content to be treated as a different radiation-string two fluid (or, perhaps, a three fluid), where, in addition to the Vaidya null (pure) radiation, we also have incoherent radiation (as the radiation fluid equation of state $`p=\rho /3`$ is now satisfied) together with the string fluid.

This note is organised as follows: In §II we introduce the metric of Glass and Krisch gk  and present the main equations to solve. We make a few comments on these equations and provide some special solutions. These solutions are used in §III to study the nature of the collapse. We show that one solution results in a black hole while the other could result in a naked singularity. The strength of the naked singularity is examined in §IV. During the analysis, we present the general forms of the Kretschmann and Ricci scalars. We investigate the collapse of Ricci flat spacetimes in §V and show that the presence of a string fluid can lead to the occurrence of locally naked singularities for radiation fluids. We compare our results to a similar metric that has recently appeared in the literature in §VI.

## II The Null Radiation–String Fluid Model

Glass and Krisch gk  start with a spherically symmetric metric in retarded time co-ordinate, $`u`$, of the form

$$\mathrm{d}s^2=A\mathrm{d}u^2+2\mathrm{d}u\mathrm{d}rr^2(\mathrm{d}\vartheta ^2+\mathrm{sin}^2\vartheta \mathrm{d}\phi ^2),$$
(1)
where

$$A=1\frac{2m(u,r)}{r}.$$
(2)
Note that they chose the mass function to be a function of both retarded time and distance along the outgoing null geodesics. As a result, they were able to extend the Vaidya vaidya  radiating metric to include a string fluid.

Remark: A spherically symmetric line element in advanced time $`u`$ is usually taken to be of the form

$$\mathrm{d}s^2=A\mathrm{d}u^2+2\mathrm{d}u\mathrm{d}r+r^2(\mathrm{d}\vartheta ^2+\mathrm{sin}^2\vartheta \mathrm{d}\phi ^2),$$
(3)
Thus, the difference between (1) and (3) is not only a signature change, but also a change of sign of $`u`$. We will comment further on this later.

Glass and Krisch gk  introduced the unit vectors $`\widehat{v}^\mu ,\widehat{r}^\mu ,\widehat{\vartheta }^\mu `$ and $`\widehat{\phi }^\mu `$ such that

$$g_{\mu \nu }=\widehat{v}_\mu \widehat{v}_\nu \widehat{r}_\mu \widehat{r}_\nu \widehat{\vartheta }_\mu \widehat{\vartheta }_\nu \widehat{\phi }_\mu \widehat{\phi }_\nu ,$$
(4)
where the unit vectors have their usual definition from (1). The Einstein tensor then took on the form

$$G_{\mu \nu }=\frac{2\dot{m}}{r^2}l_\mu l_\nu \frac{2m^{}}{r^2}(\widehat{v}_\mu \widehat{v}_\nu \widehat{r}_\mu \widehat{r}_\nu )+\frac{m^{\prime \prime }}{r}(\widehat{\vartheta }_\mu \widehat{\vartheta }_\nu +\widehat{\phi }_\mu \widehat{\phi }_\nu ),$$
(5)
where  and  denote differentiation with respect to $`r`$ and $`u`$ respectively and $`l_\mu `$ is a principal null geodesic vector, from which the nature of the matter content is easily seen.

Using the Einstein field equations they identified the matter portion of the Einstein tensor as a string fluid:

$$T_{\mu \nu }=\psi l_\mu l_\nu +\rho \widehat{v}_\mu \widehat{v}_\nu +p_r\widehat{r}_\mu \widehat{r}_\nu +p_{}(\widehat{\vartheta }_\mu \widehat{\vartheta }_\nu +\widehat{\phi }_\mu \widehat{\phi }_\nu )$$
(6)
and thus obtained

$`4\pi \psi `$ $`=`$ $`{\displaystyle \frac{\dot{m}}{r^2}}`$ (7)
$`4\pi \rho `$ $`=`$ $`4\pi p_r={\displaystyle \frac{m^{}}{r^2}}`$ (8)
$`8\pi p_{}`$ $`=`$ $`{\displaystyle \frac{m^{\prime \prime }}{r}},`$ (9)
where $`\psi `$ is the energy density of the Vaidya radiation, $`\rho `$ is the string fluid energy density and $`p_r`$ and $`p_{}`$ the string fluid thermodynamic radial and tangential pressures respectively.

In addition to the field equations (7)–(9) they also assumed that the density $`\rho `$ satisfied the diffusion equation

$$\dot{\rho }=\frac{𝒟}{r^2}\frac{}{r}\left(r^2\frac{\rho }{r}\right),$$
(10)
where $`𝒟`$ is the positive coefficient of self-diffusion. Manipulating equations (7)–(10) they concluded that, once a solution to the linear pde (10) was obtained, the mass function could be found by integrating (8) and

$$\dot{m}=4\pi 𝒟r^2\frac{\rho }{r}.$$
(11)
They presented a few solutions of (10) and determined the corresponding expressions for $`m`$.

Equation (11) is a result of the integrability conditions taking (7), (8) and (10) into account. However, it would seem that an additive function of $`u`$ was omitted in (11) which should be written as (See Appendix A.)

$$\dot{m}=4\pi 𝒟r^2\frac{\rho }{r}+k(u).$$
(12)
However, in our calculations, we will still use (11) and merely add on an arbitrary $`k(u)`$.

We observe that (10) can be transformed to the heat equation with a nonconstant diffusion coefficient if we set $`R=1/r`$. Thus all the known solutions to the heat equation can be applied in this setting.

It is easy to see that, when $`\rho =R(r)U(u)`$ we obtain

$$\rho =\mathrm{e}^{u\lambda }\left(\frac{C_1}{\mathrm{e}^{r\sqrt{\frac{\lambda }{𝒟}}}r}+\frac{C_2\mathrm{e}^{r\sqrt{\frac{\lambda }{𝒟}}}}{2r\sqrt{\frac{\lambda }{𝒟}}}\right),$$
(13)
where $`C_1`$ and $`C_2`$ are constants and $`\lambda `$ is the separation constant. This leads to the following mass function:

$`m(u,r)`$ $`=`$ $`{\displaystyle \frac{4𝒟\mathrm{e}^{u\lambda r\sqrt{\frac{\lambda }{𝒟}}}\pi C_1}{\lambda }}{\displaystyle \frac{4𝒟\mathrm{e}^{u\lambda r\sqrt{\frac{\lambda }{𝒟}}}\pi r\sqrt{\frac{\lambda }{𝒟}}C_1}{\lambda }}`$ (14)
$`+{\displaystyle \frac{2𝒟\mathrm{e}^{u\lambda +r\sqrt{\frac{\lambda }{𝒟}}}\pi rC_2}{\lambda }}{\displaystyle \frac{2𝒟^2\mathrm{e}^{u\lambda +r\sqrt{\frac{\lambda }{𝒟}}}\pi \sqrt{\frac{\lambda }{𝒟}}C_2}{\lambda ^2}}.`$
If we require that $`\rho =R(r)+U(u)`$ then

$$\rho =C_0+\lambda u+\frac{C_1}{r}+\frac{\lambda r^2}{6𝒟}$$
(15)
with corresponding mass function

$$m(u,r)=\frac{2\pi r^5\lambda }{15𝒟}4𝒟\pi uC_1+\frac{2\pi r^2(2r(u\lambda +C_0)+3C_1)}{3}+C_2.$$
(16)
If we let

$$M(u)=4𝒟\pi C_1u$$
(17)
be the Vaidya mass and note that the $`4\pi C_0r^3/3`$ term makes a contribution similar to the cosmological constant, we can define

$$S(u,r)=\frac{2\pi r^5\lambda }{15𝒟}+\frac{2\pi r^2(2\lambda ru+3C_1)}{3}+C_2$$
(18)
as a mass contribution from the string fluid. We note that both $`S(u)`$ and $`M(u)`$ have a term involving $`C_1`$. Thus, in this model, it is not possible to separate the Vaidya mass from the string mass. This is a result of the relationships (8) and (11) which are a consequence of (10).

We note that, if we require the string to behave like a perfect fluid, then we must have $`p_r=p_{}=p`$ and so $`p=\rho `$ remains as the equation of state for the string fluid. Using equations (8) and (9), this means that

$$rm^{\prime \prime }2m^{}=0$$
(19)
whence

$$m=F(u)r^3+G(u).$$
(20)
If we further impose (10), we have that $`F(u)=F_0`$, a constant and the resulting part of the mass function, ie. $`F_0r^3`$ is a cosmological constant term – the metric is Vaidya-de Sitter (See wm  for a discussion of the CCC in this context.). In the resulting discussion, we will assume that the fluid is anisotropic.

## III Nature of the collapse

The geodesic equations (taking the null condition $`K^aK_a`$ into account, where $`K^a=\mathrm{d}x^a/\mathrm{d}k`$ is the tangent vector of a geodesic) for the metric (1) are

$`{\displaystyle \frac{\mathrm{d}^2u}{\mathrm{d}k^2}}+{\displaystyle \frac{mm^{}r}{r^2}}\left({\displaystyle \frac{\mathrm{d}u}{\mathrm{d}k}}\right)^2`$ $`=`$ $`0`$ (21)
$`{\displaystyle \frac{\mathrm{d}^2r}{\mathrm{d}k^2}}{\displaystyle \frac{\dot{m}}{r}}\left({\displaystyle \frac{\mathrm{d}u}{\mathrm{d}k}}\right)^2`$ $`=`$ $`0.`$ (22)
To investigate the structure of the collapse we need to consider the radial null geodesics defined by $`\mathrm{d}s^2=0`$ taking $`\dot{\varphi }=\dot{\theta }=0`$ into account. For (1) the radial null geodesics must satisfy the null condition

$$\frac{\mathrm{d}u}{\mathrm{d}r}=\frac{2}{1\frac{2m(u,r)}{r}}.$$
(23)
However, we recall that the co-ordinate $`u`$ is a retarded time co-ordinate. As remarked earlier, in most examples in the literature, an advanced time co-ordinate is used. In order to compare our results with those usually found in the literature we need to set $`uu`$ in (23) and obtain

$$\frac{\mathrm{d}u}{\mathrm{d}r}=\frac{2}{1\frac{2m(u,r)}{r}}.$$
(24)
For appropriate forms of $`m`$, this equation will have a singularity at $`r0,u0`$. In order to classify the radial and non-radial outgoing non-spacelike geodesics terminating at this singularity in the past, we need to consider the limiting values of $`X=u/r`$ along a singular geodesic as the singularity is approached joshi . Thus, for the geodetic tangent to exist uniquely at this point we must have that

$$X_0=\underset{u0r0}{lim}\frac{u}{r}=\underset{u0r0}{lim}\frac{\mathrm{d}u}{\mathrm{d}r}=\underset{u0r0}{lim}\frac{2}{1\frac{2m(u,r)}{r}}.$$
(25)
Thus the crucial aspect of metric (1) that dictates the nature of the collapse is, as expected, the mass function $`m`$.

We note that $`\rho =\rho _0`$, a constant, is a solution of (10). The corresponding mass function is $`m(u,r)=m_0+4\pi \rho _0r^3/3`$. For this form of $`m`$, we have

$$X_0=\underset{u0r0}{lim}\frac{2}{1\frac{2m_0+8\pi \rho _0r^3/3}{r}}=\frac{2}{1lim_{r0}\frac{2m_0}{r}}=0.$$
(26)
Since there is no real, positive $`X_0`$, there is no non-spacelike geodesic emanating from the singularity and so the singularity is not visible to any observer. However, this is not a realistic model as $`m`$ is independent of time and so the model is static. The reason we comment on this form of $`m`$ is to show that any constant part of $`m`$ will result in a black hole. As a result in our further analysis of the structure of the collapse we will always set $`m_0=0`$.

Let us consider solution (14). It is clear from that solution that, as $`r,u0`$, we have $`mk`$, a constant and so, from (26) that a black hole forms.

In the case of (16) (with $`C_2=0`$) we observe that (25) reduces to

$$X_0=\frac{2}{18\pi 𝒟C_1X_0}$$
(27)
from which we obtain

$$X_0=a_\pm =\frac{1\pm \sqrt{164C_1\pi 𝒟}}{16C_1\pi 𝒟}.$$
(28)
If we define $`\delta =8C_1\pi 𝒟`$ we can rewrite (28) as

$$X_0=a_\pm =\frac{1\pm \sqrt{18\delta }}{2\delta }.$$
(29)
In order for real positive solutions in (29) to exist we must have that

$$\delta \frac{1}{8}.$$
(30)
This is exactly the same form of $`X_0`$ (and restriction) obtained by Dwivedi and Joshi dj  and Wagh and Maharaj wm . That the similarity arises should not be of any surprise. It is clear that, if one is dealing with a mass function that is of the form $`\lambda u+F(u,r)`$ where $`F(u,r)/r`$ is polynomial in $`r`$ and/or $`u`$ (of degree at least one), then the only contribution to the condition for the existence of a time-like geodesic emanating from the singularity is from $`\lambda u`$, the Vaidya mass. Indeed, this is the reason that asymptotic flatness is not essential for the existence of a locally naked singularity, a conclusion first obtained in wm .

## IV Strength of the singularity

The main importance of determining the strength of the singularity is due to the fact that the CCC does not need to rule out the possibility of the occurrence of weak naked singularities nolan . This arises as one may continue the geometry of a weak naked singularity through the singularity to make it geodesically complete es ; ori  (See also nolan .).

In order to determine the strength of the singularity we utilise the ideas of Nolan nolan  as explained in gd . Thus, for a weak singularity, we require

$$\frac{\mathrm{d}r}{\mathrm{d}k}d_0rd_0k.$$
(31)
Using the definition

$$X_0=\underset{u0r0}{lim}\frac{u}{r}=\underset{u0r0}{lim}\frac{\mathrm{d}u}{\mathrm{d}r}$$
(32)
we have, asymptotically,

$$\frac{\mathrm{d}u}{\mathrm{d}k}d_0X_0ud_0X_0k.$$
(33)
The geodesic equation (21) then becomes

$$\frac{d^2u}{dk^2}\frac{2X_0^2\pi d_0}{15k𝒟}(4\lambda d_0^4k^4+20d_0^3k^3dX_0\lambda +20d_0^2k^2𝒟C_0+15d_0k𝒟C_1+30d^2X_0C_1)$$
(34)
which is of $`O(k)^1`$ and so is inconsistent with (33). Thus the singularity is gravitationally strong in the sense of Tipler tip .

The Kretschmann scalar $`K=R_{abcd}R^{abcd}`$ for (1) is given by

$$K=\frac{4}{r^6}((m^{\prime \prime })^24m^{\prime \prime }m^{}r^3+4mm^{\prime \prime }r^2+8(m^{})^2r^216mm^{}r+12m^2).$$
(35)
If we utilise the mass function (16), (35) evaluates to

$`K`$ $`=`$ $`{\displaystyle \frac{64\pi ^2}{225𝒟^2r^6}}(1200r^6𝒟^2u\lambda C_0+600r^6𝒟^2u^2\lambda ^2+900r^5𝒟^2C_0C_1+300r^8\lambda ^2𝒟u+300r^8\lambda 𝒟C_0+150r^7\lambda 𝒟C_1`$
$`+2700𝒟^4u^2C_1^2+720r^5\lambda 𝒟^2uC_1+600r^6𝒟^2C_0^2+450r^4𝒟^2C_1^2+53r^{10}\lambda ^2)`$
which clearly diverges at the naked singularity. Hence the singularity is a scalar polynomial singularity.

The Ricci scalar for (1) is given by

$$R=\frac{2}{r^2}(m^{\prime \prime }r+2m^{})$$
(36)
which, in the case of (16) evaluates to

$$R=\frac{8\pi }{𝒟r}(r^3\lambda +4𝒟ru\lambda +4𝒟rC_0+3𝒟C_1)$$
(37)
and also diverges at the naked singularity.

## V Ricci Flat Spacetimes

It is clear from (36) that none of the solutions presented above cover the case of Ricci flat spacetimes. If we relax the requirement that the string diffuses (ie. we do not impose (10)) we can now consider the case of $`R=0`$. From (36) this means that the mass function must take on the form

$$m(u,r)=\overline{M}(u)+\frac{\overline{S}(u)}{r},$$
(38)
where we can again interpret $`\overline{M}(u)`$ as the Vaidya mass and $`\overline{S}(u)`$ as the mass contribution from the matter content. Forcing the spacetime to be Ricci flat requires the matter content to obey the (incoherent) radiation fluid equation of state, ie.

$$p\frac{p_r+2p_{}}{3}=\frac{p_r}{3}=\frac{\rho }{3}$$
(39)
from (8) and (9). Thus if we begin with a metric of the form (1) (together with an energy–momentum tensor of the form (6)) and require it to be Ricci flat, then, as expected, the fluid must satisfy (39). These spacetimes must be asymptotically flat via (38).

To proceed further, we must assume some form for $`\overline{M}(u)`$ and $`\overline{S}(u)`$. Given the interpretation above, it makes sense to assume that $`\overline{M}(u)=\alpha u`$. In order for $`\overline{S}(u)`$ to contribute to the existence (or otherwise) of a naked singularity, it must be of $`O(u)^2`$. Thus we set $`\overline{S}(u)=\beta u^2`$. The null condition (24) now becomes

$$\frac{\mathrm{d}u}{\mathrm{d}r}=\frac{2}{1\frac{2\alpha u+2\beta u^2/r}{r}}.$$
(40)
This equation has a singularity at $`r0,u0`$. The condition for a geodetic tangent to exist uniquely at this point reduces to

$$X_0=\underset{u0r0}{lim}\frac{u}{r}=\underset{u0r0}{lim}\frac{\mathrm{d}u}{\mathrm{d}r}=\underset{u0r0}{lim}\frac{2}{1\frac{2\alpha u+2\beta u^2/r}{r}}=\frac{2}{12\alpha X_02\beta X_0^2}$$
(41)
This cubic equation in $`X_0`$ has (at least) one real positive root provided that, for $`\alpha >0`$, $`\beta <0`$. This can be seen from the numerical results mma  given in Table 1. Thus, in the case of Ricci flat spacetimes, we also have the possibility of a locally naked singularity.

As above, we investigate the strength of the naked singularity. Taking (31) and (33) into account, (21) (with (38) and the assumed forms of $`\overline{M}(u)`$ and $`\overline{S}(u)`$) becomes

$$\frac{\mathrm{d}^2u}{\mathrm{d}k^2}\frac{2\beta X_0\alpha }{d_0k}$$
(42)
which is inconsistent with (33) and so this singularity is also gravitationally strong.

For Ricci flat spacetimes, the Kretschmann scalar (35) is given by

$$K=\frac{16}{r^8}(14\overline{S}(u)^2+12\overline{S}(u)\overline{M}(u)r+3\overline{M}(u)^2r^2)$$
(43)
which, for the assumed forms of $`\overline{M}(u)`$ and $`\overline{S}(u)`$, evaluates to

$$K=\frac{16u^2}{r^8}(14\beta ^2u^212\beta \alpha ur+3\alpha ^2r^2).$$
(44)
Clearly (44) diverges at the naked singularity. As expected, we have a scalar polynomial singularity.

## VI Discussion

We have presented a scenario for the gravitational collapse of a string fluid together with Vaidya null radiation in which the end state becomes a naked singularity. We have also investigated the collapse of Ricci flat spacetimes. Such collapse can also give rise to naked singularities. Both scenarios provided naked singularities which were gravitationally strong and, as expected, were scalar polynomial singularities. These are two more in the long list of examples that proliferate the literature and should serve to tighten future formulations of the CCC.

We comment briefly on our results and those of Ghosh and Dadhich gd . The metrics (apart from a signature change and change of sign of time co-ordinate) are similar. As a result, the Einstein field equations also reduce to (79) (though their $`p`$ represents isotropic pressure). However, in their model, Ghosh and Dadhich gd  interpret the matter content as strange quark matter. Thus, while $`\psi `$ is still the radiation density, their $`p`$ and $`\rho `$ refer to strange quark matter thermodynamic pressure and energy density respectively. As a result of that interpretation, they impose an equation of state relating $`p`$ to $`\rho `$:

$$p=\frac{1}{n}(\rho 4B),$$
(45)
where $`B`$ is the bag constant. This differential equation can be treated as an ode and the general solution was presented in gd .

Here, following Glass and Krisch gk , the matter content has been interpreted as a string fluid which diffuses like point particle diffusion where the number density, $`n`$, diffuses from higher numbers to lower according to (10) (Note we have used the fact that $`\rho =M_sn`$, $`M_s`$ a constant, to obtain (10).). The string equation of state is always satisfied as a result of the energy-momentum tensor (6). As a result, the solutions presented here (and in gk ) will not coincide with those found in gd .

That other interpretations of the metric (1) together with (6) are possible is evident from the case of Ricci flat collapse investigated in §V. Indeed, this case is of considerable importance as the matter content satisfies both the string equation of state $`p_r=\rho `$ and the (incoherent) radiation fluid equation of state $`p=\rho /3`$ (see (39)). Thus the matter content is a two-fluid, in addition to the Vaidya null (pure) radiation. We have shown that a naked singularity can arise in this scenario. In cjjs  it was shown that a naked singularity only forms if the ratio of pressure to density was less than $`1/3`$. Here, we have shown that a naked singularity can exist outside this range, ie. for a ratio of $`1/3`$. We conclude that it is the presence of the string matter that contributes to the occurrence of the locally naked singularity.

As observed in §II, an arbitrary function of $`u`$, $`k(u)`$ say, can be added to the mass function $`m`$. This additive function of $`u`$ can be seen as a further contribution from the string mass $`S`$. However, unless it is linear in $`u`$, it would not contribute to the condition (30) and even then, would only serve to redefine $`\delta `$.

We note that the usual manner in which the existence of a naked singularity is presented is to show that the singularity forms before the event horizon forms. Interestingly, it has been shown that, for some models, an apparent horizon will not form at all bd ; h1 ; k1  (See also waghetal .).

###### Acknowledgements.

We thank Naresh Dadhich, Roy Maartens and Sunil Maharaj for useful discussions and helpful suggestions. KSG thanks the University of Natal for continuing support. Many of the tensor calculations have been performed using GRTensor 1.79 grtens . We thank the referees for their constructive criticisms that have helped improve this paper.

## Appendix A Integrability conditions for $`m(u,r)`$

Here we consider the derivation of (11). Differentiating (7) and (8) with respect to $`r`$ and $`u`$ respectively results in

$`\dot{m}^{}`$ $`=`$ $`4\pi (r^2\psi )^{}`$ (46)
$`\dot{m}^{}`$ $`=`$ $`4\pi r^2\dot{\rho }.`$ (47)
Thus

$$(r^2\psi )^{}=r^2\dot{\rho }.$$
(48)
However, (10) yields

$$r^2\dot{\rho }=𝒟\frac{}{r}\left(r^2\frac{\rho }{r}\right).$$
(49)
Thus, from (46), (48) and (49) we have

$$(r^2\psi )^{}=\frac{\dot{m}^{}}{4\pi }=𝒟\frac{}{r}\left(r^2\frac{\rho }{r}\right)$$
(50)
whence

$$\dot{m}^{}=4\pi 𝒟\frac{}{r}\left(r^2\frac{\rho }{r}\right).$$
(51)
Glass and Krisch gk  presumably integrated this equation to obtain (11). However, since this is a pde, the integration should yield (12).


text:
# Comment on “Including Systematic Uncertainties in Confidence Interval Construction for Poisson Statistics”

## Abstract

The incorporation of systematic uncertainties into confidence interval calculations has been addressed recently in a paper by Conrad et al. (Physical Review D 67 (2003) 012002). In their work, systematic uncertainities in detector efficiencies and background flux predictions were incorporated following the hybrid frequentist-Bayesian prescription of Cousins and Highland, but using the likelihood ratio ordering of Feldman and Cousins in order to produce “unified” confidence intervals. In general, the resulting intervals behaved as one would intuitively expect, i.e. increased with increasing uncertainties. However, it was noted that for numbers of observed events less than or of order of the expected background, the intervals could sometimes behave in a completely counter-intuitive fashion – being seen to initially *decrease* in the face of increasing uncertainties, but only for the case of increasing *signal* efficiency uncertainty. In this comment, we show that the problematic behaviour is due to integration over the signal efficiency uncertainty while maximising the best fit alternative hypothesis likelihood. If the alternative hypothesis likelihood is determined by unconditionally maximising with respect to both the unknown signal and signal efficiency uncertainty, the limits display the correct intuitive behaviour.

In their recent paper, Conrad et al.bchd  incorporated uncertainties on signal and background efficiencies into limit calculations by applying the standard classical likelihood ratio techniquesoa , recently popularised by Feldman and Cousinsfc98 , to the hybrid Bayesian-frequentist method of Cousins and Highlandch . In this method, the fixed-but-unknown signal strength parameter $`\mu _s`$ is treated in a frequentist fashion, but the unknown experimental efficiencies ($`ϵ_s`$ and $`ϵ_b`$) are incorporated by integrating over prior probability distributions $`P(ϵ_s)`$ and $`P(ϵ_b)`$. This simplifies the confidence interval calculation by reducing the dimensionality of the problem down to only one unknown variable. The resulting intervals showed the correct intuitive behaviour in most cases, i.e. the confidence intervals were seen to increase as the uncertainties in the efficiencies were increased. However, it was noted that some counter-intuitive behaviour could occur for cases where the observed number of events $`n_0`$ was less than or of order of the expected background $`\mu _b`$. In these cases, the limits sometimes initially became more restrictive as the uncertainties were increased, a behaviour which anyone would agree was undesirable. This behaviour was also noted in an earlier paper by Giuntigiunti . In this comment, we show that the noted undesirable behaviour is due to the choice of likelihood ratio test implemented by Conrad et al., and show how a choice more consistent with the generalised likelihood ratio test in the presence of nuisance parameters eliminates this behaviour, leading to intervals with acceptable behaviour as uncertainties are increased. The key to the correction is the choice of treatment of the uncertainties in the likelihood ratio denominator; Conrad et al. chose to integrate over uncertainties in both the numerator and denominator, which leads to the counter intuitive behavior. We show here how integrating in the numerator, but maximising with respect to the uncertainties in the denominator leads to the correct intuitive behaviour.

We consider specifically the problematic case from Conrad et al., the determination of a confidence limit on an unknown Poisson signal strength $`\mu _s`$ in the presence of a precisely known background $`\mu _b`$ but where the signal efficiency $`ϵ_s`$ is also unknown. Before considering the approach of Conrad et al., we note that a completely frequentist treatment would involve the construction of a confidence region in the $`\mu _sϵ_s`$ plane, where the acceptability of each possible outcome $`n`$ under a null hypothesis $`\{\mu _s,ϵ_s\}`$ is determined by a likelihood ratio testsoa ; fc98 , where the likelihood ratio $`R`$

$$R=\frac{(n\mu _s,\mu _b,ϵ_s)}{(n\widehat{\mu }_s,\mu _b,\widehat{ϵ}_s)}$$
(1)
compares the likelihood under the null hypothesis to the likelihood of the alternative hypothesis $`\{\widehat{\mu }_s,\widehat{ϵ}_s\}`$ that best describes the possible observation $`n`$, i.e. $`\widehat{\mu }_s`$ and $`\widehat{ϵ}_s`$ are the values that unconditionally maximise the likelihood of observing $`n`$ events. Values of $`n`$ where $`(n\mu _s,\mu _b,ϵ_s)`$ is close to the maximum possible value $`(n\widehat{\mu }_s,\mu _b,\widehat{ϵ}_s)`$ are most acceptable and are therefore first included into acceptance intervals during the confidence interval construction. To reduce the dimensionality of the problem, the suggested proceduresoa  is to conditionally maximise the numerator with respect to the nuisance parameter $`ϵ_s`$ yielding a likelihood ratio

$$R=\frac{(n\mu _s,\mu _b,\widehat{\widehat{ϵ}}_s)}{(n\widehat{\mu }_s,\mu _b,\widehat{ϵ}_s)}$$
(2)
The Cousins-Highland approach is a hybrid frequentist-Bayesian method where the incorporation of a systematic uncertainty in the efficiency $`ϵ_s`$ proceeds by integrating over a prior probability function $`P(ϵ_s\widehat{ϵ}_s,\sigma _{ϵ_s})`$, describing our knowledge of the nuisance parameter $`ϵ_s`$, but where the unknown signal strength $`\mu _s`$ is still treated in a classical frequentist fashion. This yields a likelihood function

$`(n\mu _s,\mu _b,P(ϵ_s\widehat{ϵ}_s,\sigma _{ϵ_s}))={\displaystyle _{ϵ_s^{}}}(n\mu _s,\mu _b,ϵ_{s}^{}{}_{}{}^{})P(ϵ_s^{}\widehat{ϵ}_s,\sigma _{ϵ_s})\mathrm{d}ϵ_s^{}`$ (3)
which is used to construct classical confidence intervals in the fixed but unknown signal strength $`\mu _s`$. We still need an ordering principle to decide which values of $`n`$ are to be included in the acceptance interval for each null hypothesis $`\mu _s`$. Conrad et al. used a Feldman-Cousins style likelihood ratio ordering principle (equation 1 or 2) to choose which values of $`n`$ fall into the acceptance interval for a given $`\mu _s`$. This unifies the treatment of upper limits and two-sided confidence regions. However, it is their treatment of the uncertainties in the determination of the best alternative hypothesis in the likelihood ratio test which leads to the counter-intuitive behaviour of the confidence intervals. Conrad et al. used a likelihood ratio of the form

$$R=\frac{_{ϵ_s}(n\mu _s,\mu _b,ϵ_s^{})P(ϵ_s^{}\widehat{ϵ}_s,\sigma _{ϵ_s})dϵ_s^{}}{_{ϵ_s}(n\widehat{\mu }_s,\mu _b,ϵ_s^{})P(ϵ_s^{}\widehat{ϵ}_s,\sigma _{ϵ_s})dϵ_s^{}}$$
(4)
where the uncertainties are integrated over in both the numerator and denominator. Integration over the uncertainties has the desired effect on the numerator likelihood where values $`n`$ near the mean $`\widehat{ϵ}_s\mu _s+\mu _b`$ become less probable in favour of higher and lower values of $`n`$ as $`\sigma _{ϵ_s}`$ increases. Integration over uncertainties while finding $`\widehat{\mu }_s`$ still results in $`\widehat{\mu }_s=0.0`$ for $`n<\mu _b`$. However, for $`n>\mu _b`$, the best fit value $`\widehat{\mu }_s`$ for a given $`n`$ as well as the probability of this best fit value used in the denominator decreases as $`\sigma _{ϵ_s}`$ increases. The net effect is a shift in the peak of the likelihood ratio distribution to higher $`n`$, that then leads to the problematic behaviour in the confidence intervals by sometimes shifting the acceptance region $`\{n_{\mathrm{𝑙𝑜}}(\mu _s\sigma _{ϵ_s}=0.0),n_{\mathrm{ℎ𝑖}}(\mu _s\sigma _{ϵ_s}=0.0)\}`$ for zero uncertainties to higher values of $`n`$ where

$$n_{\mathrm{𝑙𝑜}}(\mu _s\sigma _{ϵ_s}>0.0)>n_{\mathrm{𝑙𝑜}}(\mu _s\sigma _{ϵ_s}=0.0)$$
(5)
as the uncertainty is included. Since the upper limit $`\mu _s^{\mathrm{lim}}`$ for a given observation $`n_0`$ is the highest value of $`\mu _s`$ for which $`n_0`$ is still in the acceptance interval, shifting the acceptance interval for values of $`\mu `$ at or below the limit $`\mu _s^{\mathrm{lim}}`$ will lead to a lowering of the limit. An example of this effect is shown in table 1, where we examine the specific problematic case from Conrad et al.. The 90% confidence level upper limits for observations of $`n_0=2,4`$ and 6 events on a precisely known background of $`\mu _b=6.0`$ in the presence of various percentage Gaussian uncertainties in the signal efficiency are shown, using the likelihood ratio of Conrad et al. (equation 4). We note that our limits for the Conrad et al. ordering are not exactly as given in their paper as we here test the null hypothesis in increments of 0.01, rather than 0.05. Nonetheless, the limit for the case of $`n_0=2`$ is seen to decrease as the efficiency error increases from 0 to 40% $`(1.561.45)`$, before finally increasing. For the $`n_0=4`$ and 6 cases, the limits decrease as $`\sigma _ϵ`$ goes from 0 to 10%, but then increase.

We can correct this behaviour by constructing a hybrid frequentist-Bayesian likelihood ratio test, where the numerator retains the integration over the uncertainties but where the alternative hypothesis is found by unconditionally maximising the denominator as in the pure frequentist tests in equations 1 and 2. This yields the following likelihood ratio statistic

$$R=\frac{_{ϵ_s^{}}(n\mu _s,\mu _b,ϵ_{s}^{}{}_{}{}^{})P(ϵ_s^{}\widehat{ϵ}_s,\sigma _{ϵ_s})dϵ_s^{}}{(n\widehat{\mu }_s,\mu _b,\widehat{ϵ}_s)}$$
(6)
where the inclusion into the acceptance interval of given $`n`$’s is determined by comparison of their likelihood after integration over uncertainties to the likelihood of the best fit $`\widehat{\mu }_s`$ given no uncertainty. This way, the likelihoods given any value of the uncertainty $`\sigma _{ϵ_s}`$ are always compared to the simple alternative hypothesis $`\{\widehat{\mu }_s,\widehat{ϵ}_s\}`$ which best describes the observation $`n`$, rather to one that changes with $`\sigma _{ϵ_s}`$.

Table 1 shows that the limits using the ordering from equation 6 are seen to initially remain the same (due to overcoverage from the discrete nature of the Poisson distribution) but then increase with increasing uncertainty. We also note that although the Conrad et al. limits do finally increase as the uncertainties grow bigger, they are always smaller than the limits found using the likelihood ratio in equation 6.

###### Acknowledgements.

I wish to thank Jan Conrad for useful discussions and Katherine Rawlins for helpful comments on the manuscript. This work was supported by the National Science Foundation under constract number OPP-9980474.