# Zero-field optical magnetometer based on spin-alignment

A. Meraki ID<EMAIL_ADDRESS>School of Physics and Astronomy,
University of Nottingham, University Park, Nottingham, NG7 2RD, UK  L. Elson
ID School of Physics and Astronomy, University of Nottingham, University Park,
Nottingham, NG7 2RD, UK  N. Ho ID School of Physics and Astronomy, University
of Nottingham, University Park, Nottingham, NG7 2RD, UK  A. Akbar ID School of
Physics and Astronomy, University of Nottingham, University Park, Nottingham,
NG7 2RD, UK  M. Koźbiał ID Centre for Quantum Optical Technologies, Centre of
New Technologies, University of Warsaw, Banacha 2c, 02-097 Warszawa, Poland
J. Kołodyński ID Centre for Quantum Optical Technologies, Centre of New
Technologies, University of Warsaw, Banacha 2c, 02-097 Warszawa, Poland  K.
Jensen ID<EMAIL_ADDRESS>School of Physics and Astronomy, University of
Nottingham, University Park, Nottingham, NG7 2RD, UK

###### Abstract

Optically-pumped magnetometers are an important instrument for imaging
biological magnetic signals without the need for cryogenic cooling. These
magnetometers are presently available in the commercial market and utilize the
principles of atomic alignment or orientation, enabling remarkable sensitivity
and precision in the measurement of magnetic fields. This research focuses on
utilizing a spin-aligned atomic ensemble for magnetometry at zero-field. A
novel approach is introduced, which involves evaluating how the linear
polarization of light rotates as it passes through the atomic vapor to null
the magnetic field. Analytical expressions are derived for the resulting spin
alignment and photodetection signals. Experimental results are provided,
demonstrating good agreement with the theoretical predictions. The sensitivity
and bandwidth of the magnetometer are characterized based on the detected
polarization rotation signal. Lastly, the practical utility of the
magnetometer for medical applications is demonstrated by successfully
detecting a synthetic cardiac signal.

##  I Introduction

Over the past decade, significant progress has been made in enhancing the
sensitivity of optically pumped magnetometers (OPMs) used in extremely low
magnetic field environments [1, 2, 3]. These advancements have surpassed the
sensitivity of superconducting quantum interference devices (SQUIDs) [4] while
eliminating the requirement for cryogenic cooling. Current commercially
available OPMs [5, 6, 7] are commonly operated in the spin-exchange
relaxation-free (SERF) regime, enabling measurements of one or more components
of the magnetic field near zero-field [8, 9]. OPMs have demonstrated their
exceptional ability to precisely measure small magnetic fields in a variety of
applications, some examples include biomedical imaging [10, 11, 12, 13, 14,
15], metal detection [16, 17] and material characterizations [18]. The
miniaturization and commercialization of magnetometers is highly desired,
especially for integrating them into densely packed arrays for medical imaging
applications [19, 20].

OPMs typically utilize optical pumping techniques with either circularly [3,
16] or linearly polarized light [21, 22, 23, 24]. Circularly polarized light
is frequently used in most zero-field OPM configurations to orient the total
spin of the atomic ensemble towards a specific direction [25]. Alternatively,
one can use linearly polarized light, which enables the preparation of atomic
alignment for spins larger than 1/2 [26, 27]. Linearly polarized light in the
alignment configuration can be used for pumping and probing atoms using a
single laser beam. This can be utilized in a zero-field setting [28, 29, 30].
Zeroing the magnetic field in an OPM is crucial to minimize any unintended
influences, such as ambient magnetic fields or equipment-induced disturbances.
Effectively nulling the magnetic field allows researchers to set up a reliable
reference point, enhancing the instrument’s sensitivity and precision for a
diverse range of applications.

In our work, we investigate how to use a spin-aligned atomic ensemble for
zero-field magnetometry. We introduce a novel approach for nulling the
magnetic field which involves evaluating how the linear polarization of the
light is rotated as it is transmitted through the atomic vapour. We derive
analytical expressions for the generated spin-alignment and the photodetection
signals. A room temperature paraffin coated cell containing caesium (Cs) is
used to study how the polarization rotation signal changes with small magnetic
fields near zero. By analyzing the polarization of light as it passes through
the vapor cell under different magnetic field conditions, we determine the
appropriate direction to adjust the compensation magnetic field and whether to
increase or decrease its magnitude. Having presented the method for nulling
the field, we then characterize the magnetometer’s sensitivity and bandwidth
from the detected polarization rotation signal. Our zero-field magnetometer is
modulation-free, as opposed to previous work which typically employs
oscillating magnetic fields for modulating the atomic response and where the
absorption of light is detected [29, 28]. We finally demonstrate that our
magnetometer is useful for medical applications by detecting a synthetic
cardiac signal.

##  II Theory

###  II.1 Equation of motion for alignment multipoles

Our experiments are carried out using an ensemble of room-temperature caesium
atoms in a vapour cell. Caesium has two ground states with hyperfine quantum
numbers $F=3$ and $F=4$ and the first excited state has hyperfine quantum
numbers $F^{\prime}=3$ and $F^{\prime}=4$ [31]. We use linearly polarized
light on-resonant with the $6S_{1/2}$ $\rightarrow$ $6P_{1/2}$,
$F=4\rightarrow F^{\prime}{=}3$ D1 transition at a wavelength of 895 nm to
optically pump and probe the atoms. One can describe a single atom in the
ground state with a particular hyperfine quantum number $F$ with a density
matrix $\rho$ with size $\left(2F+1\right)\times\left(2F+1\right)$ using the
basis $\left\\{m=-F,m=-F+1,...,m=F\right\\}$, where $m$ is the projection of
the total angular momentum along the $\bm{e}_{z}$ direction. Alternatively,
one can express the quantum state as

$\rho=\sum_{k=0}^{2F}\sum_{q=-k}^{k}m_{k,q}T_{q}^{k},$ (1)

where $m_{k,q}$ are multipoles with rank $k=0,1,...,2F$ and $q=-k,...,k$ and
$T_{q}^{k}$ are irreducible spherical tensor operators [21, 32, 28]. Note that
in total there are $\sum_{k=0}^{2F}\left(2k+1\right)=\left(2F+1\right)^{2}$
multipole components which is equal to the number of density matrix elements.
For a normalized state $m_{0,0}=1/\sqrt{2F+1}$. The $k=1$ multipole describes
how the spin is oriented, and they can be organised in a vector

$\mathbf{m}_{1}=\left(m_{1,-1},m_{1,0},m_{1,1}\right)^{T}.$ (2)

Here “T” denotes the transpose such that $\mathbf{m}_{1}$ is a 3-dimensional
column vector. The $k=2$ multipoles describe how the spin is aligned along
some axes, and they can be organised in the 5-dimensional column vector

$\mathbf{m}_{2}=\left(m_{2,-2},m_{2,-1},m_{2,0},m_{2,1},m_{2,2}\right)^{T}.$
(3)

It can be convenient to describe the quantum state using the multipoles
instead of the density matrix, due to the properties of the multipoles under
rotations. For example, atoms optically pumped with circular polarized light
can be well described by their orientation $\mathbf{m}_{1}$, while atoms
optically pumped with linear polarized light can be well described by their
alignment $\mathbf{m}_{2}$.

(a)

(b)

(c)

Figure 1: (a) Sketch of the polarization rotation of light transmitted through
an atomic ensemble. (b) Toy model of an atom optically pumped with linearly
polarized light. (c) Experimental setup. Red paths indicate optical pathways
and grey paths represent the routes for the electrical signals. PM Fiber,
polarization maintaining fiber; PBS, polarizing beam splitter; $\lambda/2$,
half-wave plates; BPD, balanced photodetector.

For simplicity, we will use a simpler toy-model to describe the evolution of
the single atom quantum state. Consider an atom with a $F=1$ ground state that
has three sub-levels $\left\\{m=-1,m=0,m=1\right\\}$ and a single excited
state with $F^{\prime}=0,m=0$, where the quantization axis is chosen along the
$\bm{e}_{z}$ direction. The $F=1$ ground state can be described either by a
$3\times 3$ density matrix $\rho$ or by the multipoles $m_{k,q}$ with
$k=0,1,2$ and $q=-k,...,k$.

Let us assume that the beam of light, linearly polarized along the
$\bm{e}_{z}$ direction, is propagating through an ensemble of atoms, as
illustrated in Fig. 1(a). In the absence of any magnetic field, the light will
optically pump and polarize each atom. With perfect optical pumping, where the
populations in the $m=1$ and $m=-1$ sub-levels are both equal to $1/2$, as
illustrated in Fig. 1(b), we have $m_{0,0}=1/\sqrt{3}$,
$\mathbf{m}_{1}=\left(0,0,0\right)^{T}$ and
$\mathbf{m}_{2}=\left(0,0,1/\sqrt{6},0,0\right)^{T}$. We see that the
orientation multipole is zero and that the quantum state is described with
just the alignment multipole. If the optical pumping is not perfect, then the
generated alignment will be reduced. In any case, the generated alignment can
be described by the expression

$\mathbf{m}_{2}^{\mathrm{eq}}=\left(0,0,m_{2,0}^{\mathrm{eq}},0,0\right)^{T}.$
(4)

We now assume the atom is placed in a magnetic field

$\bm{B}=B_{x}\bm{e}_{x}+B_{y}\bm{e}_{y}+B_{z}\bm{e}_{z}.$ (5)

Due to the properties of spherical tensors under rotations, it is possible to
rewrite the Liouville equation for the density matrix as a matrix equation for
the multipoles. Including decay and repopulation terms, we find a matrix
equation written compactly as

$\displaystyle\dot{\mathbf{m}}_{2}=$
$\displaystyle\left(-i\omega_{x}J_{x}^{(2)}-i\omega_{y}J_{y}^{(2)}-i\omega_{z}J_{z}^{(2)}\right)\mathbf{m}_{2}$
$\displaystyle-\Gamma\mathds{1}_{5}+\Gamma\mathbf{m}_{2}^{\mathrm{eq}}.$ (6)

Here $\omega_{x}=\gamma B_{x}$, $\omega_{y}=\gamma B_{y}$, and
$\omega_{z}=\gamma B_{z}$, where $\gamma$ is the gyromagnetic ratio of the
atom. It is possible to extend our model to include multiple relaxation rates
[21, 23]. However, for simplicity, we assume just a single ground-state
relaxation rate $\Gamma$. The matrices $J_{x}^{(2)}$, $J_{y}^{(2)}$, and
$J_{z}^{(2)}$ are the $k=2$ representation of the angular momentum operators
and $\mathds{1}_{5}$ is the $5\times 5$ identity matrix. Eq. (II.1) has been
explicitly written out in Appendix A in the Supplemental Material. We can find
the steady-state solution to Eq (II.1) by setting the time-derivative to zero
$\dot{\mathbf{m}}_{2}^{\mathrm{ss}}=0$ and solve the system of linear
equations. We find

$\displaystyle\begin{pmatrix}m_{2,-2}^{\mathrm{ss}}\\\
m_{2,-1}^{\mathrm{ss}}\\\ m_{2,0}^{\mathrm{ss}}\\\ m_{2,1}^{\mathrm{ss}}\\\
m_{2,2}^{\mathrm{ss}}\end{pmatrix}=$
$\displaystyle\frac{m_{2,0}^{\mathrm{eq}}}{\left(\Gamma^{2}+\omega_{x}^{2}+\omega_{y}^{2}+\omega_{z}^{2}\right)\left[\Gamma^{2}+4\left(\omega_{x}^{2}+\omega_{y}^{2}+\omega_{z}^{2}\right)\right]}$
$\displaystyle\times$
$\displaystyle\begin{pmatrix}-\sqrt{\frac{3}{2}}\left(\omega_{x}+i\omega_{y}\right)^{2}\left(\Gamma^{2}+\omega_{x}^{2}+\omega_{y}^{2}+3i\Gamma\omega_{z}-2\omega_{z}^{2}\right)\\\
-\sqrt{\frac{3}{2}}\left(\omega_{x}+i\omega_{y}\right)\left(i\Gamma+2\omega_{z}\right)\left(\Gamma^{2}+\omega_{x}^{2}+\omega_{y}^{2}+3i\Gamma\omega_{z}-2\omega_{z}^{2}\right)\\\
\Gamma^{4}+\left(\omega_{x}^{2}+\omega_{y}^{2}-2\omega_{z}^{2}\right)^{2}+\Gamma^{2}\left[2\left(\omega_{x}^{2}+\omega_{y}^{2}\right)+5\omega_{z}^{2}\right]\\\
\sqrt{\frac{3}{2}}\left(\omega_{x}-i\omega_{y}\right)\left(-i\Gamma+2\omega_{z}\right)\left(\Gamma^{2}+\omega_{x}^{2}+\omega_{y}^{2}-3i\Gamma\omega_{z}-2\omega_{z}^{2}\right)\\\
-\sqrt{\frac{3}{2}}\left(\omega_{x}-i\omega_{y}\right)^{2}\left(\Gamma^{2}+\omega_{x}^{2}+\omega_{y}^{2}-3i\Gamma\omega_{z}-2\omega_{z}^{2}\right)\end{pmatrix}.$
(7)

Note that $\left(m_{2,-2}^{\mathrm{ss}}\right)^{*}=m_{2,2}^{\mathrm{ss}}$ and
$\left(m_{2,-1}^{\mathrm{ss}}\right)^{*}=-m_{2,1}^{\mathrm{ss}}$, where ∗
denotes the complex conjugate. We also see that
$m_{2,0}^{\mathrm{ss}}=m_{2,0}^{\mathrm{eq}}$ when the field is zero
($\omega_{x}=\omega_{y}=\omega_{z}$).

###  II.2 Magnetometry signals

We can calculate the observable signals from the steady-state multipoles in
Eq. (II.1). We first consider a measurement of how much light is transmitted
through a vapour cell with length $l$. We have $I(l)=I(0)e^{-\kappa l}$, where
$I(0)$ is the light intensity before the vapour cell, $I(l)$ is the light
intensity after the vapour cell, and $\kappa$ is the absorption coefficient.
The transmitted light can be measured with a single photodetector.
Alternatively, we can measure how the linear polarization of the light is
rotated as it traverses the vapour cell. This polarization rotation signal $S$
can be measured using a balanced photodetector. One can calculate formulas for
the absorption $\kappa$ and the polarization rotation signal $S$. Such
formulas depend on how the coordinate system is defined. Above we assumed that
the light is polarized in the $\bm{e}_{z}$ direction and we now furthermore
assume that the light is propagating in the $\bm{e}_{x}$ direction. In this
case, one can show (see Appendix B in the Supplemental Material) that the
absorption and polarization-rotation signals generally obey

$\kappa\propto\left(m_{0,0}-\chi\cdot m_{2,0}\right),$ (8)

and

$S\propto i\left(m_{2,1}+m_{2,-1}\right),$ (9)

respectively, where only the parameter $\chi$ depends on the hyperfine quantum
number $F$ considered.

For instance, in the case of $F=1$ we have that $\chi=\sqrt{2}$ and
$m_{0,0}=1/\sqrt{3}$ (for a normalized state). For an un-polarized state,
where each $m$-sublevel is occupied with $1/3$ probability, all multipoles are
zero except $m_{0,0}$. Therefore, the absorption for an un-polarized state
will be $\kappa\propto m_{0,0}=1/\sqrt{3}$. On the other hand, for a fully
polarized state, the absorption will be
$\kappa\propto\left(m_{0,0}-\sqrt{2}m_{2,0}\right)=\frac{1}{\sqrt{3}}-\sqrt{2}\cdot\frac{1}{\sqrt{6}}=0$
at zero-field. I.e., the absorption is zero, and all light is transmitted.
This makes sense as the atoms are optically pumped into a dark state. For a
non-zero magnetic field, we insert the steady state solution given by Eq.
(II.1) into Eqs. (8) and (9) and obtain

$\kappa\propto
m_{0,0}-\chi\cdot\frac{m_{2,0}^{\mathrm{eq}}\left[\Gamma^{4}+\left(\omega_{x}^{2}+\omega_{y}^{2}-2\omega_{z}^{2}\right)^{2}+\Gamma^{2}\left\\{2\left(\omega_{x}^{2}+\omega_{y}^{2}\right)+5\omega_{z}^{2}\right\\}\right]}{\left(\Gamma^{2}+\omega_{x}^{2}+\omega_{y}^{2}+\omega_{z}^{2}\right)\left[\Gamma^{2}+4\left(\omega_{x}^{2}+\omega_{y}^{2}+\omega_{z}^{2}\right)\right]},$
(10)

and

$S\propto\frac{\sqrt{6}m_{2,0}^{\mathrm{eq}}\left[\Gamma^{3}\omega_{x}-\Gamma^{2}\omega_{y}\omega_{z}+2\omega_{y}\omega_{z}\left(\omega_{x}^{2}+\omega_{y}^{2}-2\omega_{z}^{2}\right)+\Gamma\omega_{x}\left(\omega_{x}^{2}+\omega_{y}^{2}+4\omega_{z}^{2}\right)\right]}{\left(\Gamma^{2}+\omega_{x}^{2}+\omega_{y}^{2}+\omega_{z}^{2}\right)\left[\Gamma^{2}+4\left(\omega_{x}^{2}+\omega_{y}^{2}+\omega_{z}^{2}\right)\right]}.$
(11)

If we assume two components of the magnetic field are zero (e.g. $\omega_{y}$
and $\omega_{z}$) and that the third component is large (e.g.
$\omega_{x}\gg\Gamma$), then Eq. (10) simplifies to

$\kappa\propto m_{0,0}-\chi\cdot m_{2,0}^{\mathrm{eq}}\cdot\frac{1}{4}.$ (12)

This means that the absorption in this case is $1/4$ compared to that of an
un-polarized state.

We finally emphasize that in the theoretical description of our magnetometer
we assumed a toy model atom with a $F=1$ ground state, which obviously is an
approximation for a real cesium atom which have hyperfine structure and two
ground states with $F=3$ and $F=4$. We note that all numbered equations are
still valid for a cesium atom, if we assume that the multipoles correspond to
the density matrix in the $F=4$ sub-space. Although a state with $F=4$ in
general can have multipoles with $k=0,1,...,8$, only the $k=2$ multipoles are
relevant for our work. This is because the light only probes the $k\leq 2$
multipoles (when assuming the light power is low), and because multipoles with
different $k$’s do not couple to each other under rotations (e.g. during the
evolution in a magnetic field). For $F=4$ we have
$\chi=\tfrac{1}{2}\sqrt{\tfrac{11}{7}}$ which can be calculated from Eq. (8)
by noting that the state $\rho=\left(1/2\right)|F=4,m=4\rangle\langle
F=4,m=4|+\left(1/2\right)|F=4,m=-4\rangle\langle F=4,m=-4|$ is a dark state
for light linearly polarized along $\bm{e}_{z}$.

###  II.3 Magnetometry near zero-field

It is instructive to Taylor expand the full formulas for the absorption and
polarization rotation given by Eqs. (10) and (11) to obtain simpler
expressions for the signal. To keep notation short, we define
$x=\omega_{x}/\Gamma$, $y=\omega_{y}/\Gamma$, and $z=\omega_{z}/\Gamma$ which
are proportional to the $B_{x}$, $B_{y}$, and $B_{z}$ components of the
magnetic field, respectively.

We first consider the absorption measurement. For example, if sweeping the
$B_{x}$ magnetic field, and assuming the $B_{y}$ and $B_{z}$ are close to zero
(i.e. $y,z\ll 1$), then we find

$\kappa\propto m_{0,0}-\chi\cdot
m_{2,0}^{\mathrm{eq}}\cdot\left[\frac{1}{1+x^{2}}\right].$ (13)

I.e., the absorption signal has a Lorentzian lineshape as a function of $x$,
and the absorption does not depend on $y$ and $z$ to first order. Similarly,
if sweeping the $B_{y}$ magnetic field, and assuming the $B_{x}$ and $B_{z}$
are close to zero (i.e. $x,z\ll 1$), then we find

$\kappa\propto m_{0,0}-\chi\cdot
m_{2,0}^{\mathrm{eq}}\cdot\left[\frac{1}{1+y^{2}}\right].$ (14)

Finally, if sweeping the $B_{z}$ magnetic field, and assuming the $B_{x}$ and
$B_{y}$ are close to zero (i.e. $x,y\ll 1$), then we find

$\kappa\propto m_{0,0}-\chi\cdot m_{2,0}^{\mathrm{eq}}\cdot 1.$ (15)

In this case the absorption is constant and independent of $x$, $y$ and $z$ to
first order.

We now consider the polarization rotation signal. If sweeping the $B_{x}$
magnetic field, and assuming the $B_{y}$ and $B_{z}$ are close to zero (i.e.
$y,z\ll 1$), then we find to first order in $y$ and $z$ that

$S\propto\sqrt{6}m_{2,0}^{\mathrm{eq}}\cdot\left[\frac{x}{1+4x^{2}}\right].$
(16)

This implies that signal does not depend on $y$ and $z$ to the first order and
it has a dispersive lineshape as a function of $x$. We can also calculate the
slope (i.e. the derivative with respect to $x$) of Eq. (11) to second order in
$y$ and $z$. We find

$\frac{\partial S}{\partial
x}_{|x=0}\propto\sqrt{6}m_{2,0}^{\mathrm{eq}}\cdot\left[\frac{1+y^{2}+4z^{2}}{\left(1+y^{2}+z^{2}\right)\left(1+4y^{2}+4z^{2}\right)}\right],$
(17)

which means that the slope is maximal when $y=z=0$.

On the other hand, if sweeping the $B_{y}$ magnetic field, and assuming the
$B_{x}$ and $B_{z}$ are close to zero (i.e. $x,z\ll 1$), then we find to first
order in $x$ and $z$ that

$S\propto\sqrt{6}m_{2,0}^{\mathrm{eq}}\cdot\left[x\cdot\frac{1}{1+4y^{2}}+z\cdot\frac{y\left(-1+2y^{2}\right)}{\left(1+y^{2}\right)\left(1+4y^{2}\right)}\right].$
(18)

Finally, if sweeping the $B_{z}$ magnetic field, and assuming the $B_{x}$ and
$B_{y}$ are close to zero (i.e. $x,y\ll 1$), then we find to first order in
$x$ and $y$ that

$S\propto\sqrt{6}m_{2,0}^{\mathrm{eq}}\cdot\left[x\cdot\frac{1}{1+z^{2}}-y\cdot\frac{z}{1+z^{2}}\right].$
(19)

From Eqs. (16), (17) (18), and (19) we see that the polarization rotation
signal to first order is a measure of the $B_{x}$ component of the field, i.e.
the component along the light propagation direction.

We also see that it is possible to use the polarization rotation signal given
by Eqs. (18) and (19) as a sensitive method for nulling the magnetic field.
Besides the measured polarization rotation signal, one needs a coil system for
applying fields in the $\bm{e}_{x}$, $\bm{e}_{y}$ and $\bm{e}_{z}$ directions.
For example, as seen from Eq. (18), if sweeping the $B_{y}$ magnetic field,
one can adjust the $B_{x}$ and $B_{z}$ magnetic field components using the
coil system until the signal is zero. When that is the case, those field
components will be nulled $B_{x}=B_{z}=0$. Similarly, as seen from Eq. (19),
if sweeping the $B_{z}$ magnetic field, one can adjust the $B_{x}$ and $B_{y}$
magnetic field components until the signal is zero, and when that is the case,
$B_{x}=B_{y}=0$.

##  III Experimental Setup and Procedure

We implemented the zero-field magnetometry based on spin-alignment to a table-
top, laboratory setup, see schematics in Fig. 1(c). Light with wavelength 895
nm resonant with the $F=4\rightarrow F^{\prime}=3$ Cs D1 transition was
emitted from a fiber pigtailed diode laser (Thorlabs DBR895PN). To adjust the
intensity of light, a half-wave plate was employed in combination with a
polarizer that aligned the polarization along the $\bm{e}_{z}$-axis. The light
power (measured before the vapour cell) was 3 $\mu$W for the data presented in
Figs. 2, 3, 4, and 5, and 10 $\mu$W for the data presented in Figs. 6 and 7.
The linearly $z$-polarized light entered a $\left(5~{}\mathrm{mm}\right)^{3}$
cubic vapor cell which was located in a four-layer mu-metal magnetic shield
(Twinleaf MS-1L). The vapor cell is hand-blown, contains caesium atomic
vapour, and is coated with paraffin. This leads to a long spin-coherence time.
Following that, polarimetry is performed by employing a half-wave plate, a
polarizing beam splitter and a balanced photodetector (Thorlabs PDB210A/M).
The output from the balanced photodetector is subsequently directed through a
2 kHz low-pass filter and recorded with a data acquisition card within the
computer for further analysis. For the absorption measurements (presented in
Sec. IV.2), the transmitted light intensity was measured directly with a
single photodetector (Thorlabs PDA36A2) placed immediately after the vapour
cell outside the magnetic shield. The caesium cell is positioned within three
sets of coils, where each coil set can be used to apply fields in three
($\bm{e}_{x}$, $\bm{e}_{y}$, $\bm{e}_{z}$) directions. With no current applied
to the coils, there is a small residual magnetic field
$\bm{B}^{\mathrm{residual}}$ inside the shield which is on the order of 100
nT. The outer coil set (which comes as standard with the Twinleaf MS-1L
magnetic shield) is used to apply a static field $\bm{B}^{\mathrm{DC}}$ to
compensate the residual magnetic field in order for us to completely null the
total magnetic field. The middle coil set is a three-axis cubic Helmholtz coil
set which can be used to generate a linearly sweeping magnetic field
$\bm{B}^{\mathrm{sweep}}$ which is used in the procedure for nulling the
static field. The inner coil set, which is also a three-axis cubic Helmholtz
coil, is used for generating different signals $\bm{B}^{\mathrm{signal}}$,
e.g. a sinusoidal magnetic field for characterizing the magnetometer response
and sensitivity or a synthetic cardiac signal. The total magnetic field can be
written as

$\bm{B}=\bm{B}^{\mathrm{residual}}+\bm{B}^{\mathrm{DC}}+\bm{B}^{\mathrm{sweep}}+\bm{B}^{\mathrm{signal}}.$
(20)

##  IV Results and Discussions

###  IV.1 Method for nulling the static field

(a)

(b)

(c)

(d)

Figure 2: Sweeping a magnetic field $B_{z}^{\mathrm{sweep}}$ applied along the
$\bm{e}_{z}$ direction near zero-field condition. Solid lines are experimental
data, and dashed lines are fits to Eq. (19). (a) and (b) show signals for
various settings of magnetic field applied along the $\bm{e}_{x}$ and
$\bm{e}_{y}$ directions, respectively. (c) and (d) show the corresponding fit
parameters
$x\propto\left(B_{x}^{\mathrm{residual}}+B_{x}^{\mathrm{DC}}\right)$ and
$y\propto\left(B_{y}^{\mathrm{residual}}+B_{y}^{\mathrm{DC}}\right)$. The DC
fields are $B_{z}^{\mathrm{DC}}=-97$nT, $B_{y}^{\mathrm{DC}}=-44$nT in (a) and
$B_{x}^{\mathrm{DC}}=14$nT in (b). The half width at half maximum of the
resonances in (a) equals $\Gamma/\gamma$, which was fitted for each resonance
in (a) and (b) and was found to be within the range 17 ($\pm$1) nT.

(a)

(b)

(c)

(d)

Figure 3: Sweeping a magnetic field $B_{y}^{\mathrm{sweep}}$ applied along the
$\bm{e}_{y}$ direction near zero-field condition. Solid lines are experimental
data, and dashed lines are fits to Eq. (18). (a) and (b) show signals for
various settings of an static magnetic field applied along the $\bm{e}_{x}$
and $\bm{e}_{z}$ directions, respectively. (c) and (d) show the corresponding
fit parameters
$x\propto\left(B_{x}^{\mathrm{residual}}+B_{x}^{\mathrm{DC}}\right)$ and
$z\propto\left(B_{z}^{\mathrm{residual}}+B_{z}^{\mathrm{DC}}\right)$. The DC
fields are $B_{y}^{\mathrm{DC}}=-44$nT, $B_{z}^{\mathrm{DC}}=-97$nT in (a) and
$B_{x}^{\mathrm{DC}}=14$nT in (b). The parameter $\Gamma/\gamma$ was fitted
for each resonance and was found to be within the range 12 ($\pm$ 1) nT for
the resonances in (a) and 15 ($\pm$1) nT for the resonances in (b).

We now experimentally demonstrate how to utilize the polarization rotation
signal $S$ to null the total static magnetic field
$\bm{B}^{\mathrm{residual}}+\bm{B}^{\mathrm{DC}}$. This involves
systematically sweeping magnetic fields $\bm{B}^{\mathrm{sweep}}$ applied in
the $\bm{e}_{y}$ and $\bm{e}_{z}$ directions, and then changing the applied
static field $\bm{B}^{\mathrm{DC}}$ until a zero signal $S\approx 0$ is
obtained. In that case, the total static field should be
$\bm{B}^{\mathrm{residual}}+\bm{B}^{\mathrm{DC}}\approx 0$ as seen by Eqs.
(18) and (19). The nulling procedure can be done iteratively, and the results
presented here were all recorded in conditions close to zero total static
magnetic field.

Figure 2(a) presents the measured signals with a linearly sweeping applied
magnetic field $B^{\mathrm{sweep}}_{z}$ in the $\bm{e}_{z}$ direction. For
these measurements, the static field component $B_{x}^{\mathrm{DC}}$ was
varied within the range of 12 nT to 16 nT, while the other components were
fixed at $B_{y}^{\mathrm{DC}}=-44$ nT and $B_{z}^{\mathrm{DC}}=-97$ nT. These
settings correspond to a total static magnetic field close to zero. The
observed signals looks like Lorentzian lineshapes as predicted by our model,
see Eq. (19), when $x\neq 0$ and $y=0$. Note that for $B_{x}^{\mathrm{DC}}=14$
nT, the signal transforms into a flat line as predicted by Eq. (19) when $x=0$
and $y=0$. For this setting, the applied static magnetic field precisely
cancels the residual magnetic field leading to a net zero static magnetic
field $\bm{B}^{\mathrm{residual}}+\bm{B}^{\mathrm{DC}}\approx 0$.

Measurements of the signal while sweeping the applied field
$B^{\mathrm{sweep}}_{z}$ in the $\bm{e}_{z}$ direction are also obtained when
varying the applied static field $B^{\mathrm{DC}}_{y}$ component. Those
results are presented in Fig. 2(b). The lineshapes now look dispersive, as
predicted by Eq. (19) with $x=0$ and $y\neq 0$.

The obtained data in Figs. 2(a) and 2(b) are fitted to Eq. (19). Figure 2(c)
shows the parameter
$x\propto\left(B_{x}^{\mathrm{residual}}+B_{x}^{\mathrm{DC}}\right)$ obtained
by fitting the data in Fig. 2(a). As expected, we see a linear relationship
between the applied static field and the total static field. The parameter $y$
(not shown) was also fitted for the data in Fig. 2(a) and was found to be much
smaller than $x$. Figure 2(d) shows the parameter
$y\propto\left(B_{y}^{\mathrm{residual}}+B_{y}^{\mathrm{DC}}\right)$ obtained
for the data in Fig. 2(b). Again we see a linear relationship. The parameter
$x$ (not shown) was also fitted for the data in Fig.2(b) and was found to be
much smaller than $y$.

The procedure is now repeated but with instead sweeping a field
$B_{y}^{\mathrm{sweep}}$ applied in the $\bm{e}_{y}$ direction. Results are
shown in Fig. 3. We see that it is possible to to adjust the
$B_{x}^{\mathrm{DC}}$ and $B_{z}^{\mathrm{DC}}$ magnetic field components
until the measured signal reaches zero. As mentioned, the nulling procedure
can be repeated iteratively, until one obtains a signal very close to zero,
for some specific settings of $B_{x}^{\mathrm{DC}}$, $B_{y}^{\mathrm{DC}}$,
and $B_{z}^{\mathrm{DC}}$. The total static field has then been nulled to
better than a small fraction of a nT, as seen in Figs. 2(c), 2(d), 3(c) and
3(d) which is much smaller than the linewidth of the magnetic resonances, and
therefore sufficient for zero-field magnetometry. A larger range where
$B^{\mathrm{residual}}+B^{\mathrm{DC}}$ is varied from -40 nT to +40 nT can be
seen in Appendix C in the Supplemental Material where the data is fitted to
Eq. (11).

###  IV.2 Zero-field resonances measured by absorption

Having nulled the total static magnetic field, we now measure the zero-field
resonances using absorption, i.e. by measuring the intensity of the
transmitted light using a single detector placed immediately after the vapour
cell. If a magnetic field is applied, it affects the quantum state of the
caesium atoms as calculated in Eq. (II.1) and will lead to less light being
transmitted through the vapour cell. To first order, the atoms are only
sensitive to the magnetic field components $B_{x}$ and $B_{y}$ which are
transverse to the light propagation direction. Figure 4 shows the sequential
sweeping of the three components of the magnetic field using triangular ramps
of ±300 (nT). We see that the signal is constant when sweeping $B_{z}$, and
that the signal has a Lorentzian lineshape when sweeping either $B_{x}$ or
$B_{y}$, as expected from Eqs. (13), (14) and (15). Due to the constant /
Lorentzian lineshapes, a direct measurement of the intensity of the
transmitted light is not useful for determining the magnetic field components,
as several field values, e.g. $B_{x}$ & $-B_{x}$, or $B_{y}$ & $-B_{y}$ give
rise to the same transmitted light intensity, which means that one cannot
uniquely infer a field component from the signal. However, we note that it is
possible to measure all components of the magnetic field using modulation
techniques, where several auxiliary oscillating magnetic fields are applied
[28, 29].

Figure 4: When sweeping the magnetic field from positive to negative values
for $B_{x}$, $B_{y}$, and $B_{z}$, zero-field resonances can be observed in a
direct measurement of the transmitted light intensity. The parameter
$\Gamma/\gamma$ corresponding to the half width at half maximum of the
resonances was fitted and found to be 18.44 ($\pm$ 0.01) nT and 18.48 ($\pm$
0.01) nT when sweeping $B_{x}$ and $B_{y}$, respectively.

(a)

(b)

Figure 5: Sweeping a magnetic field applied along the $\bm{e}_{x}$ direction
near zero-field condition. Solid lines are experimental data, and dashed lines
are fits to Eq. (16). (a) and (b) show signals for various settings of an
magnetic field applied along the $\bm{e}_{y}$ and $\bm{e}_{z}$ directions,
respectively. In (a) $B_{z}^{\mathrm{DC}}=-65$nT and in (b)
$B_{y}^{\mathrm{DC}}=-49$nT. The parameter $\Gamma/\gamma$ was fitted for each
resonance in (a) and (b) and was found to be within the range 11 ($\pm$ 1) nT.

###  IV.3 Zero-field magnetometry measured by polarization rotation

####  IV.3.1 Magnetic resonance

We now continue with demonstrating zero-field magnetometry by detecting the
change in linear polarization of the transmitted light. As mentioned, this is
done using polarimitry and a balanced photodetector. After having nulled the
total static magnetic field, we apply a linearly sweeping magnetic field in
the $\bm{e}_{x}$ direction $B_{x}^{\mathrm{sweep}}$ with strength ranging from
-150 nT to 150 nT and measure the signal for various applied
$B^{\mathrm{DC}}_{y}$ and $B^{\mathrm{DC}}_{z}$ static fields. These
measurements are presented in Figs. 5(a) and 5(b). As seen, the signal
exhibits a dispersive lineshape and is more or less unaffected by variations
in the static magnetic field components as predicted by our theoretical model,
see Eq. (16). I.e., in zero-field conditions, the magnetometer is only
sensitive to changes in the $\bm{e}_{x}$ component of the magnetic field.

####  IV.3.2 Bandwidth

We now stop sweeping any magnetic fields. The bandwidth of the zero-field
magnetometer is then assessed by sequentially applying 14 different sinusoidal
magnetic fields with constant amplitude 3.16 $\mathrm{nT_{rms}}$ and varying
frequency ranging from 8 Hz to 1818 Hz. Two additional measurements were
conducted in which the sinusoidal field was not applied (denoted “No RF”) and
where the balanced photodetector was disconnected from the data-acquisition
card (denoted “Electronic Noise Floor”). These two measurements were performed
to assess the noise level of the magnetometer and of the electronic noise of
the data acquisition card. For each of the $14+2=16$ measurements, the power
spectral density of the signal was calculated and plotted in Fig. 6. Each of
the 14 spectra have a sharp peak at the frequency corresponding to the
respective oscillation. Those peak heights were fitted as a function of
frequency to the response of a first order low-pass filter, i.e. the function
$K/\sqrt{1+\left(f/f_{c}\right)^{2}}$. The low-pass filter response function
has been derived for another type of zero-field OPM based on a spin-oriented
atomic ensemble [10] but is here used phenomenologically to describe the
response of our alignment based magnetometer. The fit allows us to extract the
cut-off frequency $f_{c}=\left(97\pm 3\right)$ Hz which is the magnetometer
bandwidth. We would expect this number to equal the half width at half maximum
of the magnetic resonances in Fig. 5 which is 11 nT corresponding to 39 Hz
(when multiplying with the gyromagnetic ratio), however, we do experimentally
find that the numbers differ by a factor of 2-3.

####  IV.3.3 Sensitivity

The noise floor of the magnetometer obtained with “No RF” depends on frequency
as seen from Fig. 6. It is $\approx 2$ $\mathrm{pT}/\sqrt{\mathrm{Hz}}$ at 1
Hz and then improving to $\approx 0.3$ $\mathrm{pT}/\sqrt{\mathrm{Hz}}$
between 10 Hz and 100 Hz. Our experimentally obtained sensitivity is mainly
limited by the electronic noise of our balanced photodetector (data not shown)
together with a contribution from the electronic noise of the data acquisition
card.

The sensitivity of an optically pumped magnetometer is fundamentally limited
by quantum fluctuations in the atomic spin and the light shot noise such that
the total quantum noise equals
$\delta\mathrm{B_{quantum}}=\sqrt{\delta\mathrm{B^{2}_{spin}}+\delta\mathrm{B^{2}_{shot}}}$.
The spin-projection noise can be written as [22]

$\delta\mathrm{B_{spin}}=\frac{2\hbar}{g_{F}\mu_{B}\sqrt{nVT_{2}}},$ (21)

where $\mu_{B}$ is the Bohr magneton, $g_{F}$ = 1/4 for the $F=4$ caesium
ground state, $n\approx 2.2\times 10^{16}\mathrm{m^{-3}}$ is the number
density of caesium atoms at room temperature (18.5∘C) [24] and
$V=(5\mathrm{mm})^{3}$ is the volume of our vapor cell. $T_{2}$ is the
transverse relaxation time and is calculated as $T_{2}=1/\Gamma\approx
1/\left[2\pi(11\mathrm{nT}\cdot 3.5\mathrm{Hz/nT})\right]\approx 4.1$ ms,
where we used the $\Gamma/\gamma\approx 11$ nT fitted in Fig. 5 and the
gyromagnetic ratio for caesium which can be expressed as $3.5$ Hz/nT. Ideally,
an optically pumped magnetometer has
$\delta\mathrm{B_{shot}}=\delta\mathrm{B_{spin}}$ such that the total quantum
noise $\delta\mathrm{B_{quantum}}=\sqrt{2}\delta\mathrm{B_{spin}}$. Hence the
quantum noise limit for our magnetometer is $\delta\mathrm{B_{quantum}}\approx
40~{}\mathrm{fT/\sqrt{Hz}}$.

Figure 6: Calculated power spectral density (averaged 10 times) using
$10\times 2$ s time traces of the magnetometer signal for 14 different
measurements where sinusoidal magnetic fields were applied with various
frequencies ranging from 8 Hz to 1818 Hz, and for two additional measurements
without a sine wave applied (“No RF”), and with the photodetector disconnected
from the data acquisition card (“Electronic Noise Floor”).

####  IV.3.4 Synthetic magnetocardiography signal

In order to determine if the performance of the magnetometer meets the
standards for biomagnetic measurements, particularly magnetocardiography
(MCG), a synthetic cardiac signal is applied in the $\bm{e}_{x}$ direction
with an amplitude of 100 pT peak-to-peak. This amplitude corresponds to the
magnetic field of an adult person’s heart [33, 34] The waveform was generated
using the stored electrocardiogram (ECG) waveform in a function generator
(RIGOL DG1032Z) connected to one of the coils. Figure 7a shows the raw
magnetometer signal, which is not averaged, and was recorded with a 10 kHz
sampling rate. The signal is relatively noisy, but simple signal post-
processing done by first a 49-51 Hz bandstop filter (to block 50 Hz noise),
see Fig. 7b, and then a 17 Hz low-pass filter (to block high-frequency noise),
see Fig. 7c, significantly improves the signal quality. We see that the
filtered magnetometer signal looks very similar to the ECG waveform from the
function generator. In particular, the real-time measurements captured by the
magnetometer reveal the presence of the QRS complex and P and T waves which
are features in the cardiogram which are important for medical doctors when
assessing the health of a heart.

Figure 7: Magnetometer time-trace with a synthetic cardiogram waveform
applied. (a) Raw data, (b) data processed with a 49-51 Hz bandpass filter, and
(c) with additionally a 17 Hz low-pass filter.

##  V CONCLUSIONS

In this study, we investigate the use of spin-aligned atomic ensembles for
zero-field magnetometry. A novel method is introduced, which involves
assessing the polarization rotation in transmitted light through an atomic
vapor. We derive analytical expressions for the generated spin alignment and
photodetection signals. By analyzing the obtained polarization rotation
signal, we null the magnetic field components through an iterative procedure
that adjusts the applied magnetic fields until a zero signal is achieved,
indicating the absence of any static magnetic field. Achieving a zero-field
environment is important for several applications, including nuclear magnetic
resonance (NMR) gyroscopes [35, 36], and for biomedical measurements, such as
magnetoencephalography (MEG) and magnetocardiography (MCG) where non-zero
field components leads to artifacts or systematic errors [37]. Our method for
nulling the field could potentially be employed in such applications, as
demonstrated here by successfully detecting a synthetic cardiac signal. As
opposed to other methods, e.g. using circular polarized light with several
modulation fields, our method is modulation-free. After nulling the field, our
magnetometer measures one component of the field using a single laser beam. If
instead of using a paraffin-coated vapour cell (where all atoms inside the
cell are probed), one uses a buffer-gas cell [24] (where only atoms inside the
beam volume are probed), a vector magnetometer could in principle be
implemented using three non-overlapping laser beams. Future work could
investigate how to use only one or two laser linearly polarized beams for
vector measurements.

We demonstrated a magnetic field sensitivity of $\approx 2$
$\mathrm{pT}/\sqrt{\mathrm{Hz}}$ at 1 Hz, which further improved to about 0.3
$\mathrm{pT}/\sqrt{\mathrm{Hz}}$ between 10 Hz and 100 Hz. The quantum limited
sensitivity was calculated to $\delta\mathrm{B_{quantum}}\approx
40~{}\mathrm{fT/\sqrt{Hz}}$ meaning that further technical improvements are
possible. Notably, some commercial optically pumped magnetometers have reached
a sensitivity level of around $\sim 20$ fT/$\sqrt{\mathrm{Hz}}$ [5], however
they are operated at an elevated temperature $\sim 150^{\circ}$C, as opposed
to our magnetometer which was operated at room temperature. To enhance our
sensitivity, various strategies can be explored, including the utilization of
photodetectors with reduced electronic noise, multiple passes of the laser
beam through the vapor cell, employing a larger vapor cell, increasing atom
density through vapor cell heating, and using a repump laser beam to transfer
atoms from the “dark” $F=3$ ground state back to the $F=4$ ground state, which
is not probed by the laser light.

## ACKNOWLEDGEMENTS

This work was supported by the QuantERA grant C’MON-QSENS! funded by the
Engineering and Physical Sciences Research Council (EPSRC) (Grant No.
EP/T027126/1), the Novo Nordisk Foundation (Grant No. NNF20OC0064182), and the
University of Nottingham Hermes Award. Project C’MON-QSENS! is supported by
the National Science Centre (2019/32/Z/ST2/00026), Poland under QuantERA,
which has received funding from the European Union’s Horizon 2020 research and
innovation programme under grant agreement no. 731473.

## References

  * Budker and Romalis [2007] D. Budker and M. Romalis, Nat. Phys. 3, 227 (2007).
  * Wilson  _et al._ [2019] N. Wilson, P. Light, A. Luiten, and C. Perrella, Phys. Rev. Applied 11, 044034 (2019).
  * Kominis  _et al._ [2003] I. Kominis, T. Kornack, J. Allred, and M. V. Romalis, Nature 422, 596 (2003).
  * Fagaly [2006] R. L. Fagaly, Rev. Sci. Instrum. 77, 101101 (2006).
  * Qus [2023] “Quspin,” https://quspin.com (2023).
  * fie [2023] “Fieldline inc.” https://fieldlineinc.com/ (2023).
  * twi [2023] “Twinleaf,” https://twinleaf.com/ (2023).
  * Li  _et al._ [2023] R. Li, Y. Liu, L. Cao, S. Li, J. Li, and Y. Zhai, Measurement 214, 112733 (2023).
  * Le Gal and Palacios-Laloy [2022] G. Le Gal and A. Palacios-Laloy, Phys. Rev. A 105, 043114 (2022).
  * Jensen  _et al._ [2018] K. Jensen, M. A. Skarsfeldt, H. Stærkind, J. Arnbak, M. V. Balabas, S.-P. Olesen, B. H. Bentzen, and E. S. Polzik, Sci. Rep. 8, 16218 (2018).
  * Marmugi and Renzoni [2016] L. Marmugi and F. Renzoni, Sci. Rep. 6, 23962 (2016).
  * Brookes  _et al._ [2022] M. J. Brookes, J. Leggett, M. Rea, R. M. Hill, N. Holmes, E. Boto, and R. Bowtell, Trends in Neurosciences (2022).
  * Xia  _et al._ [2006] H. Xia, A. Ben-Amar Baranga, D. Hoffman, and M. Romalis, Appl. Phys. Lett. 89, 211104 (2006).
  * Sander  _et al._ [2012] T. Sander, J. Preusser, R. Mhaskar, J. Kitching, L. Trahms, and S. Knappe, Biomed. Opt. Express 3, 981 (2012).
  * Wyllie  _et al._ [2012] R. Wyllie, M. Kauer, R. T. Wakai, and T. G. Walker, Opt. Lett. 37, 2247 (2012).
  * Rushton  _et al._ [2022] L. M. Rushton, T. Pyragius, A. Meraki, L. Elson, and K. Jensen, Rev. Sci. Instrum. 93, 125103 (2022).
  * Bevington  _et al._ [2019] P. Bevington, R. Gartman, and W. Chalupczak, J. Appl. Phys. 125 (2019).
  * Romalis and Dang [2011] M. V. Romalis and H. B. Dang, Materials today 14, 258 (2011).
  * Boto  _et al._ [2018] E. Boto, N. Holmes, J. Leggett, G. Roberts, V. Shah, S. S. Meyer, L. D. Muñoz, K. J. Mullinger, T. M. Tierney, S. Bestmann, _et al._ , Nature 555, 657 (2018).
  * Borna  _et al._ [2020] A. Borna, T. R. Carter, A. P. Colombo, Y.-Y. Jau, J. McKay, M. Weisend, S. Taulu, J. M. Stephen, and P. D. Schwindt, Plos one 15, e0227684 (2020).
  * Weis  _et al._ [2006] A. Weis, G. Bison, and A. S. Pazgalev, Phys. Rev. A 74, 033401 (2006).
  * Ledbetter  _et al._ [2007] M. P. Ledbetter, I. M. Savukov, V. M. Acosta, D. Budker, and M. Romalis, Phys. Rev. A 75, 032107 (2007).
  * Di Domenico  _et al._ [2007] G. Di Domenico, H. Saudan, G. Bison, P. Knowles, and A. Weis, Phys. Rev. A 76, 023407 (2007).
  * Rushton  _et al._ [2023] L. M. Rushton, L. Elson, A. Meraki, and K. Jensen, Phys. Rev. Applied 19, 064047 (2023).
  * Savukov  _et al._ [2005] I. M. Savukov, S. Seltzer, M. Romalis, and K. Sauer, Phys. Rev. Lett. 95, 063004 (2005).
  * Zigdon  _et al._ [2010] T. Zigdon, A. Wilson-Gordon, S. Guttikonda, E. Bahr, O. Neitzke, S. Rochester, and D. Budker, Opt. Express 18, 25494 (2010).
  * Ingleby  _et al._ [2018] S. J. Ingleby, C. O’Dwyer, P. F. Griffin, A. S. Arnold, and E. Riis, Phys. Rev. Applied 10, 034035 (2018).
  * Beato  _et al._ [2018] F. Beato, E. Belorizky, E. Labyt, M. Le Prado, and A. Palacios-Laloy, Phys. Rev. A 98, 053431 (2018).
  * Bertrand  _et al._ [2021] F. Bertrand, T. Jager, A. Boness, W. Fourcault, G. Le Gal, A. Palacios-Laloy, J. Paulet, and J. M. Léger, Rev. Sci. Instrum. 92, 105005 (2021).
  * Wang  _et al._ [2022] B. Wang, X. Peng, H. Wang, W. Xiao, and H. Guo, Sensors 22, 4184 (2022).
  * Steck [2023] D. A. Steck, “Cs D Line Data,” https://steck.us/alkalidata/cesiumnumbers.1.6.pdf (2023).
  * Auzinsh  _et al._ [2014] M. Auzinsh, D. Budker, and S. M. Rochester,  _Optically polarized atoms: understanding light-atom interactions_ (Oxford University Press, 2014).
  * Bison  _et al._ [2009] G. Bison, N. Castagna, A. Hofer, P. Knowles, J.-L. Schenker, M. Kasprzak, H. Saudan, and A. Weis, Appl. Phys. Lett. 95, 173701 (2009).
  * Sander  _et al._ [2020] T. Sander, A. Jodko-Władzińska, S. Hartwig, R. Brühl, and T. Middelmann, Advanced Optical Technologies 9, 247 (2020).
  * Walker and Larsen [2016] T. Walker and M. Larsen (Academic Press, 2016) pp. 373–401.
  * Eklund [2008] E. J. Eklund,  _Microgyroscope based on spin-polarized nuclei_ , Ph.D. thesis, University of California, Irvine (2008).
  * Borna  _et al._ [2022] A. Borna, J. Iivanainen, T. R. Carter, J. McKay, S. Taulu, J. Stephen, and P. D. Schwindt, NeuroImage 247, 118818 (2022).


# Long-term stability and generalization of observationally-constrained
stochastic data-driven models for geophysical turbulence

Ashesh Chattopadhyay
Rice University
Houston, TX 77005
&Jaideep Pathak
NVIDIA Corporation
Santa Clara, CA 95051
&Ebrahm Nabizadeh
Rice University
Houston, TX 77005
&Wahid Bhimji
Lawrence Berkeley
National Laboratory
Berkeley, CA 94720
Pedram Hassanzadeh
Rice University
Houston, TX 77005

###### Abstract

Recent years have seen a surge in interest in building deep learning-based
fully data-driven models for weather prediction. Such deep learning models if
trained on observations can mitigate certain biases in current state-of-the-
art weather models, some of which stem from inaccurate representation of
subgrid-scale processes. However, these data-driven models, being over-
parameterized, require a lot of training data which may not be available from
reanalysis (observational data) products. Moreover, an accurate, noise-free,
initial condition to start forecasting with a data-driven weather model is not
available in realistic scenarios. Finally, deterministic data-driven
forecasting models suffer from issues with long-term stability and unphysical
climate drift, which makes these data-driven models unsuitable for computing
climate statistics. Given these challenges, previous studies have tried to
pre-train deep learning-based weather forecasting models on a large amount of
imperfect long-term climate model simulations and then re-train them on
available observational data. In this paper, we propose a convolutional
variational autoencoder-based stochastic data-driven model that is pre-trained
on an imperfect climate model simulation from a 2-layer quasi-geostrophic flow
and re-trained, using transfer learning, on a small number of noisy
observations from a perfect simulation. This re-trained model then performs
stochastic forecasting with a noisy initial condition sampled from the perfect
simulation. We show that our ensemble-based stochastic data-driven model
outperforms a baseline deterministic encoder-decoder-based convolutional model
in terms of short-term skills while remaining stable for long-term climate
simulations yielding accurate climatology.

##  1 Introduction

A surge of interest in building deep learning-based data-driven models for
complex systems such as chaotic dynamical systems [1, 2], fully turbulent flow
[3], and weather and climate models [4, 5, 6, 7, 8] has been seen in the
recent past. This interest in data-driven modeling stems from the hope that if
these data-driven models are trained on observational data, (i) they would not
suffer from some of the biases of physics-based numerical climate and weather
models [9, 10], e.g., due to ad-hoc parameterizations, (ii) they can be used
to generate large ensemble forecasts for data assimilation [8], and (iii) they
might seamlessly be integrated to perform climate simulations [11] that would
allow for generating large synthetic datasets to study extreme events [12].

Despite the promise of data-driven weather prediction (DDWP) models, most DDWP
models cannot compete with state-of-the-art numerical weather prediction (NWP)
models [13, 7, 4, 5, 8] although more recently FourCastNet [14] has shown
promise of being quite competitive with the best NWP models. There have also
been some previous work, where pre-training models on climate simulations have
also resulted in improved short-term forecasts on reanalysis datasets [7].
Moreover, these DDWP models, while comparable with NWP in terms of short-term
forecasts, are unstable when integrated for longer time scales. These
instabilities either appear in the form of blow-ups or unphysical climate
drifts [13, 3].

In this paper, we consider the following realistic problem set-up, wherein, we
have long-term climate simulations from an imperfect 2-layer quasi-geostrophic
(QG) flow (we will call it "imperfect system" from here on) at our disposal to
pre-train a DDWP. Following that, we have a few noisy observations from a
perfect 2-layer QG flow ("perfect system" from here on) on which we are able
to fine-tune our DDWP. This observationally-constrained DDWP model is then
expected to generalize to the perfect system with a noisy initial condition
(obtained from the available observations) and perform stochastic short-term
forecasting. It is also expected to be seamlessly time-integrated to obtain
long-term climatology. We show that in this regard, a stochastic forecasting
approach with a generative model such as a convolutional variational
autoencoder (VAE) outperforms a deterministic encoder-decoder-based
convolutional architecture in terms of short-term forecasts and remains stable
without unphysical climate drift for long-term integration. A deterministic
encoder-decoder would fail to do so and shows unphysical climate drift
followed by numerical blow-up.

##  2 Methods

###  2.1 Imperfect and Perfect Systems

In this paper, we consider 2-layer QG flow as the geophysical system on which
we intend to perform both short-term forecasting and compute long-term
statistics.

The dimensionless dynamical equations of the 2-layer QG flow are derived
following Lutsko et al. [15], and Nabizadeh et al. [16]. The system consists
of two constant density layers with a $\beta$-plane approximation in which the
meridional temperature gradient is relaxed towards an equilibrium profile. The
equation of the system is described as:

$\displaystyle\begin{split}\frac{\partial q_{j}}{\partial
t}+J(\psi_{j},q_{j})=-\frac{1}{\tau_{d}}(-1)^{j}\left(\psi_{1}-\psi_{2}-\psi_{R}\right)\\\
-\frac{1}{\tau_{f}}\delta_{k2}\nabla^{2}\psi_{j}-\nu\nabla^{8}q_{j}.\end{split}$
(1)

Here, $q$ is potential vorticity and is expressed as

$\displaystyle
q_{j}=\nabla^{2}\psi_{j}+(-1)^{j}\left(\psi_{1}-\psi_{2}\right)+\beta y,$ (2)

where $\psi_{j}$ is the stream function of the system. In Eqs. (1) and (2),
$j$ denotes the upper ($j=1$) and lower ($j=2$) layers. $\tau_{d}$ is the
Newtonian relaxation time scale while $\tau_{f}$ is the Rayleigh friction time
scale, which only acts on the lower layers represented by the Kronecker
$\delta$ function ($\delta_{k2}$). $J$ denotes the jacobian, $\beta$ is the
$y-$gradient of the coriolis parameter, and $\nu$ denotes the hyperdiffusion
coefficient. We have induced a baroclinically unstable jet at the center of a
horizontally (zonally) periodic channel by setting $\psi_{1}-\psi_{2}$ to be
equal to a hyperbolic secant centered at $y=0$. When eddy fluxes are absent,
$\psi_{2}$ is identically zero, making zonal velocity in the upper layer
$u(y)=-\frac{\partial\psi_{1}}{\partial y}=-\frac{\partial\psi_{R}}{\partial
y}$ where we set

$\displaystyle-\frac{\partial\psi_{R}}{\partial
y}=\operatorname{sech}^{2}\left(\frac{y}{\sigma}\right),$ (3)

$\sigma$ being the width of the jet. Parameters of the model are set following
the previous studies [15, 16], in which $\beta=0.19$, $\sigma=3.5$,
$\tau_{f}=15$, and $\tau_{d}=100$.

To non-dimensionalize the equations, we have used the maximum strength of the
equilibrium velocity profile as the velocity scale ($U$) and the deformation
radius ($L$) for the length scale. The system’s time scale ($L/U$) is referred
to as the “advective time scale” ($\tau_{adv}$).

The spatial discretization is spectral in both $x$ and $y$ where we have
retained $96$ and $192$ Fourier modes, respectively. The length and width of
the domain is equal to $46$ and $68$ respectively after non-dimensionalizing
the numbers. The spurious waves on the northern and southern boundaries are
damped by applying sponge layers. Note that the domain is wide enough that the
sponges do not affect the dynamics. Here $5\tau_{adv}\approx$ $1$ Earth day
$\approx 200\Delta t_{n}$, where $\Delta t_{n}=0.025$ is the time step of the
leapfrog time integrator used in the numerical scheme.

The imperfect system in this paper has an increased value of $\beta$, given by
$\beta^{*}=3\beta$ and a decreased size of the jet, given by
$\sigma*=\frac{4}{5}\sigma$. We assume that we have long-term climate
simulations of the imperfect system while we have only a few noisy
observations of the perfect system. The difference in the long-term averaged
zonal-mean velocity in the upper layer of the perfect and imperfect systems is
shown in Fig. 1. It is to be noted that the motivation behind having an
imperfect and perfect system is the fact that we assume that the imperfect
system represents a climate model with biases from which we have long
simulations at our disposal, while the perfect system is analogous to actual
observations from our atmosphere, where the sample size is small and
contaminated with measurement noise.

Figure 1: Time-averaged zonal-mean velocity (over $20000$ days),
$\left<\bar{u}\right>$, of the imperfect and perfect system. The difference
between $\left<\bar{u}\right>$ of the perfect and imperfect systems indicates
a challenge for DDWP models to seamlessly generalize from one system to the
other

###  2.2 VAE and Transfer Learning Frameworks

Here, we propose a generative modeling approach to perform both short- and
long-term forecasting for the perfect QG system after being trained on the
imperfect QG system. We assume that we have a long-term climate simulation
from the imperfect system, with states $x^{m}(k\Delta t)$, where
$k\in{1,2,\cdots K}$ and $x^{m}\in\mathcal{R}^{2\times 96\times 192}$. A
convolutional VAE is trained on the states of this imperfect system to predict
an ensemble of states with a probability density function,
$p\left(x^{m}\left(t+\Delta t\right)\right)$, from $x^{m}(t)$. Here, $\Delta
t=40\Delta t_{n}$, which is also the sampling interval for training the model.
Extensive trial and error-based search has been performed for hyperparameter
optimization for the VAE. The details on the VAE architecture is given in
Table 1.

This pre-trained VAE would then undergo transfer learning on a small number of
noisy observations of the states of the perfect system $x^{o}(n\Delta t)$
where $n\in{1,2,\cdots N}$ and $N\ll K$. The VAE would then stochastically
forecast an ensemble of states of the perfect system with a noisy initial
condition from the perfect system. The noise in the initial condition is
sampled from a Gaussian normal distribution with $0$ mean and standard
deviation being $\eta\sigma_{Z}$, where $\sigma_{Z}$ is the standard deviation
of $\psi_{k}$ and $\eta$ is a fraction that determines the amplitude of the
noise vector. A schematic of this framework is shown in Fig. 2. In this work,
we consider the stream function $\psi_{k}$ ($k=1,2$) to be the states of the
system on which we would train our VAE.

The VAE is trained on $9$ independent ensembles of the imperfect system with
$1400$ consecutive days of climate simulation each with a sampling interval of
$\Delta t$. Hence, the training size for the VAE is about $12600$ days of
data. For transfer learning, we assume that only $10\%$ of the number of
training samples is available as noisy observations from the perfect system.

Number | Layer |  Number of Filters
---|---|---
1 |  $5\times 5$ 2D Convolution |  32
2 |  $2\times 2$ Max Pooling |  –
3 |  $5\times 5$ 2D Convolution |  32
4 |  $2\times 2$ Max Pooling |  –
5 |  $5\times 5$ 2D Convolution |  32
6 |  $2\times 2$ Max Pooling |  –
7 | flatten |  –
8 | 128 neurons for mean and standard deviation |  –
9 | Dense layer |  –
10 |  $5\times 5$ 2D Convolution |  32
11 | Up-sampling |  –
12 |  $5\times 5$ 2D Convolution |  32
13 | Up-sampling |  –
14 |  $5\times 5$ 2D Convolution |  32
15 | Up-sampling |  –
15 |  $5\times 5$ 2D Convolution |  32
16 | Up-sampling |  –
17 |  $5\times 5$ 2D Convolution |  2
Table 1: Number of layers and filters in the convolutional VAE architecture
used in this paper. The baseline convolutional encoder-decoder model would
have the exact same architecture and the same latent space size as that of the
VAE.  Figure 2: A schematic of the transfer learning framework with a VAE pre-
trained on imperfect simulations and transfer learned on noisy observations
from the perfect system. Here, $x^{m}\left(k\Delta t\right)$ are states
obtained from the imperfect system and $x^{o}\left(n\Delta t\right)$ are noisy
observations from the perfect system. Here, $\Delta t=40\Delta t_{n}$. Note
that, our proposed VAE is convolutional and the schematic is just
representative. Details about the architecture are given in Table 1

##  3 Results

###  3.1 Short-term Stochastic Forecasting on the Imperfect System

Figure 3: Performance of the stochastic convolutional VAE as compared to
baseline convolutional encoder-decoder-based model when trained on the
imperfect system and predicts from a noisy initial condition sampled from the
imperfect system. (a) Anomaly correlation coefficient (ACC) between predicted
$\psi_{1}$ and true $\psi_{1}$ [17] from the imperfect system with $\eta=5\%$
(noise level) for the initial condition. (b) Same as (a) but for RMSE. (c)
Prediction horizon (number of $\Delta t_{n}$ until ACC $\leq 0.60$) for
different noise levels added to the initial condition. (d) Averaged error over
$2$ days for different noise levels of initial condition. The shading shows
the standard deviation across $100$ ensembles generated by the VAE model
during inference

In this section, we show how well the stochastic VAE performs in terms of
short-term forecasting on the imperfect system when the initial condition is
sampled from the imperfect system. Our baseline model is a convolutional
encoder-decoder-based deterministic model that has the exact same architecture
as that of the convolutional VAE with the same size of the latent space.
During inference, the stochastic VAE generates $100$ ensembles at every
$\Delta t$ and the mean of these ensembles is fed back into the VAE
autoregressively to predict future time steps. We see from Fig. 3 that the
stochastic VAE model outperforms the deterministic model. Our experiments
suggest that increasing the noise level of the initial condition does not
affect the prediction horizon of VAE if we consider ACC $=0.60$ as the limit
of prediction (although it is an ad-hoc choice) as shown in Fig. 3(c). This
maybe attributed to the reduction of initial condition error (due to noise)
with ensembling. The averaged RMSE error over $2$ days grows faster with an
increase in the noise level in the deterministic model as compared to the
stochastic VAE model as shown in Fig. 3(d).

###  3.2 Short-term Stochastic Forecasting on the Perfect System

In this section, we show how well the pre-trained VAE (as compared to the
baseline encoder-decoder and imperfect numerical model) on the imperfect
climate simulations performs on the perfect system when constrained by noisy
observations from the perfect system on which it is re-trained. Only two
layers before and after the bottleneck layer are re-trained for both the VAE
and the baseline model. Figure 4 (a) and (b) show that for short-term
forecasting, VAE outperforms both the baseline deterministic encoder-decoder
model as well as the imperfect numerical model. Figure 4(c) shows a similar
behaviour for VAE with increase in noise level in the initial condition as
Fig. 3(c). Figure 4(d) shows the effect of the number of re-training samples
used as noisy observations for transfer learning on the prediction horizon.
While the numerical model is unaffected by the impact of noisy observations
(since it does not undergo any learning), an increase in noise level to the
initial condition makes it more susceptible to error as compared to the data-
driven models which are inherently more robust to noise (due to re-training on
noisy observations). In all cases, the VAE outperforms both the deterministic
encoder-decoder model and the imperfect numerical model for short-term
forecasts.

Figure 4: Performance of the stochastic convolutional VAE compared to that of
the baseline convolutional encoder-decoder-based model and the imperfect
numerical model when trained on the imperfect system, re-trained on noisy
observations from the perfect system, and initialized with a noisy initial
condition sampled from the perfect system. (a) Anomaly correlation coefficient
(ACC) between predicted $\psi_{1}$ and true $\psi_{1}$ from the perfect system
with $\eta=5\%$ (noise level) for the initial condition. (b) Same as (a) but
for RMSE. (c) Prediction horizon (number of $\Delta t_{n}$ until ACC $\leq
0.60$) for different noise levels added to initial condition. (d) Prediction
horizon of the models with different sample sizes of noisy observation data
(as percentages of original training sample size) from perfect system. Note
that the numerical model is not trained on data and hence would show no
effect.

###  3.3 Long-term Climate Statistics

In this section, we show the comparison between long-term climatology obtained
from seamlessly integrating the VAE for 20000 days and the true climatology of
the imperfect system. It must be noted that the deterministic encoder-decoder
model is not stable and would show unphysical drifts in the climate similar to
other deterministic data-driven models within a few days of seamless
integration [3]. Figure. 5 shows stable and non-drifting physical climate
obtained from the VAE whose mean of $\psi_{j}$, $u_{j}$, and EOF1 of $u_{1}$
closely match those of the true system.

Figure 5: Long-term climatology obtained from integrating the models for
$20000$ days. $\left<\bar{\psi}_{1}\right>$ is the time mean zonally-averaged
$\psi_{1}$, $\left<\bar{u}\right>$ is the time mean zonally averaged upper-
layer velocity, and EOF1, is the first EOF (empiricial orthogonal function) of
$\left<\bar{u}\right>$. The VAE shows non-drifting physical climate. Baseline
encoder-decoder model is not shown since it becomes unstable within $100$ days
of seamless integration.

##  4 Conclusion

In this paper, we have developed a convolutional VAE for stochastic short-term
forecasting and long-term climatology of the stream function of a 2-layer QG
flow. The VAE is trained on climate simulations from an imperfect QG system
before being fine-tuned on noisy observations from a perfect QG system on
which it performs both short- and long-term forecasting. The VAE outperforms
both the baseline convolutional encoder-decoder model as well as the imperfect
numerical model.

One of the main advantages of using the stochastic VAE with multiple ensemble
members is the reduction of the effect of noise in the initial condition for
short-term forecasting. Moreover, the VAE remains stable through a seamless
$20000$ days integration yielding a physical climatology that matches the
numerical model. Deterministic data-driven models do not remain stable for
long-term integration and hence a stochastic generative modeling approach may
be a potential candidate when developing data-driven weather forecasting
models that can seamlessly be integrated to yield long-term climate
simulations.

Despite the improvement in performance in short-term skills as well as long-
term stability, one of the caveats of the VAE is that it is less interpretable
as compared to a deterministic model. Moreover, the VAE is more difficult to
scale as compared to a deterministic model, especially for high-dimensional
systems. This is because one would need to evolve a large number of ensembles
for high-dimensional systems, and the exact relationship between short-term
skills and long-term stability with the number of ensembles is not well known.
Finally, while the application of VAE has resulted in long-term stability, the
underlying causal mechanism by which instability is exhibited in a
determinisitc data-driven model is still largely unknown. Future work needs to
be undertaken to understand the cause of this instability such that
determinisitc models that may be more scalable than VAE can also be used for
data-driven weather and climate prediction.

#### Acknowledgments

This work was started as an internship project by A.C. in Lawrence Berkeley
National Laboratory in the Summer of 2021 under J.P. and W.B. It was then
continued under P.H. as a part of A.C.’s PhD at Rice University. This research
used resources of NERSC, a U.S. Department of Energy Office of Science User
Facility operated under Contract No. DE-AC02-05CH11231. A.C., E.N., and P.H.
were supported by ONR grant N00014-20-1-2722, NASA grant 80NSSC17K0266, and
NSF CSSI grant OAC-2005123. A.C. also thanks the Rice University Ken Kennedy
Institute for a BP HPC Graduate Fellowship. We are grateful to Kamyar
Azizzadenesheli, Mustafa Mustafa, and Sanjeev Raja for providing invaluable
advice for this work.

#### Competing Interests

None

#### Data Availability Statement

The codes for this work would be available to the public upon publication in
this Github link: https://github.com/ashesh6810/Stochastic-VAE-for-Digital-
Twins.git.

#### Ethical Standards

The research meets all ethical guidelines, including adherence to the legal
requirements of the study country.

## References

  * Pathak et al. [2018] Jaideep Pathak, Brian Hunt, Michelle Girvan, Zhixin Lu, and Edward Ott.  Model-free prediction of large spatiotemporally chaotic systems from data: A reservoir computing approach.  _Physical Review Letters_ , 120(2):024102, 2018\.
  * Chattopadhyay et al. [2020a] Ashesh Chattopadhyay, Pedram Hassanzadeh, and Devika Subramanian.  Data-driven predictions of a multiscale Lorenz 96 chaotic system using machine-learning methods: reservoir computing, artificial neural network, and long short-term memory network.  _Nonlinear Processes in Geophysics_ , 27(3):373–389, 2020a.
  * Chattopadhyay et al. [2020b] Ashesh Chattopadhyay, Mustafa Mustafa, Pedram Hassanzadeh, and Karthik Kashinath.  Deep spatial transformers for autoregressive data-driven forecasting of geophysical turbulence.  In _Proceedings of the 10th International Conference on Climate Informatics_ , pages 106–112, 2020b.
  * Weyn et al. [2019] Jonathan A Weyn, Dale R Durran, and Rich Caruana.  Can machines learn to predict weather? using deep learning to predict gridded 500-hPa geopotential height from historical weather data.  _Journal of Advances in Modeling Earth Systems_ , 11(8):2680–2693, 2019.
  * Weyn et al. [2020] Jonathan A Weyn, Dale R Durran, and Rich Caruana.  Improving data-driven global weather prediction using deep convolutional neural networks on a cubed sphere.  _Journal of Advances in Modeling Earth Systems_ , 12(9):e2020MS002109, 2020.
  * Weyn et al. [2021] Jonathan A Weyn, Dale R Durran, Rich Caruana, and Nathaniel Cresswell-Clay.  Sub-seasonal forecasting with a large ensemble of deep-learning weather prediction models.  _arXiv preprint arXiv:2102.05107_ , 2021.
  * Rasp and Thuerey [2021] Stephan Rasp and Nils Thuerey.  Data-driven medium-range weather prediction with a resnet pretrained on climate simulations: A new model for weatherbench.  _Journal of Advances in Modeling Earth Systems_ , page e2020MS002405, 2021.
  * Chattopadhyay et al. [2022] Ashesh Chattopadhyay, Mustafa Mustafa, Pedram Hassanzadeh, Eviatar Bach, and Karthik Kashinath.  Towards physics-inspired data-driven weather forecasting: integrating data assimilation with a deep spatial-transformer-based U-NET in a case study with ERA5.  _Geoscientific Model Development_ , 15(5):2221–2237, 2022.
  * Balaji [2021] V Balaji.  Climbing down charney’s ladder: machine learning and the post-dennard era of computational climate science.  _Philosophical Transactions of the Royal Society A_ , 379(2194):20200085, 2021.
  * Schultz et al. [2021] MG Schultz, C Betancourt, B Gong, F Kleinert, M Langguth, LH Leufen, Amirpasha Mozaffari, and S Stadtler.  Can deep learning beat numerical weather prediction?  _Philosophical Transactions of the Royal Society A_ , 379(2194):20200097, 2021.
  * Watson-Parris [2021] Duncan Watson-Parris.  Machine learning for weather and climate are worlds apart.  _Philosophical Transactions of the Royal Society A_ , 379(2194):20200098, 2021.
  * Chattopadhyay et al. [2020c] Ashesh Chattopadhyay, Ebrahim Nabizadeh, and Pedram Hassanzadeh.  Analog forecasting of extreme-causing weather patterns using deep learning.  _Journal of Advances in Modeling Earth Systems_ , 12(2):e2019MS001958, 2020c.
  * Scher and Messori [2019] Sebastian Scher and Gabriele Messori.  Weather and climate forecasting with neural networks: using general circulation models (GCMs) with different complexity as a study ground.  _Geoscientific Model Development_ , 12(7):2797–2809, 2019.
  * Pathak et al. [2022] Jaideep Pathak, Shashank Subramanian, Peter Harrington, Sanjeev Raja, Ashesh Chattopadhyay, Morteza Mardani, Thorsten Kurth, David Hall, Zongyi Li, Kamyar Azizzadenesheli, et al.  Fourcastnet: A global data-driven high-resolution weather model using adaptive Fourier neural operators.  _arXiv preprint arXiv:2202.11214_ , 2022.
  * Lutsko et al. [2015] Nicholas J Lutsko, Isaac M Held, and Pablo Zurita-Gotor.  Applying the fluctuation–dissipation theorem to a two-layer model of quasigeostrophic turbulence.  _Journal of the Atmospheric Sciences_ , 72(8):3161–3177, 2015.
  * Nabizadeh et al. [2019] Ebrahim Nabizadeh, Pedram Hassanzadeh, Da Yang, and Elizabeth A Barnes.  Size of the atmospheric blocking events: Scaling law and response to climate change.  _Geophysical Research Letters_ , 46(22):13488–13499, 2019.
  * Murphy and Epstein [1989] Allan H Murphy and Edward S Epstein.  Skill scores and correlation coefficients in model verification.  _Monthly weather review_ , 117(3):572–582, 1989\.


# Fermion Tunneling Beyond Semiclassical Approximation

Bibhas Ranjan Majhi
S. N. Bose National Centre for Basic Sciences,
JD Block, Sector III, Salt Lake, Kolkata-700098, India
E-mail<EMAIL_ADDRESS>
Abstract:
Applying the Hamilton-Jacobi method beyond the semiclassical approximation
prescribed in [12] for the scalar particle, Hawking radiation as tunneling of
Dirac particle through an event horizon is analysed. We show that, as before,
all quantum corrections in the single particle action are proportional to the
usual semiclassical contribution. We also compute the modifications to the
Hawking temperature and Bekenstein-Hawking entropy for the Schwarzschild black
hole. Finally, the coefficient of the logarithmic correction to entropy is
shown to be related with the trace anomaly.

Semiclassical methods of modeling Hawking radiation as a tunneling effect were
developed over the past decade and have generated a lot of interest [1, 2, 3,
4, 5, 6, 7, 8, 9, 10]. From this approach an alternative (intuitive) way of
understanding black hole radiation emerged. However, most of the calculations
in the literature [1, 2, 4, 5, 6, 7, 8, 9, 10] have been performed just for
scalar particles. Since a black hole can radiate all types of particles like a
black body, the emission spectrum should contain particles of all spins.
Therefore a detailed study of spin one-half particle emission is necessary.
Although there exist some computations [11] in this context, these are
confined to the semiclassical approximation and do not consider quantum
corrections.

In our previous work [12], we formulated the Hamilton-Jacobi method of
tunneling beyond semiclassical approximation by considering all the terms in
the expansion of the one particle action for a scalar particle. We showed that
the higher order terms are proportional to the semiclassical contribution.
This result, together with properties of conformal transformations, eventually
led to corrected expressions for thermodynamic variables of a black hole. It
is not obvious whether a similar analysis is valid for the case of spin-half
fermion tunneling. This issue is addressed here.

In this paper we will discuss the Dirac particle tunneling beyond
semiclassical approximation employing the Hamilton-Jacobi method suggested in
[12]. We will explicitly show that the higher order terms in the single
particle action are again proportional to the semiclassical contribution. By
dimensional argument the form of these proportionality constants, upto some
dimensionless parameters, are determined. In particular for Scwarzschild
spacetime, these are given by the inverse powers of the square of the mass of
the black hole, because in this case, the only macroscopic parameter is mass.
Using the principle of “detailed balance” [2, 9] the modified Hawking
temperature is identified. Then the corrections to the Bekenstein-Hawking area
law are derived by using the Gibbs form of first law of thermodynamics.
Interestingly, the leading order correction to the entropy is the logarithmic
of the semiclassical entropy which was found earlier in [13, 14, 15, 16, 17,
18, 19, 20, 21, 22]. Finally, using a constant scale transformation to the
metric, we show that the coefficient of the logarithmic correction is related
to trace anomaly.

Our method involves calculating the imaginary part of the action for the
(classically forbidden) process of s-wave emission across the horizon which in
turn is related to the Boltzmann factor for emission at the Hawking
temperature. We consider a massless Dirac particle in a general class of
static, spherically symmetric spacetime of the form

$\displaystyle ds^{2}=-f(r)dt^{2}+\frac{dr^{2}}{g(r)}+r^{2}d\Omega^{2}$ (1)

where the horizon $r=r_{H}$ is given by $f(r_{H})=g(r_{H})=0$. The massless
Dirac equation is given by

$\displaystyle i\gamma^{\mu}\nabla_{\mu}\psi=0$ (2)

where for this case the $\gamma$ matrices are defined as,

$\displaystyle\gamma^{t}$ $\displaystyle=$
$\displaystyle\frac{1}{\sqrt{f(r)}}\left(\begin{array}[]{c c}i&0\\\
0&-i\end{array}\right);\,\,\ \gamma^{r}=\sqrt{g(r)}\left(\begin{array}[]{c
c}0&\sigma^{3}\\\ \sigma^{3}&0\end{array}\right)$ (7)
$\displaystyle\gamma^{\theta}$ $\displaystyle=$
$\displaystyle\frac{1}{r}\left(\begin{array}[]{c c}0&\sigma^{1}\\\
\sigma^{1}&0\end{array}\right);\,\,\,\
\gamma^{\phi}=\frac{1}{r\textrm{sin}\theta}\left(\begin{array}[]{c
c}0&\sigma^{2}\\\ \sigma^{2}&0\end{array}\right).$ (12)

The covariant derivative is given by,

$\displaystyle\nabla_{\mu}=\partial_{\mu}+\frac{i}{2}\Gamma{{}^{\alpha}}{{}_{\mu}}{{}^{\beta}}\Sigma_{\alpha\beta};\,\,\
\Gamma{{}^{\alpha}}{{}_{\mu}}{{}^{\beta}}=g^{\beta\nu}\Gamma^{\alpha}_{\mu\nu};\,\,\
\Sigma_{\alpha\beta}=\frac{i}{4}\Big{[}\gamma_{\alpha},\gamma_{\beta}\Big{]}$
$\displaystyle\\{\gamma^{\mu},\gamma^{\nu}\\}=2g^{\mu\nu}$ (13)

Since for radial trajectories only the $(r-t)$ sector of the metric (1) is
important, (2) can be expressed as

$\displaystyle
i\gamma^{\mu}\partial_{\mu}\psi-\frac{1}{2}\Big{(}g^{tt}\gamma^{\mu}\Gamma^{r}_{\mu
t}-g^{rr}\gamma^{\mu}\Gamma^{t}_{\mu r}\Big{)}\Sigma_{rt}\psi=0.$ (14)

Here the required nonvanishing connections are

$\displaystyle\Gamma^{r}_{tt}=\frac{f^{\prime}g}{2};\,\,\
\Gamma^{t}_{tr}=\frac{f^{\prime}}{2f}.$ (15)

Therefore, under the metric (1), the Dirac equation (2) reduces to

$\displaystyle
i\gamma^{t}\partial_{t}\psi+i\gamma^{r}\partial_{r}\psi+\frac{f^{\prime}g}{2f}\gamma^{t}\Sigma_{rt}\psi=0$
(16)

and the matrix form of $\Sigma_{rt}$ from (12) and (13) is given by

$\displaystyle\Sigma_{rt}=\frac{i}{2}\left(\begin{array}[]{c c c
c}0&0&i\sqrt{\frac{f}{g}}&0\\\ 0&0&0&-i\sqrt{\frac{f}{g}}\\\
-i\sqrt{\frac{f}{g}}&0&0&0\\\ 0&i\sqrt{\frac{f}{g}}&0&0\end{array}\right).$
(21)

To solve (16) we employ the following ansatz for the spin up (i.e. +ve
$r$-direction) and spin down (i.e. -ve $r$-direction) $\psi$ as

$\displaystyle\psi_{\uparrow}(t,r)=\left(\begin{array}[]{c}A(t,r)\\\ 0\\\
B(t,r)\\\
0\end{array}\right){\textrm{exp}}\Big{[}\frac{i}{\hbar}I_{\uparrow}(t,r)\Big{]}$
(26) $\displaystyle\psi_{\downarrow}(t,r)=\left(\begin{array}[]{c}0\\\
C(t,r)\\\ 0\\\ D(t,r)\\\
\end{array}\right){\textrm{exp}}\Big{[}\frac{i}{\hbar}I_{\downarrow}(t,r)\Big{]}$
(31)

where $I(r,t)$ is the one particle action which will be expanded in powers of
$\hbar$. Here we will only solve the spin up case explicitly since the spin
down case is fully analogous. Substituting the ansatz (26) in (16), we obtain
the following two equations:

$\displaystyle\Big{(}\frac{iA}{\sqrt{f}}\partial_{t}I_{\uparrow}+B\sqrt{g}\partial_{r}I_{\uparrow}\Big{)}+\hbar\Big{(}\frac{1}{\sqrt{f}}\partial_{t}A-i\sqrt{g}\partial_{r}B+i\frac{f^{\prime}\sqrt{g}}{4f}B\Big{)}=0$
(32)
$\displaystyle\Big{(}-\frac{iB}{\sqrt{f}}\partial_{t}I_{\uparrow}+A\sqrt{g}\partial_{r}I_{\uparrow}\Big{)}+\hbar\Big{(}-\frac{1}{\sqrt{f}}\partial_{t}B-i\sqrt{g}\partial_{r}A+i\frac{f^{\prime}\sqrt{g}}{4f}A\Big{)}=0.$
(33)

Since the last terms within the first bracket of the above equations do not
involve the single particle action, they will not contribute to the
thermodynamic entities of the black hole. Therefore we will drop these two
terms. Now taking $I_{\uparrow}=I$ and expanding $I,A$ and $B$ in powers of
$\hbar$, we find,

$\displaystyle I(r,t)=I_{0}(r,t)+\displaystyle\sum_{i}\hbar^{i}I_{i}(r,t)$
$\displaystyle
A=A_{0}+\displaystyle\sum_{i}\hbar^{i}A_{i};\,\,\,B=B_{0}+\displaystyle\sum_{i}\hbar^{i}B_{i}.$
(34)

where $i=1,2,3,......$. In these expansions the terms from ${\cal{O}}(\hbar)$
onwards are treated as quantum corrections over the semiclassical value
$I_{0}$, $A_{0}$ and $B_{0}$ respectively. Substituting (34) in (32) and (33)
and then equating the different powers of $\hbar$ on both sides, we obtain the
following two sets of equations:

$\displaystyle{\textrm{Set I}}:~{}\hbar^{0}:$
$\displaystyle\frac{i}{\sqrt{f}}A_{0}\partial_{t}I_{0}+\sqrt{g}B_{0}\partial_{r}I_{0}=0$
(35) $\displaystyle\hbar^{1}:$
$\displaystyle\frac{i}{\sqrt{f}}A_{0}\partial_{t}I_{1}+\frac{i}{\sqrt{f}}A_{1}\partial_{t}I_{0}+\sqrt{g}B_{0}\partial_{r}I_{1}+\sqrt{g}B_{1}\partial_{r}I_{0}=0$
(36) $\displaystyle\hbar^{2}:$
$\displaystyle\frac{i}{\sqrt{f}}A_{0}\partial_{t}I_{2}+\frac{i}{\sqrt{f}}A_{1}\partial_{t}I_{1}+\frac{i}{\sqrt{f}}A_{2}\partial_{t}I_{0}$
$\displaystyle+$
$\displaystyle\sqrt{g}B_{0}\partial_{r}I_{2}+\sqrt{g}B_{1}\partial_{r}I_{1}+\sqrt{g}B_{2}\partial_{r}I_{0}=0$
. . . and so on. $\displaystyle{\textrm{Set II}}:~{}\hbar^{0}:$
$\displaystyle-\frac{i}{\sqrt{f}}B_{0}\partial_{t}I_{0}+\sqrt{g}A_{0}\partial_{r}I_{0}=0$
(38) $\displaystyle\hbar^{1}:$
$\displaystyle-\frac{i}{\sqrt{f}}B_{0}\partial_{t}I_{1}-\frac{i}{\sqrt{f}}B_{1}\partial_{t}I_{0}+\sqrt{g}A_{0}\partial_{r}I_{1}+\sqrt{g}A_{1}\partial_{r}I_{0}=0$
(39) $\displaystyle\hbar^{2}:$
$\displaystyle-\frac{i}{\sqrt{f}}B_{0}\partial_{t}I_{2}-\frac{i}{\sqrt{f}}B_{1}\partial_{t}I_{1}-\frac{i}{\sqrt{f}}B_{2}\partial_{t}I_{0}$
$\displaystyle+$
$\displaystyle\sqrt{g}A_{0}\partial_{r}I_{2}+\sqrt{g}A_{1}\partial_{r}I_{1}+\sqrt{g}A_{2}\partial_{r}I_{0}=0$
. . . and so on.

Equations (35) and (38) are collectively known as the semiclassical Hamilton-
Jacobi equations for a Dirac particle. Since the metric (1) is stationary it
has timelike Killing vectors. Thus we will look for solutions of (35) and (38)
which behave as

$\displaystyle I_{0}=\omega t+W(r),$ (41)

where $\omega$ is the energy of the particle. Substituting this in (35) and
(38) we obtain,

$\displaystyle\frac{iA_{0}}{\sqrt{f}}\omega+B_{0}\sqrt{g}W^{\prime}(r)=0$
$\displaystyle-\frac{iB_{0}}{\sqrt{f}}\omega+A_{0}\sqrt{g}W^{\prime}(r)=0.$
(42)

These two equations have two possible solutions:

$\displaystyle A_{0}=iB_{0};\,\,\
W_{+}(r)=\omega\int_{0}^{r}\frac{dr}{\sqrt{f(r)g(r)}}$ $\displaystyle
A_{0}=-iB_{0};\,\,\,W_{-}(r)=-\omega\int_{0}^{r}\frac{dr}{\sqrt{f(r)g(r)}}$
(43)

where $W_{+}$($W_{-}$) corresponds to ingoing (outgoing) solutions. The limits
of the integration are chosen such that the particle goes through the horizon
$r=r_{H}$. Therefore the solution for $I_{0}(r,t)$ is

$\displaystyle I_{0}(r,t)=\omega
t\pm\omega\int_{0}^{r}\frac{dr}{\sqrt{f(r)g(r)}}.$ (44)

Now, it is interesting to note that using (43) and (44) in the equations of
Set I and Set II simultaneously and then solving we get relations connecting
different orders in the expansion of $A$ with those of $B$:

$\displaystyle A_{a}=\pm iB_{a}$ (45)

where $a=0,1,2,3,....$. These lead to a simplified form of all the equations
in Set I and Set II as,

$\displaystyle\partial_{t}I_{a}=\pm\sqrt{fg}\partial_{r}I_{a}$ (46)

i.e. the functional form of the above individual linear differential equations
is same and is identical to the usual semiclassical Hamilton-Jacobi equations
(35) and (38). Therefore the solutions of these equations are not independent
and $I_{i}$’s are proportional to $I_{0}$. A similar situation happened for
scalar particle tunneling [12]. Since $I_{0}$ has the dimension of $\hbar$ the
proportionality constants should have the dimension of inverse of $\hbar^{i}$.
Again in the units $G=c=k_{B}=1$ the Planck constant $\hbar$ is of the order
of square of the Planck Mass $M_{P}$ and so from dimensional analysis the
proportionality constants have the dimension of $M^{-2i}$ where $M$ is the
mass of black hole. Specifically, for Schwarzschild type black holes having
mass as the only macroscopic parameter, these considerations show that the
most general expression for $I$, following from (34), valid for (46), is given
by,

$\displaystyle
I(r,t)=\Big{(}1+\sum_{i}\beta_{i}\frac{\hbar^{i}}{M^{2i}}\Big{)}I_{0}(r,t).$
(47)

where $\beta_{i}$’s are dimensionless constant parameters.

The above analysis shows that to obtain a solution for $I(r,t)$ it is
therefore enough to solve for $I_{0}(r,t)$ which has the solution of the form
(44). In fact the standard Hamilton-Jacobi solution determined by this
$I_{0}(r,t)$ is just modified by a prefactor to yield the complete solution
for $I(r,t)$. Substituting (44) in (47) we obtain

$\displaystyle
I(r,t)=\Big{(}1+\sum_{i}\beta_{i}\frac{\hbar^{i}}{M^{2i}}\Big{)}\Big{[}\omega
t\pm\omega\int_{0}^{r}\frac{dr}{\sqrt{f(r)g(r)}}\Big{]}.$ (48)

Therefore the ingoing and outgoing solutions of the Dirac equation (2) under
the background metric (1) is given by exploiting (26) and (48),

$\displaystyle\psi_{{\textrm{in}}}\sim{\textrm{exp}}\Big{[}-\frac{i}{\hbar}(1+\sum_{i}\beta_{i}\frac{\hbar^{i}}{M^{2i}})\Big{(}\omega
t+\omega\int_{0}^{r}\frac{dr}{\sqrt{f(r)g(r)}}\Big{)}\Big{]}$ (49)

and

$\displaystyle\psi_{{\textrm{out}}}\sim{\textrm{exp}}\Big{[}-\frac{i}{\hbar}(1+\sum_{i}\beta_{i}\frac{\hbar^{i}}{M^{2i}})\Big{(}\omega
t-\omega\int_{0}^{r}\frac{dr}{\sqrt{f(r)g(r)}}\Big{)}\Big{]}.$ (50)

Now for the tunneling of a particle across the horizon the nature of the
coordinates change. The sign of the metric coefficients in the $(r-t)$ sector
is altered. This indicates that ‘$t$’ coordinate has an imaginary part for the
crossing of the horizon of the black hole and correspondingly there will be a
temporal contribution to the probabilities for the ingoing and outgoing
particles. This has similarity with[23] where they show for the Schwarzschild
metric that two patches across the horizon are connected by a discrete
imaginary amount of time.

The ingoing and outgoing probabilities of the particle are, therefore, given
by,

$\displaystyle
P_{{\textrm{in}}}=|\psi_{{\textrm{in}}}|^{2}\sim{\textrm{exp}}\Big{[}\frac{2}{\hbar}(1+\sum_{i}\beta_{i}\frac{\hbar^{i}}{M^{2i}})\Big{(}\omega{\textrm{Im}}~{}t+\omega{\textrm{Im}}\int_{0}^{r}\frac{dr}{\sqrt{f(r)g(r)}}\Big{)}\Big{]}$
(51)

and

$\displaystyle
P_{{\textrm{out}}}=|\psi_{{\textrm{out}}}|^{2}\sim{\textrm{exp}}\Big{[}\frac{2}{\hbar}(1+\sum_{i}\beta_{i}\frac{\hbar^{i}}{M^{2i}})\Big{(}\omega{\textrm{Im}}~{}t-\omega{\textrm{Im}}\int_{0}^{r}\frac{dr}{\sqrt{f(r)g(r)}}\Big{)}\Big{]}$
(52)

Now the ingoing probability $P_{\textrm{in}}$ has to be unity in the classical
limit (i.e. $\hbar\rightarrow 0$) - when there is no reflection and everything
is absorbed - instead of zero or infinity [12].Thus, in the classical limit,
(51) leads to,

$\displaystyle{\textrm{Im}}~{}t=-{\textrm{Im}}\int_{0}^{r}\frac{dr}{\sqrt{f(r)g(r)}}.$
(53)

From the above one can easily show that ${\textrm{Im}}~{}t=-2\pi M$ for the
Schwarzschild spacetime which is precisely the imaginary part of the
transformation $t\rightarrow t-2i\pi M$ when one connects the two regions
across the horizon as shown in [23]. Therefore the probability of the outgoing
particle is

$\displaystyle
P_{{\textrm{out}}}\sim{\textrm{exp}}\Big{[}-\frac{4}{\hbar}\omega\Big{(}1+\sum_{i}\beta_{i}\frac{\hbar^{i}}{M^{2i}}\Big{)}{\textrm{Im}}\int_{0}^{r}\frac{dr}{\sqrt{f(r)g(r)}}\Big{]}.$
(54)

Now using the principle of “detailed balance” [2, 9]

$\displaystyle
P_{{\textrm{out}}}={\textrm{exp}}\Big{(}-\frac{\omega}{T_{h}}\Big{)}P_{\textrm{in}}={\textrm{exp}}\Big{(}-\frac{\omega}{T_{h}}\Big{)}$
(55)

we obtain the temperature of the black hole as

$\displaystyle
T_{h}=T_{H}\Big{(}1+\sum_{i}\beta_{i}\frac{\hbar^{i}}{M^{2i}}\Big{)}^{-1}$
(56)

where

$\displaystyle
T_{H}=\frac{\hbar}{4}\Big{(}{\textrm{Im}}\int_{0}^{r}\frac{dr}{\sqrt{f(r)g(r)}}\Big{)}^{-1}$
(57)

is the standard semiclassical Hawking temperature of the black hole and other
terms are the corrections due to the quantum effect. Using this expression and
knowing the metric coefficients $f(r)$ and $g(r)$ one can easily find out the
temperature of the corresponding black hole. The same result was also obtained
in [12] for scalar particle tunneling.

For the Schwarzschild black hole the metric coefficients are

$\displaystyle f(r)=g(r)=(1-\frac{r_{H}}{r});\,\,\,r_{H}=2M.$ (58)

Therefore using (56) and (57) it is easy to write the corrected Hawking
temperature:

$\displaystyle T_{h}=\frac{\hbar}{8\pi
M}\Big{(}1+\sum_{i}\beta_{i}\frac{\hbar}{M^{2i}}\Big{)}^{-1}.$ (59)

Now use of the Gibbs form of first law of thermodynamics gives the corrected
form of the Bekenstein-Hawking entropy:

$\displaystyle S_{\textrm{bh}}$ $\displaystyle=$
$\displaystyle\int\frac{dM}{T_{h}}=\frac{4\pi M^{2}}{\hbar}+8\pi\beta_{1}\ln
M-\frac{4\pi\hbar\beta_{2}}{M^{2}}+{\textrm{higher order terms in $\hbar$}}$
(60) $\displaystyle=$ $\displaystyle\frac{\pi
r_{H}^{2}}{\hbar}+8\pi\beta_{1}\ln
r_{H}-\frac{16\pi\hbar\beta_{2}}{r_{H}^{2}}+{\textrm{higher order terms in
$\hbar$}}$

The area of the event horizon is

$\displaystyle A=4\pi r_{H}^{2}$ (61)

so that,

$\displaystyle S_{\textrm{bh}}$ $\displaystyle=$
$\displaystyle\frac{A}{4\hbar}+4\pi\beta_{1}\ln
A-\frac{64\pi^{2}\hbar\beta_{2}}{A}+......................$ (62)

It is noted that the first term is the usual semiclassical contribution to the
area law $S_{\textrm{BH}}=\frac{A}{4\hbar}$ [24, 25]. The other terms are the
quantum corrections. Now it is possible to express the quantum corrections in
terms of $S_{\textrm{BH}}$ by eliminating $A$:

$\displaystyle S_{\textrm{bh}}=S_{\textrm{BH}}+4\pi\beta_{1}\ln
S_{\textrm{BH}}-\frac{16\pi^{2}\beta_{2}}{S_{\textrm{BH}}}+.........$ (63)

Interestingly the leading order correction is logarithmic in $A$ or
$S_{\textrm{BH}}$ which was found earlier in [13, 14] by field theory
calculations and later in [15, 20] by quantum geometry method. The higher
order corrections involve inverse powers of $A$ or $S_{\textrm{BH}}$.

To determine the value of the coefficients $\beta_{1},\beta_{2}$ etc we will
adopt the following steps. The point is that nonzero values for these
coefficients are related to quantum corrections (loop effects). Such
corrections, in a field theoretical approach, are manifested by the presence
of anomalies. Now it is a well known fact that it is not possible to
simultaneously preserve general coordinate (diffeomorphism) invariance and
conformal invariance. Retaining general coordinate invariance, one finds the
breakdown of conformal invariance leading to the presence of nonvanishing
trace of the stress tensor. We now show that the coefficients appearing in
(63) are related to this trace anomaly.

We begin by studying the behaviour of the action (47) upto order $\hbar^{2}$
under an infinitesimal constant scale transformation, parametrised by $k$, of
the metric coefficients,

$\displaystyle\tilde{g}{{}_{\mu\nu}}=kg_{\mu\nu}\simeq(1+\delta k)g_{\mu\nu}.$
(64)

Under this the metric coefficients of (1) change as
$\tilde{f}=kf,\tilde{g}=k^{-1}g$. Also, in order to preserve the scale
invariance of the Dirac equation (2), the field $\psi$ should transform as
$\tilde{\psi}=k^{\frac{1}{2}}\psi$. On the other hand, $\psi$ has the
dimension of ${(\textrm{mass})}^{\frac{3}{2}}$ and since in our case the only
mass parameter is the black hole mass $M$, the infinitesimal change of it is
given by,

$\displaystyle\tilde{M}=k^{\frac{1}{3}}M\simeq(1+\frac{1}{3}\delta k)M.$ (65)

Now from (54) the imaginary part of the semiclassical contribution of the
single particle action is

$\displaystyle\textrm{Im}I{{}_{0}}{{}_{(\textrm{out})}}=-2\omega{\textrm{Im}}~{}\int_{0}^{r}\frac{dr}{\sqrt{f(r)g(r)}}$
(66)

where $\omega$ gets identified with the energy (i.e. mass $M$) of a stable
black hole [8]. Therefore $\omega$ and $\hbar$ transforms like (65) and
$M^{2}$ respectively under (64).

Considering only the $\hbar$ and $\hbar^{2}$ order terms in (47) and using
(65) we obtain, under the scale transformation,

$\displaystyle{\tilde{I}}_{(1+2)}$ $\displaystyle\equiv$
$\displaystyle\hbar\textrm{Im}\tilde{I}{{}_{1}}{{}_{(\textrm{out})}}+\hbar^{2}\textrm{Im}\tilde{I}{{}_{2}}{{}_{(\textrm{out})}}$
(67) $\displaystyle=$
$\displaystyle\Big{(}\frac{{\tilde{\hbar}}\beta_{1}}{\tilde{M}^{2}}+\frac{{\tilde{\hbar}}^{2}\beta_{2}}{\tilde{M}^{4}}\Big{)}\textrm{Im}\tilde{I}{{}_{0}}{{}_{(\textrm{out})}}$
$\displaystyle\simeq$
$\displaystyle\Big{(}\frac{\hbar\beta_{1}}{{M}^{2}}+\frac{{{\hbar}}^{2}\beta_{2}}{{M}^{4}}\Big{)}(1+\frac{1}{3}\delta
k)\textrm{Im}I{{}_{0}}{{}_{(\textrm{out})}}$ $\displaystyle=$ $\displaystyle
I_{(1+2)}+\Big{(}\frac{\hbar\beta_{1}}{{M}^{2}}+\frac{{{\hbar}}^{2}\beta_{2}}{{M}^{4}}\Big{)}\frac{1}{3}\delta
k\textrm{Im}I{{}_{0}}{{}_{(\textrm{out})}}$

Therefore

$\displaystyle\delta I_{(1+2)}$ $\displaystyle=$
$\displaystyle{\tilde{I}}_{(1+2)}-I_{(1+2)}$ (68) $\displaystyle\simeq$
$\displaystyle\Big{(}\frac{\hbar\beta_{1}}{{M}^{2}}+\frac{{{\hbar}}^{2}\beta_{2}}{{M}^{4}}\Big{)}\frac{1}{3}\delta
k\textrm{Im}I{{}_{0}}{{}_{(\textrm{out})}}$

leading to,

$\displaystyle\frac{\delta I_{(1+2)}}{\delta
k}=\Big{(}\frac{\hbar\beta_{1}}{{M}^{2}}+\frac{{{\hbar}}^{2}\beta_{2}}{{M}^{4}}\Big{)}\frac{1}{3}\textrm{Im}I{{}_{0}}{{}_{(\textrm{out})}}$
(69)

Now use of the definition of the energy-momentum tensor and (69) yields,

$\displaystyle\textrm{Im}\int d^{4}x\sqrt{-g}T_{\mu}^{\mu}=\frac{2\delta
I_{(1+2)}}{\delta
k}=\Big{(}\frac{\hbar\beta_{1}}{{M}^{2}}+\frac{{{\hbar}}^{2}\beta_{2}}{{M}^{4}}\Big{)}\frac{2}{3}\textrm{Im}I{{}_{0}}{{}_{(\textrm{out})}}$
(70)

Thus, in the presence of a trace anomaly, the action is not invariant under
the scale transformation. Since for the Schwarzschild black hole $f(r)$ and
$g(r)$ are given by (58), from (66) we obtain
$\textrm{Im}S_{0}^{(\textrm{out})}=-4\pi\omega M$. Substituting this in (70)
we obtain for $\omega=M$ as,

$\displaystyle\hbar\beta_{1}+\frac{\hbar^{2}\beta_{2}}{M^{2}}=-\frac{3}{8\pi}\textrm{Im}\int
d^{4}x\sqrt{-g}T_{\mu}^{\mu}$ (71)

where $T_{\mu}^{\mu}$ is calculated upto two loops. Starting from the action
(47) and following the identical steps as above, a similer relation among all
$\beta$’s with the right hand side of (71) can be established. In this case
$T_{\mu}^{\mu}$ is due to all loop expansions.

Since the higher loop calculations to get $T_{\mu\nu}$ (from which
$T_{\mu}^{\mu}$ is obtained) is very much complicated, usually in literature
[28] only one loop calculation for $T_{\mu\nu}$ is discussed. Thus, comparing
only the $\hbar^{1}$ order on both sides of (71), we obtain,

$\displaystyle\beta_{1}=-\frac{3}{8\pi}{\textrm{Im}}\int
d^{4}x\sqrt{-g}{T^{\mu}_{\mu}}^{(1)}$ (72)

This relation clearly shows that $\beta_{1}$ is connected to the trace
anomaly. A similar relation is given in [26] where it has been shown that the
coefficient $\beta_{1}$ is related to trace anomaly for the scalar particle
tunneling. The only difference is the factor before the integration. This
agrees well with the earlier conclusion [27, 13] where using conformal field
theory technique it was shown that this $\beta_{1}$ is related to trace
anomaly and is given by,

$\displaystyle\beta_{1}=-\frac{1}{360\pi}\Big{(}-N_{0}-\frac{7}{4}N_{\frac{1}{2}}+13N_{1}+\frac{233}{4}N_{\frac{3}{2}}-212N_{2}\Big{)}$
(73)

‘$N_{s}$’ denotes the number of fields with spin ‘$s$’. In our case
$N_{\frac{1}{2}}=1$ and $N_{0}=N_{1}=N_{\frac{3}{2}}=N_{2}=0$.

To conclude, we have successfully extended our approach [12] of scalar
particle tunneling beyond semiclassical approximation to the model of fermion
tunneling. We have considered all orders in the single particle action for
fermion tunneling through the event horizon of the black hole. We showed that
higher order correction terms of the action are proportional to the
semiclassical contribution. A similar result was shown earlier in [12] for the
scalar particle tunneling. By dimensional argument and principle of “detailed
balance” the same form of the modified Hawking temperature, as in the scalar
case, was recovered. The logarithmic and inverse powers of area corrections to
the Bekenstein-Hawking area law were also reproduced. Finally, we showed that
the coefficient of the logarithmic term of entropy is related to trace
anomaly. However, the prefactor appearing in this term is different from the
scalar particle example, a result that is supported by earlier works [27, 13].

Here we have only told about $\beta_{1}$. Discussions on other coefficients
can also be given from (71), but since no information about $T_{\mu\nu}$ due
to multi-loops is available in the literature, we cannot say anything about
them at this moment.

Acknowledgment:
I wish to thank Prof. Rabin Banerjee for suggesting this investigation and
constant encouragement.

## References

  * [1] M.K.Parikh and F.Wilczek, Phys. Rev. Lett. 85, 5042 (2000) [arXiv:hep-th/9907001].
M.K.Parikh, Int. J. Mod. Phys. D 13, 2351 (2004) [arXiv:hep-th/0405160].

  * [2] K.Srinivasan and T.Padmanabhan, Phys. Rev. D 60, 024007 (1999) [arXiv:gr-qc/9812028].
S. Shankaranarayanan, K. Srinivasan and T. Padmanabhan, Mod. Phys. Lett. A 16,
571 (2001) [arXiv:gr-qc/0007022].
S.Shankaranarayanan, T.Padmanabhan and K.Srinivasan, Class. Quantum Grav. 19,
2671 (2002) [arXiv:gr-qc/0010042].
S.Shankaranarayanan, Phys. Rev. D 67, 084026 (2003) [arXiv:gr-qc/0301090].

  * [3] V.A.Berezin, A.M.Boyarsky and A.Yu.Neronov, Gravitation and Cosmology 5, 16 (1999) [arXiv:gr-qc/0605099].
  * [4] M.Arzano, A.J.M.Medved and E.C.Vagenas, JHEP 0509, 037 (2005) [arXiv:hep-th/0505266].
A.J.M.Medved and E.C.Vagenas, Mod. Phys. Lett. A 20, 2449 (2005) [arXiv:gr-
qc/0504113].

  * [5] Qing-Quan Jiang, Shuang-Qing Wu and Xu Cai, Phys. Rev. D 73 064003 (2006) [arXiv:hep-th/0512351].
Yapeng Hu, Jingyi Zhang and Zheng Zhao, Mod. Phys. Lett. A 21 2143 (2006)
[arXiv:gr-qc/0611026].
Zhibo Xu and Bin Chen, Phys. Rev. D 75 024041 (2007) [arXiv:hep-th/0612261].
Cheng-Zhou Liu and Jian-Yang Zhu, “Hawking radiation as tunneling from
Gravity’s rainbow”, [arXiv:gr-qc/0703055].
Qing-Quan Jiang and Shuang-Qing Wu, Phys. Lett. B 635, 151 (2006) [arXiv:hep-
th/0511123].
Yapeng Hu, Jingyi Zhang and Zheng Zhao, Int.J.Mod.Phys. D 16, 847 (2007)
[arXiv:gr-qc/0611085].
G.E.Volovik, “On de Sitter radiation via quantum tunneling”,
[arXiv:0803.3367].
Jingyi Zhang, “Black hole quantum tunnelling and black hole entropy
correction” [arXiv:0806.2441].

  * [6] E.T.Akhmedov, V.Akhmedova and D.Singleton, Phys. Lett. B 642, 124 (2006) [arXiv:hep-th/0608098].
E.T.Akhmedov, V.Akhmedova, D.Singleton and T.Pilling, Int.J.Mod.Phys. A 22,
1705 (2007) [arXiv:hep-th/0605137].
T.Pilling, Phys. Lett. B 660, 402 (2008) [arXiv:0709.1624].
T.K.Nakamura, “Factor two discrepancy of Hawking radiation temperature”,
[arXiv:0706.2916].
B.D.Chowdhury, Pramana 70, 593 (2008) [arXiv:hep-th/0605197].

  * [7] M.Angheben, M.Nadalini, L.Vanzo and S.Zerbini, JHEP 0505, 014 (2005) [arXiv:hep-th/0503081].
R.Kerner and R.B.Mann, Phys. Rev. D 73, 104010 (2006) [arXiv:gr-qc/0603019].
P.Mitra, Phys. Lett. B 648, 240 (2007) [arXiv:hep-th/0611265].

  * [8] R.Banerjee and B.R.Majhi, Phys. Lett. B 662, 62 (2008) [arXiv:0801.0200].
  * [9] R.Banerjee, B.R.Majhi and S.Samanta, Phys. Rev. D 77 124035 (2008) [arXiv:0801.3583].
  * [10] R.Banerjee, B.R.Majhi and S.K.Modak, “Area Law in Noncommutative Schwarzschild Black Hole”, [arXiv:0802.2176].
  * [11] R.Kerner and R.B.Mann, Class. Quant. Grav. 25, 095014 (2008) [arXiv:0710.0612].
Ran Li and Ji-Rong Re, Phys. Lett. B 661 370 (2008) [arXiv:0802.3954].
R.D.Criscienzo and L.Vanzo, “Fermion Tunneling from Dynamical Horizons”,
[arXiv:0803.0435].
Ran Li and Ji-Rong Ren, Class. Quant. Grav. 25 125016 (2008)
[arXiv:0803.1410].
R.Kerner and R.B.Mann, Phys. Lett. B 665 277 (2008) [arXiv:0803.2246].
De-You Chen, Qing-Quan Jiang, Shu-Zheng Yang and Xiao-Tao Zu, “Fermions
tunneling from the charged dilatonic black holes”, [arXiv:0803.3248].
D.Y.Chen, Q.Q.Jiang and X.T.Zu, Phys. Lett. B 665 106 (2008)
[arXiv:0804.0131].
Shiwei Zhou, Wenbiao Liu, Phys. Rev. D 77 104021 (2008).
Q.Q.Jiang, Phys.Rev. D 78 044009, (2008).

  * [12] R.Banerjee and B.R.Majhi, JHEP 0806 095 (2008) [arXiv:0805.2220].
  * [13] D.V.Fursaev, Phys. Rev. D 51, R5352 (1995) [arXiv:hep-th/9412161].
  * [14] R.B.Mann and S.N.Solodukhin, Nucl. Phys. B 523, 293 (1998) [arXiv:hep-th/9709064].
  * [15] R.K.Kaul and P.Majumdar, Phys. Rev. Lett. 84, 5255 (2000) [arXiv:gr-qc/0002040].
  * [16] T.R.Govindarajan, R.K.Kaul and V.Suneeta, Class. Quantum Grav. 18, 2877 (2001) [arXiv:gr-qc/0104010].
  * [17] S.Das, P.Majumdar and R.K.Bhaduri, Class. Quantum Grav. 19, 2355 (2002) [arXiv:hep-th/0111001].
  * [18] S.S.More, Class. Quantum Grav. 22, 4129 (2005) [gr-qc/0410071].
  * [19] S.Mukherji and S.S.Pal, JHEP 0205, 026 (2002) [arXiv:hep-th/0205164].
  * [20] A.Ghosh and P.Mitra, Phys. Lett. B 616, 114 (2005) [arXiv:gr-qc/0411035].
  * [21] For a review and a complete list of papers on logarithmic corrections, see D.N.Page, New Journal of Phys. 7, 203 (2005) [arXiv:hep-th/0409024].
  * [22] S.K.Modak, “Corrected entropy of BTZ black hole in tunneling approach”, [arXiv:0807.0959].
  * [23] E.T.Akhmedov, T.Pilling and D.Singleton, “Subtleties in the quasi-classical calculation of Hawking radiation”, [arXiv:0805.2653].
  * [24] J.D.Bekenstein, PhD Thesis Princeton University, Princeton, NJ (1972).
J.D.Bekenstein, Lett. Nuovo Cimento 4 737 (1972).
J.D.Bekenstein, Phys. Rev. D 7, 2333 (1973).
J.D.Bekenstein, Phys. Rev. D 9 3292 (1974).

  * [25] J.M.Bardeen, B.Carter and S.W.Hawking, Commun. Math. Phys. 31, 161 (1973).
  * [26] R.Banerjee and B.R.Majhi, “Quantum Tunneling and Trace Anomaly”, [arXiv:0808.3688].
  * [27] S.M.Christensen and M.J.Duff, Phys. Lett. B 76 571 (1978).
  * [28] B.S.DeWitt, Phys. Rep. 19, 295 (1975).