# Scaling of Entanglement-Assisted Communication in Amplified Fiber Links

Sekavčnik, Simon and Nötzel Janis

(November 2022)

###### Abstract

Quantum communication technology offers several advanced strategies. However,
their practical use is often times still not well understood. In this work we
outline the theoretical communication capacity scaling in amplified fiber
links. We present a scenario in which the assistance via pre-shared
entanglement offers an arbitrary capacity, given enough bandwidth and spatial
modes are provided by the fiber. The future capacity advantage over
conventional classical techniques as well as non-assisted quantum techniques
is potentially infinite. We discuss this theoretical observation in connection
to current trends in fiber development.

##  I Introduction

The ever growing complexities of present communication networks are a true
testament to human ingenuity. Our reliance on the global network is growing by
the day. With the amount of network traffic poised to increase, coupled with
high demand for ultra low-latency communication and limit in transmission
power allowed by fiber, we are presented with a challenge to reduce the energy
consumption while accommodating the required needs.

Wast majority of network traffic is serviced by optical fiber links. With
current classical technologies the communication capacity is limited by the
Shannon capacity. The eventual emergence of Quantum Information Processing
(QIP) technologies will - in addition to providing new, never seen before
applications and communication protocols - move the theoretical communication
capacity bounds.

The fact that capacity bounds are higher for QIP enabled communication
channels is well documented [1, 2]. Finding use-cases where effects would
translate into significant practical communication advantage is a difficult
task. Especially difficult is fitting entanglement assisted communication into
a real world communication problem. In [3, 4] the authors examined the use of
QIP for the transceiver design of optical fiber networks and pointed out the
possibility to reduce energy consumption of the optical amplifiers. In [3]
different capacity limits of amplified optical fiber limits were calculated as
well, which accounted for different relations between transmission power and
amplifier gains. The work [5] studied the performance under more realistic
conditions, taking into account the Kerr-induced phase noise based on the
model [6] and a specific design of the transceiver based on the work [7].
Later [8] tried to showcase the advantages offered by entanglement-assisted
QIP in high-density, ultra low latency networks. The advantages were clearly
laid out, however, it was also pointed out that the parameter region required
for meaningful entanglement-assisted QIP advantages are difficult to arrive at
and hard to argue for.

In this work we focus on outlining practical communication scenario in which
the use of QIP, particularly with entanglement assistance, translates into
capacity advantage over purely classical and non-assisted QIP method. In order
to showcase the different capacity bounds we consider an optical fiber link
with three different transmitter and receiver types.

The differences in the above three capacities stem from the fact that when the
product $\tau P/M$ is small enough (for example below one), the QIP enabled
system will display an advantage that scales logarithmically with the $\tau
P/M$. Additionally the entanglement assisted links will display an advantage
over non-assisted QIP links if there is a sufficient amount of noise in the
system $\nu\gg 1$, combined with transmitters operating with low power per
pulse $P/M\ll 1$.

The questions we are concerned with here are two-fold: First, what the
technological meaning of these relations is today, and second, how they might
influence future communication system design. Therefore we root our modelling
into the current fiber transmission reality by assuming communication at
$1550nm$. Here, the fiber nonlinearities dictate a power limit of
approximately $100mW$ per fiber to current system design, which leads to a
value of around $10^{18}$ photons being emitted at the sender per second. We
assume this energy is spread over a wide bandwidth $B$ by utilizing several
orthogonal channels. The most extreme number $B$ one could think of (based on
fibers available today) would be in the order of $B=10^{13}$, leaving us with
$10^{5}$ photons per slot. It is thus difficult to motivate the operating
regime $P/M\ll 1$ based solely on setting $M=B$. However, the recent
literature points out the possibility to utilizing multimode fiber to transmit
several spatially orthogonal channels over the same fiber, and a number of
$N=10^{3}$ of such channels has already been demonstrated as technically
feasible today [9]. Taking such possibilities into account we may set
$M=N\cdot B$ to argue for the very low power per pulse needed to see
advantages of entanglement-assisted communication. Since optical amplification
is the standard technology enabling high throughput data transmission over
more than $100km$, and since optical amplification not only regenerates the
signal but inevitably also adds noise to it, it is therefore possible to
motivate communication in a setting where both $\nu\gg 1$ and $P/M\ll 1$. Our
analysis clarifies that the condition $\nu\gg 1$ might be relaxed to $\nu>0$
without changing the general scaling of entanglement-assisted capacity with
the value $M$ when amplifiers are used. In our model we assume an amplifier
setting which _completely_ regenerates the signal. This is in contrast to the
previous work [4] and more in line with the analysis of the quantum limit as
in [3], where it has been shown that under the rule of constant amplification
(in the sense of being independent of the signal energy), both Shannon- and
Holevo capacity approach a limiting value when $M\to\infty$. Our results show
that, under the same amplifier rule, the entanglement-assisted capacity can
still grow unlimited when $M\to\infty$.

For the interpretation of our results it is important to note: the key
assumption to the analysis is that future fiber networks will need to obey a
per-fiber power limit $P$, while at the same time being able to transport an
increasing number of orthogonal channels, so that the assumption $M\to\infty$
is justified.

##  II System model

Following the exposition of [10] we use a modelling approach here which starts
with a narrowband, linearly polarized optical signal coming in the form of
uniformly spaced pulses in temporal slots of duration $B^{-1}$. Due to the
modulation being independently applied to these slots, the Fourier transform
shows that such a system occupies a spectral bandwidth proportional to the
slot rate $B$. Such a system can quantum-mechanically be modelled by assigning
a quantum state $|\alpha\rangle\in\mathcal{F}$ to every slot; $\mathcal{F}$
being the Fock space which has a countably-infinite orthonormal basis
$\\{|n\rangle\\}_{n=0}^{\infty}$ called the number-state basis. A system in
state $|n\rangle$ consists of exactly $n$ photons, and these number states can
be used to express a system in a coherent state $|\alpha\rangle$ as

$\displaystyle|\alpha\rangle=\exp^{-|\alpha|^{2}/2}\sum_{n=0}^{\infty}\frac{\alpha^{n}}{\sqrt{n!}}|n\rangle.$
(1)

The expected energy of a coherent state $|\alpha\rangle$ is then equal to
$|\alpha|^{2}$. If these pulses are detected and then converted to electronic
values in a slot-by-slot fashion, and if thermal noise with an average number
of $\nu$ photons per slot is present at the receiver, one recovers the well-
known Shannon formula for the capacity of the channel:

$\displaystyle C_{S}(M,P,\tau,\nu)$
$\displaystyle=M\cdot\log{\Big{(}1+\frac{\tau\cdot P}{M(1+\nu)}\Big{)}}$ (2)

Motivated by the work [11, 12] of Holevo, researchers have also studied the
capacity of such an optical data transmission system when arbitrary (quantum)
measurements are allowed at the receiver. When it comes to concrete
realization of this concept, one speaks of the so-called Joint Detection
Receiver (JDR), where a first proposal was made in [7] and a second one in
[13]. Both proposals use an interferometric setup to interfere a series of
consecutive pulses (slots) with each other prior to detection. This setup is
ideal to harness the advantage of JDRs, which is in the region where the
received signal energy _per slot_ is very weak. Due to the relation

$\displaystyle P=\mathbb{E}(|\alpha_{j}|^{2})Bhf_{c}$ (3)

where $P$ is the average optical power carried by the signal, $h=6.626\cdot
10^{-34}J\cdot s$, $f_{c}$ is the carrier frequency and the complex numbers
$\alpha_{j}$ specify the symbols that are generated at the sender, it is
possible to reach the limit of small received energy per slot by keeping $P$
constant while increasing $B$. The ultimate capacity for a communication
system employing a JDR is given by the Holevo capacity,

$\displaystyle C_{J}(M,P,\tau,\nu)$ $\displaystyle=M\cdot\Big{(}g(\tau\cdot
P\cdot M^{-1}+\nu)-g(\nu)\Big{)},$ (4)

where $g(x)=(x+1)\log(x+1)-x\log x$. Note that we assumed a noise level $\nu$
which is _independent_ of the power $P/M$ of the individual pulses both in Eq.
(2) and in Eq. (4). This assumption is actually the method by which a
separation between the data transmission concepts outlined so far on the one
side and entanglement-assisted data transmission on the other can be achieved
when $M\to\infty$.

Finally, let us introduce the third major data transmission concept that
emerged from the work [14]. It is the concept of entanglement-assisted
communication. Here, it is assumed that at the point in time when classical
information is encoded into a quantum state at the sender side, this quantum
state is not a simple coherent state $|\alpha\rangle$ but rather an arbitrary
quantum state shared between the sender and receiver. This implies the act of
sharing said entangled state is to be taken into account in the system design.
The work [15] described such data transmission system as one where idle
periods are used to execute the task of sharing entangled states, which is
then consumed when data needs to be transmitted by an application. The states
that can be utilized for this task can for example be so-called squeezed
states

$\displaystyle|s_{r}\rangle=\frac{1}{\cosh
r}\sum_{n=0}^{\infty}\tanh(r)^{n}|n\rangle_{I}|n\rangle_{D},$ (5)

where the index $I$ stands for the idler mode which is sent prior to data
transmission, and the index $D$ for the data mode which is used to encode the
data when it arrives at the sender and which is then transmitted to the
receiver. The capacity of the above described system when the idler modes are
shared and it is used in entanglement-assisted mode is given by

$\displaystyle C_{E}(M,P,\tau,\nu)$
$\displaystyle=M\sum_{x=0}^{1}g(\tfrac{\tau
P}{M}+x\nu)-g(d_{x}(\tau,\tfrac{P}{M},\nu)).$ (6)

Here, the functions $d_{x}$ are defined via

$\displaystyle d_{x}(\tau,n,\nu)$
$\displaystyle=(d(\tau,n,\nu)-1+(-1)^{x}((\tau-1)n+\nu))/2$ (7) $\displaystyle
d(\tau,n,\nu)$ $\displaystyle=\sqrt{((1+\tau)n+\nu+1)^{2}-4\tau n(n+1)}.$ (8)

We choose to divide a link into $K$ segments of equal length $L$, as presented
in Figure 1. At the end of each segment the signal is amplified by amplifier.
At this level of analysis the amplifier is completely modeled by it’s
amplification factor $G$. The transmittivity of one segment can thus be
computed by $\tau_{L}=e^{-\alpha\cdot L}$, where for numerical examples we use
the value $\alpha=0.05$ which models the transmittivity of SMF-28 fiber [6].

Figure 1: Visualization of amplified segmented link without pre-amplification
prior to detection (above) and with pre-amplification prior to detection
(below).

There are two different possible implementations which we study here. In both
implementations, there is a receiver at the end of the last segment. The
implementations differ in the sense that one uses a pre-amplifier directly in
front of the receiver and the other does not. The two corresponding capacities
for the three receiver types ($X\in\\{S,J,D\\}$) can then be computed as

$\displaystyle C_{X,np}(M,P,\tau)$
$\displaystyle=C_{X}(M,P,\tau,1-\tau^{L-1})$ (9)

when no pre-amplifier is used in front of the receiver and as

$\displaystyle C_{X,p}(M,P,\tau)$
$\displaystyle=C_{X}(M,P,1,\tau^{-1}-\tau^{L-1})$ (10)

if pre-amplification is used. These formulas arise from the calculation of
overall loss for the setting without pre-amplification as

$\displaystyle\tau_{L}=\tau\prod_{i=1}^{K-1}G\tau=\tau$ (11)

where $G=1/\tau$ was used, and as $\tau_{L}=1$ if pre-amplification is used.
Moreover, as the noise photons injected by the amplification process are
subject to loss as well, the overall noise (for $K\geq 2$) is calculated using
the geometric series as

$\displaystyle\nu_{K}$ $\displaystyle=\sum_{i=1}^{K-1}\tau^{K-i}(G-1)$ (12)
$\displaystyle=(1-\tau)\frac{1-\tau^{K-1}}{1-\tau}$
$\displaystyle=1-\tau^{K-1}.$

When pre-amplification is used, the last amplifier injects an additional
number of $\tau^{-1}-1$ noise photons. In this case, the noise level satisfies
$\nu\gg 1$. As one can see we have, in both cases, an amount of noise which is
strictly larger than $0$. We discuss our approach in Subsection III-A.

##  III Results and Discussion

###### Theorem 1

The scaling of the entanglement-assisted capacity for $M\to\infty$ is given by
the formula

$\displaystyle C_{E}(M,\tau,P,\nu)=\mathcal{O}(T\cdot P\cdot\log M)$ (13)

where $T=\frac{1+\nu-\tau}{1+\nu}$.

The logarithmic scaling of the capacity with the signal energy has already
been observed early on in [16]. The added value of this work is thus in the
concrete phrasing in terms of conventionally used system parameters allowing
for comparison with the other two technologies, and in the interpretation of
formula 13 in Subsection III-A.

Using formula (13) we may take any constant $c>0$ modelling - for example the
capacity limit of a given amplified link, specified by values $P$ and $\nu$ \-
and ask what values of $M$ are needed such that, for example,
$C_{E}(M,\tau,P,\nu)\geq 10\cdot c$. The approximate answer, which is
increasingly exact as $P/M\to 0$, is then given by

$\displaystyle\log M=10\cdot c\cdot\frac{1+\nu}{P(1+\nu-\tau)}.$ (14)

To show the applicability of this formula, we compute its solution for two
particular cases of a $160km$ link and a $10.000km$ link, both with
$P=10^{16}$.

$\displaystyle{10^{1}}$$\displaystyle{10^{2}}$$\displaystyle{10^{3}}$$\displaystyle{10^{4}}$$\displaystyle{10^{5}}$$\displaystyle{10^{6}}$N$\displaystyle{10^{15}}$$\displaystyle{10^{16}}$$\displaystyle{10^{17}}$$\displaystyle{10^{18}}$Communication
capacityPreamplified final receiver$\displaystyle C_{S,p}$$\displaystyle
C_{J,p}$$\displaystyle C_{E,p}$ Figure 2: Scaling of the capacity in the
setting when pre-amplification is used at the receiver. The bandwidth $B$ is
set to $10^{13}$ and only the number $N$ of spatial modes is varied. Starting
from a value of $N=10^{5}$, the system operates under the condition
$PM^{-1}<1$.
$\displaystyle{10^{1}}$$\displaystyle{10^{2}}$$\displaystyle{10^{3}}$$\displaystyle{10^{4}}$$\displaystyle{10^{5}}$$\displaystyle{10^{6}}$N$\displaystyle{10^{15}}$$\displaystyle{10^{16}}$$\displaystyle{10^{17}}$$\displaystyle{10^{18}}$Communication
capacityNon-preamplified final receiver$\displaystyle C_{S,np}$$\displaystyle
C_{J,np}$$\displaystyle C_{E,np}$ Figure 3: Scaling of the capacity in the
setting with no pre-amplifier. The bandwidth $B$ is set to $10^{13}$ and only
the number $N$ of spatial modes is varied. Starting from a value of
$N=10^{5}$, the system operates under the condition $PM^{-1}<1$.

###  III-A Discussion

There are four major themes which impact the technological plausibility of our
analysis. First, this is the amplification model. Second, this is the assumed
growth of the value $B$. Third, this is the assumed power limit $P$ imposed on
the entire optical link. Fourth, the growth of $N$ is highly relevant. We will
discuss these aspects separately.

####  III-A1 Amplifier Model

The quantum model of amplification as such predicts that, in order to amplify
with a gain $G$, an amount of $G-1$ noise photons has to be tolerated.
Inspection of the condition $PM^{-1}\ll 1$ which defines the regions where
entanglement-assisted communication shows its advantage, indicates that, at
the first amplifier stage, the received signal has a photon number $\tau
PM^{-1}\ll 1$ so that amplification with $G=\tau^{-1}$ adds a number
$\tau^{-1}-1$ of noise photons. If $\tau\ll 1$ as is easily the case in a
situation where $\tau=e^{-0.05\cdot L}$ and $L>1km$, the number of added noise
photons will be much higher than the number of signal photons. This leads to
an incredibly bad signal to noise ratio, thereby imposing a limit on the
growth of capacity in Eq. (2). Since the amount of noise that gets induced
_per channel_ is in this model _independent_ from the energy per channel, and
therefore does not change when $M\to\infty$, the quantity Eq. (4) gets limited
as well. We note that this behaviour changes when the amplifier gain is
adjusted to match the power per pulse $P/M$ [3]. The question as to what
degree the amplifier operating mode or even the amplifier model itself deserve
modifications in a setting with a high number of modes both in spatial and in
frequency domain will need to be addressed separately.

From the perspective of communication system design, the major question is
whether a supposed quantum advantage is kept when moving to continuous
amplification (Ramann amplification) as opposed to the discrete version (using
e.g. Erbium-doped amplifiers) - as described by us. From our discussion of the
noise model Eq. (12) the limit of Ramann amplification leads to a situation
where exactly one photon per slot is added as noise. As one can see from
Theorem 1 and the discussion below, this does not change the general argument
that both conventional detection methods and joint detection methods deliver
asymptotically constant capacity under a power limit under the proposed
amplification rule, while entanglement-assisted communication yields unbounded
(logarithmic) growth.

####  III-A2 Growth of $M$ through $B$

In the work [17], a growth of $B$ over time has been observed which seems to
be well approximated by the formula $B=10^{10}(1+j-2000)$, where $j$ is the
year and the current experiments operating around a value of $200\cdot 10^{9}$
pulses per second are taken as the value of today. The growth of $B$ is
hindered by many practical challenges, among which is the conversion of light
between different spectral domains. Fortunately the recent literature [18]
demonstrates that new ways of accessing more optical spectrum in a cost-
effective way have been found. Taking only the $C$,$L$ and $O$ band in the
existing SMF-28 fiber therefore already justifies that a value of $B=10^{13}$
might become accessible. In addition, the development of the hollow-core fiber
[19] brings about new fiber types with large transmission windows which are
essentially homogeneous in terms of the physical properties affecting
communications, thereby simplifying system design. However, even the very high
value of $B=10^{13}$ is not sufficient to motivate the condition $PM^{-1}\ll
1$. Therefore, in this work we also consider a growth of $M$ induced by other
factors, such as spatial multiplexing.

###  III-B Power Limit $P$

Optical fiber systems are subject to a power limit induced by the nonlinear
response of the fiber to the energy which is emitted into the fiber per second
[20]. In practice, this results in systems utilizing around $P=100mW$ spread
over large numbers of individual channels (the ITU lists $144$ channels for
communication in the C-band). While the Shannon formula, Eq. (2), predicts a
logarithmic capacity growth with the transmission power, the nonlinear
response of the fiber leads to a _decrease_ of the system capacity when the
emitted power rises above a threshold.

This fact leads to the prediction of a so-called _capacity crunch_ where fiber
networks are not able to support the (exponentially) growing demand for data
transmission anymore. As a solution to this problem, one can simply put more
fiber, use more clever ways of handling the data to prevent exponential growth
in the demand or change the fiber design. Among the activities in the latter
direction is the above-mentioned development of hollow-core fibers [19].
Another direction of research pursues multi-mode fiber in general, where one
of the research directions pursues approaches to increase the core diameter of
the fiber so that it supports multiple spatial modes. The different approaches
to multi-mode fiber design are well described in [9].

Due to the novelty of the newly proposed fiber designs, analytical treatments
from the perspective of (quantum) communication systems design are scarce. The
analytical prediction that entanglement-assisted communication may at some
point outperform conventional and even joint detection receiver technologies
however rests on the presence of a power limit.

While the question of an existence of a power limit for fiber transmission is
not on the table, the question of how that power limit depends on the number
$N$ of spatial modes is largely open [9]. Although it might decrease the value
of this work, we would like to note the following: some results [21] seem to
indicate that nonlinearities (which are the major reason for the presence of a
global power limit such as 15 in optical transmission systems) might couple
the power limit to the number $N$ in a way where an increase from $N$ to
$x\dot{N}$ would not inevitably necessitate a reduction of power per channel
by a factor of $x$ anymore.

####  III-B1 Growth of $M$ through $N$

As pointed out in [9], the number of spatial modes $N$ in a fiber of
sufficient core diameter can already today easily reach values of $N=1000$.
This new degree of freedom may separate future network design from the current
one, where standard single-mode fiber is the technology standard for long-haul
transmission systems.

##  IV Derivation of Approximating Formula for $C_{E}$

Our goal here is to derive an approximation to $C_{E}$ which allows us to
study, with high accuracy, its asymptotic scale when the number $M$ of
orthogonal channels, which are subject to a global power constraint

$\displaystyle\textstyle{\sum_{m=1}^{M}P_{m}\leq P},$ (15)

increases. The scaling of the overall system capacity with $M$ is, in this
case, described by the quantities $C_{X,x}(M,P,\tau)$. For $X=C$ and $X=J$ the
limit of these quantities has been stated in [3], where the formula
$\lim_{M\to\infty}C_{C,x}(M,P,\tau)=P\cdot\tau\cdot\log(1+\nu)\nu^{-1}$ can be
found with $\tau$ being the transmittivity of the last segment (which is equal
to $1$ if pre-amplification is used) and $\nu$ being equal to $1-\tau^{K-1}$
for $x=np$ and $\nu_{K}=\tau^{-1}-\tau^{K-1}$ for $x=p$. This value sets the
true benchmark for entanglement-assisted communication to provide an advantage
in terms of data transmission capacity. The question is thus for which values
of $(M,P,\tau)$ the entanglement-assisted capacity can outperform the non-
assisted ones by a given factor.

As one can see from (6), it is helpful to understand the individual scaling of
the involved four terms. In this analysis, it is helpful to study the terms
$(g(\tau\cdot P\cdot M^{-1})-g(d_{+}(\tau,P\cdot M^{-1},\nu)))M$ and
$(g(\tau\cdot P\cdot M^{-1}+\nu)-g(d_{-}(\tau,P\cdot M^{-1},\nu)))M$
separately.

Our approximation rests on the following observations: For every
$\tau\in[0,1]$ and $\nu\geq 0$ it holds $d(\tau,0,\nu)=\nu+1$, which then
implies

$\displaystyle d_{-}(\tau,0,\nu)=0.$ (16)

Obviously it holds $P\cdot M^{-1}\to 0$ as $M\to\infty$. Our second
observation is that the derivative of $d_{-}(1,n,\nu)$ at $n=0$ satisfies

$\displaystyle\partial_{n}d_{-}(\tau,n,\nu)|_{n=0}=(1+\nu-\tau)\cdot(1+\nu)^{-1}.$
(17)

Since $\lim_{\epsilon\to 0}g(\nu+\epsilon)=g(\nu)$,
$\lim_{M\to\infty}d_{+}(\tau,P\cdot M^{-1},\nu)=\nu$ and
$\partial_{n}d_{+}(\tau,P\cdot M^{-1},\nu)|_{n=0}=\tau\nu/(1+\nu)$ we conclude
by L’Hospital’s rule that the term $(g(\tau\cdot P\cdot
M^{-1}+\nu)-g(d_{+}(\tau,P\cdot M^{-1},\nu)))M$ asymptotically approaches a
constant given by $(1-\tau\nu/(1+\nu))\log(1+\nu)/\nu$. Since we are aiming to
study the asymptotic growth of $C_{E}$ we neglect this constant. Thus, the
asymptotic growth of $C_{E}$ must be well approximated by

$\displaystyle C_{E}(M,\tau,P,\nu)\approx M\left(g(\tau
PM^{-1})-g(d_{-}(\tau,PM^{-1},\nu))\right).$ (18)

For small values of $\epsilon>0$ it holds
$g(\epsilon)\approx(1+\epsilon)\epsilon-\epsilon\log(\epsilon)$. Therefore,
setting $T:=(1+\nu-\tau)/(1+\nu)$ and abbreviating $PM^{-1}$ as $\epsilon$ we
arrive at

$\displaystyle C_{E}$ $\displaystyle(M,\tau,P,\nu)\approx
P\left(g(\epsilon)-g(T\epsilon)\right)/\epsilon$ (19)
$\displaystyle\approx\tfrac{P}{\epsilon}\left((1+\epsilon)\epsilon-(1+T\epsilon)T\epsilon-\epsilon\log\epsilon+T\epsilon\log
T\epsilon\right)$ $\displaystyle\approx
P\left(-\epsilon\log\epsilon+T\epsilon\log T\epsilon\right)/\epsilon$
$\displaystyle=P\left(-\epsilon\log\epsilon+\epsilon\log
T\epsilon+(1-T)\epsilon\log T\epsilon\right)/\epsilon$
$\displaystyle=P\left(\epsilon\log T+(1-T)\epsilon\log
T\epsilon\right)/\epsilon$ $\displaystyle\approx P(1-T)\log T\epsilon$
$\displaystyle\approx P(1-T)\log M.$

Here we subsequently neglected all terms converging to constants.

This formula applies to both the pre-amplified and the not pre-amplified case,
where the difference between the two cases is handled by the term $1-T$. It
can be set in comparison to existing ones, and it shows that entanglement-
assisted communication can enable infinite growth even for amplified optical
networks. However, when numbers for existing networks are inserted, this
formula also shows that the parameter regions where entanglement-assisted
communication might be helpful have not yet been reached.

$\displaystyle{10^{4}}$$\displaystyle{10^{6}}$$\displaystyle{10^{8}}$$\displaystyle{10^{10}}$$\displaystyle{10^{12}}$$\displaystyle{10^{14}}$$\displaystyle{10^{16}}$$\displaystyle{10^{18}}$N$\displaystyle{10^{0}}$$\displaystyle{10^{1}}$Communication
capacityApproximation Ratio$\displaystyle C^{A}_{EA}/C^{NA}_{EA}$ Figure 4:
Quality of approximation formula. As the power per pulse per mode $P/(BN)$
approaches $0$, the quality of approximation improves, with the approximation
function $C_{EA}^{A}$ converging the formula $C_{EA}^{NA}$ which we are trying
to approximate.

##  V Conclusions

We produced a formula to calculate the threshold after which a desired factor
(for example a factor of $10$) of an improvement in capacity can be expected
from the use of entanglement-assisted links, where the formula involves total
transmission power $P$ over all channels, number of spatial channels $M$ and
the bandwidth $B$ consumed by the system.

## References

  * [1] Saikat Guha, Quntao Zhuang and Boulat A Bash  “Infinite-fold enhancement in communications capacity using pre-shared entanglement”  In _2020 IEEE International Symposium on Information Theory (ISIT)_ , 2020, pp. 1835–1839  IEEE
  * [2] Vittorio Giovannetti, Seth Lloyd, Lorenzo Maccone and Peter W Shor  “Broadband channel capacities”  In _Physical Review A_ 68.6 APS, 2003, pp. 062323
  * [3] Janis Nötzel and Matteo Rosati  “Operating Fiber Networks in the Quantum Limit”, 2022  arXiv:2201.12397 [quant-ph]
  * [4] Marcin Jarzyna, Raul Garcia-Patron and Konrad Banaszek  “Ultimate capacity limit of a multi-span link with phase-insensitive amplification”  In _45th European Conference on Optical Communication (ECOC 2019)_ , 2019, pp. 1–4  DOI: 10.1049/cp.2019.0742
  * [5] Zuhra Amiri, Boulat A. Bash and Janis Nötzel  “Performance of Quantum Preprocessing under Phase Noise”  arXiv, 2022  DOI: 10.48550/ARXIV.2210.14008
  * [6] Ludwig Kunz, Matteo G.. Paris and Konrad Banaszek  “Noisy propagation of coherent states in a lossy Kerr medium”  In _J. Opt. Soc. Am. B_ 35.2 OSA, 2018, pp. 214–222  DOI: 10.1364/JOSAB.35.000214
  * [7] Saikat Guha  “Structured Optical Receivers to Attain Superadditive Capacity and the Holevo Limit”  In _Phys. Rev. Lett._ 106.24, 2011, pp. 240502  DOI: 10.1103/PhysRevLett.106.240502
  * [8] Simon Sekavčnik and Janis Nötzel  “Effects of Quantum Communication in Large-Scale Networks at Minimum Latency”  In _arXiv preprint arXiv:2210.13267_ , 2022
  * [9] René-Jean Essiambre, Robert W. Tkach and Roland Ryf  “Chapter 1 - Fiber Nonlinearity and Capacity: Single-Mode and Multimode Fibers”  In _Optical Fiber Telecommunications (Sixth Edition)_ , Optics and Photonics  Boston: Academic Press, 2013, pp. 1–43  DOI: https://doi.org/10.1016/B978-0-12-396960-6.00001-8
  * [10] Konrad Banaszek, Ludwig Kunz, Michał Jachura and Marcin Jarzyna  “Quantum limits in optical communications”  In _Journal of Lightwave Technology_ 38.10 IEEE, 2020, pp. 2741–2754
  * [11] Alexander S. Holevo  “Bounds for the Quantity of Information Transmitted by a Quantum Communication Channel”  In _Probl. Inform. Transm._ 9.3, 1973, pp. 177–183  URL: http://mi.mathnet.ru/ppi903
  * [12] A S Holevo  “Coding theorems for quantum communication channels”  In _IEEE Int. Symp. Inf. Theory - Proc._ , 1998, pp. 84  DOI: 10.1109/ISIT.1998.708669
  * [13] Konrad Banaszek and Michał Jachura  “Receiver for receiving information transmitted using very weak light pulses and a method for transmitting information by means of very weak light pulses”, 2018
  * [14] Charles H Bennett, Peter W Shor, John A Smolin and Ashish V Thapliyal  “Entanglement-assisted classical capacity of noisy quantum channels”  In _Physical Review Letters_ 83.15 APS, 1999, pp. 3081
  * [15] Janis Nötzel and Stephen DiAdamo  “Entanglement-Assisted Data Transmission as an Enabling Technology: A Link-Layer Perspective”  In _2020 IEEE International Symposium on Information Theory (ISIT)_ , 2020, pp. 1955–1960  DOI: 10.1109/ISIT44484.2020.9174366
  * [16] A.. Holevo and R.. Werner  “Evaluating capacities of Bosonic Gaussian channels”  arXiv, 1999  DOI: 10.48550/ARXIV.QUANT-PH/9912067
  * [17] Si-Ao Li et al.  “Enabling Technology in High-Baud-Rate Coherent Optical Communication Systems”  In _IEEE Access_ 8, 2020, pp. 111318–111329  DOI: 10.1109/ACCESS.2020.3003331
  * [18] G. Ronniger et al.  “Efficient Ultra-Broadband C-to-O Band Converter Based on Multi-Mode Silicon-on-Insulator Waveguides”  In _2021 European Conference on Optical Communication (ECOC)_ , 2021, pp. 1–4  DOI: 10.1109/ECOC52684.2021.9606033
  * [19] Benoît Debord et al.  “Hollow-Core Fiber Technology: The Rising of “Gas Photonics””  In _Fibers_ 7.2, 2019  DOI: 10.3390/fib7020016
  * [20] René-Jean Essiambre et al.  “Capacity Limits of Optical Fiber Networks”  In _Journal of Lightwave Technology_ 28.4, 2010, pp. 662–701  DOI: 10.1109/JLT.2009.2039464
  * [21] Fatih Yaman et al.  “Long distance transmission in few-mode fibers”  In _Opt. Express_ 18.12 Optica Publishing Group, 2010, pp. 13250–13257  DOI: 10.1364/OE.18.013250


11institutetext: Department of Physics, University of Patras, Patras, Rio,
26504, Greece

# Revisiting thermoelectric effects in the crust of neutron stars

Dionysios Gakis<EMAIL_ADDRESS>Konstantinos N. Gourgouliatos
<EMAIL_ADDRESS>
(Received ?? ?? 2024; accepted ?? ?? 2024)

###### Abstract

Context. Large thermal variations have been observed in neutron stars that
typically are not aligned with density gradients. Such terms may activate the
Biermann battery effect, leading to thermo-electric interactions and the
generation of electromotive force.

Aims. We aim to identify the impact a temperature anisotropy on a neutron
star’s crust can have in the evolution of its magnetic field, through the
thermo-electric terms.

Methods. We consider a neutron star crust with large temperature gradients,
associated with long-lived hot spots, described by a Gaussian-type function
localized. We simulate the interplay between the battery term and the Hall and
Ohmic evolution numerically, for axisymmetric systems.

Results. The results indicate that for crust temperatures of $\sim$$10^{9}$ K
the toroidal field can be amplified up to $\sim$$10^{14}-10^{15}$ G near the
points of maximum temperature gradients, and it locally changes the
architecture of the poloidal field lines. For internal crustal temperatures
around $\sim$$10^{8}$ K, the temperature gradient generates fields of about
two orders of magnitude lower. In such cases, saturation is achieved after
some hundred thousand years, after which the battery and Ohmic dissipation
balance each other.

Conclusions. We conclude that the thermoelectric effect can impact the overall
magnetic field evolution, provided that the thermal gradient is maintained for
a sufficiently long time. Neutron stars endowned with moderate strength
magnetic fields may be affected by the thermoelectric effect if the hotspots
survive for timescales of a few kiloyears.

###### Key Words.:

neutron star – thermoelectric battery mechanism – magnetic field

##  1 Introduction

Neutron stars are known to sustain the strongest magnetic fields found in the
universe, ranging from $10^{8}$ G in the case of millisecond pulsars (e.g.
Possenti et al. 2003) up to $10^{15}$ G for magnetars (e.g. Nakagawa et al.
2009). Still, the origin of these fields remains an open question. Spruit
(2008) reviews specifically the following two scenarios. The field could be
inherited by the progenitor main sequence star through flux conservation (e.g.
Makarenko et al. 2021) or generated by dynamo processes before the supernova
explosion. The bottom line is that they appear as plausible explanations only
for pulsars but not for magnetars; other interesting theories are suggested
for the latter.

Whatever the source is, we have indications that magnetic fields in neutron
stars undergo evolution with time. The majority of neutron stars have strong
magnetic fields on their surfaces ($\sim$$10^{10}$-$10^{13}$ G), while the
ones with stronger magnetic fields, host high-energy phenomena, like flares
and bursts, which are linked to changes of their magnetic field structure
(e.g. Mazets et al. 1999; Turolla et al. 2015; Kaspi & Beloborodov 2017; Coti
Zelati et al. 2018). The observed field strengths spanning many orders of
magnitude is yet another argument towards the supposed evolution of
magnetization in neutron stars.

Broadly speaking, the structure of a typical neutron star may be divided in
its crust and core, which in turn are further split into inner and outer
parts. The outer crust is physically described by a degenerate relativistic
electron gas flowing on a solid crystal Coulomb lattice of protons and ions,
whereas matter consisting of neutron-rich nuclei, at densities between nuclear
and the neutron drip point, constitute the inner crust. Moving deeper, the
outer core is a superconductive superfluid of protons and neutrons with only a
few free electrons. Regarding the innermost part of a neutron star, the inner
core, the extremely high dominating densities restrict us to only speculations
about exotic forms of matter, such as free quarks or other unknown particles
(e.g. Baym et al. 2018). The part of a neutron star between the crust and the
magnetosphere (governed by plasma) is the atmosphere, containing the ocean
(sometimes called envelope). Although it might be possible that the magnetic
field permeates both crust and core of a neutron star, simulations typically
consider solely the crustal magnetic field, and only a handful of recent
studies include the core as well (e.g. Ciolfi & Rezzolla 2013; Moraga et al.
2023; Skiathas & Gourgouliatos 2024). In our paper, we will also assume that
the magnetic flux in the core is frozen, and deal with its evolution in the
crust.

The main processes governing the magnetic field evolution in neutron stars are
typically three (Goldreich & Reisenegger 1992). The Hall drift is the
transport of magnetic flux by a moving electron fluid rising electric currents
and is conservative. The conversion of magnetic energy into heat because of
the crust’s finite conductivity is called Ohmic dissipation. Thirdly, the
ambipolar diffusion, another dissipative process, is caused by the relative
movement of magnetic field and charged particles to neutrons (e.g. Shalybkov &
Urpin 1995). The latter effect is important only in the core and inner crust
of neutron stars though, where neutrons are abundant, and hence is disregarded
in the context of this study. The inner battle between Hall and Ohmic terms in
neutron star crusts has been thoroughly analyzed in literature (e.g.
Hollerbach & Rüdiger 2002; Cumming et al. 2004; Pons & Geppert 2007; Kojima &
Kisaka 2012; Gourgouliatos & Cumming 2014; Gourgouliatos et al. 2016).
However, the magnetic and thermal interactions cannot be considered isolated
in neutron star physics, and due to their non-linear interactions, simulating
their coupling has only been possible recently, in 2-D (e.g. Viganò et al.
2013) and 3-D (e.g. Igoshev et al. 2021; De Grandis et al. 2020; Dehman et al.
2023; Ascenzi et al. 2024). Motivated by such efforts, in this paper, we
examine a yet another possible source of magnetic field production, driven the
Biermann battery mechanism (Biermann 1950).

The possible generation of magnetic field through thermoelectric instabilities
was initially studied in laboratory environments about laser-heated plasmas
(e.g. Stamper et al. 1971; Tidman & Shanny 1974). The Biermann battery
mechanism has been proposed as a mechanism for the formation of magnetic
structures by weak seed fields in astrophysics contexts as well, for instance
in cosmological scales (e.g. Gnedin et al. 2000). Later works motivated by
Biermann (1950) (e.g. Roxburgh 1966) suggested that differential rotation in
fluid stars could produce misalignment in isobars and equipotentials and
consequently, battery effect. Its effect is generally rather small compared to
the dynamo action by helical motion in fluid which is far more effective.
Nonetheless, thermoelectric effects might be important in neutron stars since
dynamo action cannot occur in solid crusts.

The thermoelectric mechanism was first studied in cores of white dwarfs by
Dolginov & Urpin (1980a, b). However, as noted by Blandford et al. (1983), the
mechanism should not be as relevant in white dwarfs as their crystallization
occurs when they are very old and little thermal energy remains. Thus, neutron
stars remain as main candidates for this phenomenon.

In neutron star crusts, the battery mechanism could work the following way.
For some reason, it is possible that thermal anisotropies are formed within
the star. Such temperature gradients, along with density gradients which
appear naturally within depth, may circulate the electron fluid, carrying both
thermal and electric currents. This electric field generated
thermoelectrically acts as a battery, generating a magnetic field in a
vertical direction from the density and temperature gradients.

The first attempt of studying such effects in neutron star crusts was by
Blandford et al. (1983) who extended the work of Urpin & Yakovlev (1980) about
thermogalvanomagnetic effects in degenerate stars111We note that in a few
cases in literature the induction of electromotive force from thermal
gradients is referred to as ”thermomagnetic” or ”thermogalvanomagnetic”. In
this work we shall refer to that as thermoelectric effect, as the term
”thermomagnetic” is related to the change of the magnetization of material due
to temperature variations and such use may lead to confusion. In particular,
they examined two mechanisms of generation of magnetic flux, in the solid
crust and in the melted envelope layer above it. Regarding the former case,
they showed that heat carried by degenerate electrons could cause
thermoelectric instability in the solid crust, resulting in horizontal
magnetic fields that grow exponentially with time, under appropriate
conditions. In the latter case, they found that since pressure and density
gradients must be parallel, the battery mechanism cannot work the same way as
previously. However, the liquid layers contain horizontal magnetic fields as
well, heat flux can drive circulation amplifying the field strength (possibly
leading to dynamo effect), and hence, flux may be convected from the liquid
layers to the solid ones. If these instabilities develop, the field rises
rapidly to $\sim$$10^{12}$ G, with additional growth saturating at
$\sim$$10^{14}$ G, the point at which non-linear perturbations describing the
field evolution must be included, or when the magnetic and the lattice yield
stresses are comparable.

Some years later, Urpin et al. (1986) considered unstable modes in shallower
depths, consisting of free degenerate electrons and fully ionized ions in
liquid or gaseous state, and unlike Blandford et al. (1983) they did not
include hydrostatic equlibrium in their considerations. They worked in the
linear frame as well and showed that thermoelectric modes are unstable if the
star rotates. They were able to deduce toroidal field growth in quite hot
neutron stars in a very narrow region close to surface. The growth rate was
found to depend of the surface temperature and the typical size of the field
fluctuation (with a maximum at around $\sim$100 m). Therefore, the generated
field has a strong non-dipolar character, as large scale fields decay in
linear approximation.

Next, the amplification of neutron star magnetic fields through thermoelectric
effects was further analyzed extensively in a series of six papers (Geppert &
Wiebicke 1991; Wiebicke & Geppert 1991, 1992, 1995; Geppert & Wiebicke 1995;
Wiebicke & Geppert 1996). In the first paper (Geppert & Wiebicke 1991) the
general formalism was presented, which considered an isolated, non-accreting
slowly rotating (thus neglecting deviations from spherical symmetry) neutron
star and examining the field growth in the limit that the Lorentz force can be
safely ignored and magnetic pressure remains negligible compared to other
forms of pressure (i.e. the magnetic field can alter the thermal but not the
density profile of the star). Thermoelectric amplification is sought in the
liquid layer, but as opposed to Urpin et al. (1986), the outer boundary
condition was slightly different, as were the mathematical methods for
simultaneously solving the system of heat transport and induction equations.

In the linear approximation Wiebicke & Geppert (1991), valid for moderate
field strengths ($\lesssim 10^{11}$ G), growth and decay rates of the magnetic
field components for different multipolar components were calculated for
several surface temperatures and gravitation intensities. They found that
small-scale ($n$ $\sim$ 1000) toroidal components grow exponentially in some
years under certain conditions within a thin layer ($\sim$20-100 m) below
surface creating sinusoidal variations with meridional wavelengths of
$\sim$100 m (like the results obtained by Urpin et al. (1986)). Large-scale
modes like dipole and quadruple slowly decay in linear approximation. These
outputs were re-infroced by the analytical approximations performed in
Wiebicke & Geppert (1996). In this last paper, simple formulae for the maximum
position and width of the temperature and field (Gaussian-assumed) functions
were derived. Interestingly, these two functions were found quite symmetric
despite the strong dependence of the parameters on density. It was shown that
the layer and maximum position of the functions loosely depend on $n$, though
with larger $n$, thickness decreases and the maximum position is shifted
outwards. Lastly, they were able to reproduce the described behavior of the
the growth rate in their analytical formula, and reached a straightforward
criterion for the occurrence of thermoelectric instability.

In the other papers, non-linear calculations were done. Wiebicke & Geppert
(1992) investigated the dominating non-linear terms on the toroidal small-
scale magnetic field rapid growth, and yielded that the rate decreasing with
increasing field strength because of some form of self-coupling.
Asymptotically, a stationary solution at $10^{11}$ G was obtained. Saturation
occurs when the characteristic time (given by the inverse of electron Larmor
frequency) becomes comparable with the relaxation time for electron-ion
collisions (in principle, when the electrons are forced to spiral around the
field lines and hence, collisions become improbable). Since the seed field is
a priori unknown, Wiebicke & Geppert (1995) introduced averaged small scale
modes to get the structure of an individual axisymmetric mode, the relative
weights of which are defined by initial conditions (assuming chaotic
distributions). Regardless of the seed field geometry, particularly the
axisymmetric toroidal large-scale modes with even multipolarity are strongly
induced, while all other modes (including poloidal ones) are amplified at the
same rate and their (absolute) field values are significantly lower.
Additionally, conversely with the linear approximation, large-scale modes can
be amplified rapidly in non-linear coupling with small-scale modes, and some
selection rules of this coupling were identified. Wiebicke & Geppert (1996)
described the evolution of the strongly induced axisymmetric large-scale
toroidal field modes with even multipolarities (without considering their
coupling with poloidal field). Their evolution seems rather insensitive to the
underlying star properties, and may be split into a phase of exponential
growth followed by inward diffusion along with further enhancement gradually
slowing down.

In spite of those advancements and the somewhat promising results,
thermoelectric effects on neutron stars have been quite overlooked in recent
years. Some encouragements to explore this direction have indeed been
formulated (for instance we refer an interested reader to Geppert (2017) for a
review of magneto-thermal evolution of neutron stars) along with individual
speculations about their application (e.g. Konenkov & Geppert 2001; De Grandis
et al. 2020; Geppert et al. 2021). Nevertheless, a thorough numerical
investigation about the battery term and its battle with the Hall and Ohmic
effects has not been done. This is the goal of the current study.

This paper is organized as follows. Section 2 contains the set of the
underlying equations describing the problem and in Section 3, we discuss the
numerical approach implemented for solving them. We continue with analyzing
the results from the simulations (Section 4) and discuss their implications
for realistic neutron stars in Section 5. A summary of this work is given in
Section 6.

##  2 Mathematical setup

###  2.1 Magnetic induction equation

We work in the one-fluid approximation of the crust, at which ions are
considered static forming a crystal Coulomb lattice, and electrons have
freedom to move, thus, they are the only carriers of electric current and
thermal energy. The typical electron velocities are low enough to resort to
non-relativistic electron magneto-hydrodynamics (e-MHD). Thus the electric
current density $\bm{j}$ is given by:

$\bm{j}=-en_{e}\bm{\upsilon},$ (1)

where $e$ is the elementary electron charge $n_{e}$ the electron number
density and $\bm{\upsilon}$ the mean velocity of electrons

In the case of finite electric conductivity $\sigma$, the electric field
$\bm{E}$ may be written through the generalized Ohm’s law:

$\bm{E}=-\frac{\bm{\upsilon}\times\bm{B}}{c}+\frac{\bm{j}}{\sigma}-\frac{1}{e}S_{e}\nabla
T,$ (2)

in which $c$ is speed of light, $\bm{B}$ the magnetic field, $T$ the
temperature and $S_{e}$ represents the electron entropy, defined as:

$S_{e}=\left(\frac{\pi^{4}}{3n_{e}}\right)^{1/3}\frac{k_{B}^{2}T}{c\hbar},$
(3)

where $k_{B}$ is Boltzmann’s constant and $\hbar$ is the reduced Planck
constant. In this paper, the last term of equation 2 is of particular
importance, which expresses an additional component of electric field, the one
formed by a temperature gradient. In its general form, this electric field
component is given by:

$\bm{E}=\hat{G}\nabla T,$ (4)

where the tensor $\hat{G}$ is the thermopower and has two parts: an isotropic
one (Seebeck term) and an anisotropic one (Ettingshausen-Nernst term). For
simplicity, we restrict ourselves to the first part only. It has the form
$G_{ij}=-S_{e}/e\delta_{ij}$ for a completely degenerate Fermi gas
(Gourgouliatos et al. 2022), as written above.

The magnetic fields have currents as sources, according to Ampère’s law:

$\bm{j}=-\frac{c}{4\pi}\nabla\times\bm{B}.$ (5)

Faraday’s law states that any change of magnetic field induces an electric
field as:

$\frac{\partial\bm{B}}{\partial t}=-c\nabla\times\bm{E}.$ (6)

Combining the above equations, we obtain the full magnetic induction equation:

$\displaystyle\frac{\partial\bm{B}}{\partial
t}=-c\nabla\times\left(\frac{1}{4\pi
en_{e}}(\nabla\times\bm{B})\times\bm{B}+\frac{c}{4\pi\sigma}\nabla\times\bm{B}-\frac{1}{e}S_{e}\nabla
T\right)\,$ (7)

In this equation, the first term expresses the Hall effect, which is dependent
on $\bm{B}^{2}$, and the second refers to Ohmic dissipation (Jones 1988;
Goldreich & Reisenegger 1992). The third term is the thermal battery, which is
independent of $\bm{B}$, and its strength is modulated by the $T$ gradient and
the thermopower.

We now closely follow the formalism by Gourgouliatos & Cumming (2014). In
axisymmetry, we decompose the magnetic field into a poloidal ($\bm{B}_{P}$)
and a toroidal ($\bm{B}_{T}$) component, and trace their evolution in terms of
two scalar quantities, $\Psi(r,\theta)$ and $I(r,\theta)$:

$\bm{B}=\bm{B}_{P}+\bm{B}_{T}=\nabla\Psi\times\nabla\phi+I\nabla\phi\,,$ (8)

where we use spherical coordinates ($r$, $\theta$, $\phi$) and
$\nabla\phi=\hat{\phi}/rsin\theta$. By construction, the above expression
fulfills Gauss’ law for the magnetic field, $\nabla\cdot\bm{B}=0$, provided
that $\Psi$ and $I$ are differentiable quantities. The function $\Psi$ is
referred as flux function (identical in form to Stokes stream function), since
$2\pi\Psi$ is the poloidal magnetic flux passing through a spherical cap of
radius $r$ and opening angle $\theta$. $I$ is the poloidal current function,
as $cI/2$ expresses the poloidal current passing through the same spherical
cap. The magnetic field lines lie on the surfaces where $\Psi$ is constant,
while the current or electron flow lines lie on surfaces of constant $I$.

We now express equation 7 in terms of the quantities $\Psi$ and $I$. As in
Reisenegger et al. (2007), we define the quantity

$\chi={\frac{c}{4\pi{\rm e}n_{\rm e}r^{2}\sin^{2}\theta}}\,,$ (9)

the toroidal current as

$\bm{j}_{T}={c\over 4\pi}\nabla\times\bm{B}_{P}=-{c\over
4\pi}\Delta^{*}\Psi\nabla\phi\,,$ (10)

and the angular velocity of the electrons is

$\Omega=-{j_{T}\over n_{\rm e}{\rm e}r\sin\theta}=\chi\Delta^{*}\Psi\,,$ (11)

where Grad-Shafranov operator is given by

$\Delta^{*}=\frac{\partial^{2}}{\partial
r^{2}}+\frac{\sin\theta}{r^{2}}\frac{\partial}{\partial\theta}\left(\frac{1}{\sin\theta}\frac{\partial}{\partial\theta}\right).$
(12)

After some vector calculus manipulation, the coupled differential equations
for $\Psi$ and $I$ are as follows. The one related to poloidal field evolution
is

$\displaystyle\frac{\partial\Psi}{\partial t}-r^{2}\sin^{2}\theta\chi(\nabla
I\times\nabla\phi)\cdot\nabla\Psi=\frac{c^{2}}{4\pi\sigma}\Delta^{*}\Psi\,,$
(13)

and the one concerning toroidal field evolution is

$\displaystyle\frac{\partial I}{\partial
t}+r^{2}\sin^{2}\theta(\left(\nabla\Omega\times\nabla\phi\right)\cdot\nabla\Psi+I\left(\nabla\chi\times\nabla\phi\right)\cdot\nabla
I)$
$\displaystyle=\frac{c^{2}}{4\pi\sigma}\left(\Delta^{*}I-\frac{1}{\sigma}\nabla
I\times\nabla\sigma\right)+\frac{c}{e}r\sin\theta(\nabla S_{e}\times\nabla
T)\cdot\hat{\phi}\,.$ (14)

hence, the toroidal and poloidal components of the field couple strongly to
each other in their evolution.

###  2.2 Crust properties

We now outline our model of the microphysical parameters within the crust,
used to solve equations 13 and 2.1. Firstly, the crust is modeled as a
spherical shell, covering the 5% outermost part of star in terms of radius.
Hence, the crust extends from the crust-core boundary, $r_{cc}$, up to the
neutron drip point, $r_{nd}$, and $r_{cc}/r_{nd}=0.95$. As the radius of our
modeled star is 11.3 km, we are only considering the last $\sim$500 m of it.

We assume that the density throughout the crust varies as a function of the
depth with the following formula deduced by Chamel & Haensel (2008) (adopted
by Karageorgopoulos et al. (2019) as well):

${\rho}(r)=\left(1+\left(\frac{r_{nd}-r}{r_{nd}-r_{cc}}\right)\frac{\rho_{cc}}{\rho_{out}}\right)\rho_{out},$
(15)

where $\rho_{cc}=1.3\times 10^{14}$ g cm-3 is the density at the base of the
crust and $\rho_{out}$ is the density at the neutron drip point.

Crystallization occurs when the electrostatic potential energy between nuclei
dominates over thermal energy. Their ratio $\Gamma$ is the so-called Coulomb
parameter:

$\Gamma=\frac{Z^{2}e^{2}}{a_{I}k_{B}T}.$ (16)

The critical value of $\Gamma$ is 175, which corresponds to the melting point,
i.e. for larger values crystallization happens (Potekhin & Chabrier 2000).
Assuming a pure iron crust ($Z=26$, $A=56$), the ion sphere radius (also
referred in literature as interionic spacing) is:

$a_{I}=\left(\frac{4\pi n_{I}}{3}\right)^{-1/3}.$ (17)

The density $\rho_{out}$ at the bottom of the ocean may be found using the
numerical density $n_{I}$ (for a fully ionized crust):

$\rho_{out}=n_{I}mA,$ (18)

which by construction is a function of temperature. The liquid part of a
neutron star (ocean) exists for lower densities, but is not included in our
considerations due to its comparably small size ($\sim$ some cm) and the fact
that e-MHD fails in this region. We also note that we ignore any small
fluctuations that may occur on the surface of the star and assume that
$\rho_{out}$, throughout the star.

In general, the electron number density $n_{e}$ is given by:

$n_{e}=\frac{Z}{A}(1-x_{fn})\frac{\rho}{m_{b}},$ (19)

in which $x_{fn}$ is the fraction of free neutron outside ions and $m_{b}$ the
atomic mass unit. Lander & Gourgouliatos (2019) performed a polynomial fitting
of $Z$, $A$, $x_{fn}$ in terms of $\rho$ utilizing the equation of state by
Douchin & Haensel (2001) and their result may be approximated accurately
enough by the following straightforward expression:

$n_{e}=10^{36}(1.5\tilde{\rho}^{2/3}+1.9\tilde{\rho}^{2})\ cm^{-3},$ (20)

which was also used later by Gourgouliatos & Lander (2021). The scaled density
is $\tilde{\rho}=\rho/\rho_{cc}$.

Regarding conductivity, we assume that it scales as $\sigma\propto
n_{e}^{2/3}$, based on phonon scattering or impurity scattering (Cumming et
al. 2004), being $\sigma(r=r_{cc})=10^{24}s^{-1}$ at the crust-core boundary.
It is not considered a function of temperature in order to keep our code
manageable (as it would require additional solving of the heat transfer
equations) and focus on the physics of the magnetic field evolution.

The typical timescale for the Ohmic evolution is:

$t_{\rm Ohm}\sim{4\pi\sigma L^{2}\over c^{2}}=566\ {\rm Myr}\ \left({L\over
11.3\ {\rm km}}\right)^{2}\left({\sigma\over 10^{24}\ {\rm s^{-1}}}\right),$
(21)

whereas the Hall timescale is

$t_{\rm Hall}\sim{4\pi{\rm e}L^{2}n_{\rm e}\over cB}={81.2\ {\rm Myr}\over
B_{14}}\ \left({L\over 11.3\ {\rm km}}\right)^{2}\left({n_{\rm e}\over
10^{36}\ {\rm cm^{-3}}}\right),$ (22)

and the timescale of the battery is

$t_{\rm battery}\sim{eL^{2}B\over cTS_{e}}=33.6\ {\rm Myr}{B_{14}\over
T_{9}^{2}}\ \left({L\over 11.3\ {\rm km}}\right)^{2}\left({n_{\rm e}\over
10^{36}\ {\rm cm^{-3}}}\right),$ (23)

where $B_{14}$ is the strength of the magnetic field expressed in units of
$10^{14}$ G and $T_{9}$ is the temperature in $10^{9}$ K. The relative
importance of the Hall and Ohm terms is gauged by the ratio of these
timescales, which is called magnetic Reynolds number (or Hall parameter):

$R_{M}=\frac{\sigma|B|}{cen_{e}}.$ (24)

$R_{M}$ does vary through crust, but in our computations we account for its
maximum value. In our simulations, we keep it at maximum a few times of 102,
because when it reaches some 103 simulations become unstable.

##  3 Numerical Setup

###  3.1 Numerical implementation

Our simulations are performed in a Fortran 90 environment. We consider the
magnetic field of the star as axisymmetric, thus our calculations are
restricted in 2-D. We discretized the numerical domain in $r$ and
$\mu=cos\theta$ (in order to avoid infinities at the poles). The resolution is
typically $100\times 100$ but we have experimented in several cases with a
double analysis to confirm the numerical convergence of our outputs.

Derivatives are numerically calculated using finite differences; a central
difference scheme is used to calculate spatial derivatives, and a three- and
five-point stencil for second and third derivatives, respectively. The
numerical method adopted for the time intergations of equations 13 and 2.1 is
the second order Adams-Bashford (AB2). The timestep is adaptable, according to
Courant condition (Courant et al. 1952):

$\Delta t\leq\frac{L}{\upsilon_{max}}$ (25)

in which $L$ is axis size and we set as $\upsilon_{max}$ the maximum electron
plasma velocity, as we work in e-MHD.

Three sets of boundary conditions are adopted, enclosing the numerical domain
of our simulation world. Regarding the surface, we assume a vacuum magnetic
field outside of star using a multipole fit. This implies that $B_{\phi}=0$ at
the surface ($r=r_{nd}$), giving a multipolar current-free poloidal magnetic
field, $I(r_{nd},\theta)=0$, so that no current flows from the star to
magnetosphere. The continuity of $B_{\theta}=0$ comes as follows. In the
absence of current outside of the star, we may write $\nabla\times\bm{B}=0$,
hence the vacuum field can be expressed as $\bm{B}=\nabla V$, where $V$ is a
scalar function. Along with $\nabla\cdot\bm{B}=0$, we get Laplace equation,
$\nabla^{2}V=0$. Working in axisymmetry, we decompose the external field in a
sum of multipolar expansion and get (Marchant et al. 2011):

$V(r,\theta)=\sum_{l=1}^{\infty}\frac{a_{l}}{r^{l+1}}P_{l}(\cos\theta),$ (26)

where $P_{l}(\cos\theta)$ are the Legendre polynomials. We have:

$B_{r}(r,\theta)=\frac{\partial V}{\partial
r}=-\sum_{l=1}^{\infty}(l+1)\frac{a_{l}}{r^{l+2}}P_{l}(\cos\theta),$ (27)

and hence the coefficients are calculated as:

$a_{l}=-r^{l+2}\frac{2l+1}{2l+2}\int_{0}^{\pi}P_{l}(\cos\theta)B_{r}(r,\theta)\sin\theta\,d\theta.$
(28)

Thus, $B_{\theta}$ is given by:

$B_{\theta}(r,\theta)=\frac{1}{r}\frac{\partial
V}{\partial\theta}=\sum_{l=1}^{\infty}\frac{a_{l}}{r^{l+2}}\frac{dP_{l}(\cos\theta)}{d\theta}=-\frac{1}{r\sin\theta}\frac{\partial\Psi}{\partial
r},$ (29)

giving the boundary condition:

$\frac{\partial\Psi(r_{nd},\theta)}{\partial
r}=-\sin\theta\sum_{l=1}^{\infty}\frac{a_{l}}{r_{nd}^{l+1}}\frac{dP_{l}(\cos\theta)}{d\theta}.$
(30)

Field lines are not allowed to cross the symmetry axis, that is $\Psi$ and $I$
functions are set to zero on the axis of the star, i.e.
$\Psi(r,0)=\Psi(r,\pi)=I(r,0)=I(r,\pi)=0$. Concerning the inner boundary,
$dI/dr=0$, whereas we consider two cases for the boundary condition for
$\Psi$. For magnetic field confined in the crust, it is
$\Psi(r_{cc},\theta)=0$ or equal to a function $f$,
$\Psi(r_{cc},\theta)=f(\theta)$, when allowed to permeate to the core.
However, under the Meissner effect assumption, which states that any flux is
expelled out of core faster than thermal or magnetic timescales (e.g. Lander
2014; Passamonti et al. 2017), no magnetic field evolution is allowed in the
core itself or at the crust-core boundary, as is done in most relevant
studies.

###  3.2 Initial conditions

We examine two scenarios for the initial poloidal magnetic field. In the
option at which the field is restricted in the crust, the initial poloidal
flux is given by the expression:

$\Psi_{1}(t=0)=\Psi_{0}\left(1.05-\frac{r}{r_{nd}}\right)\left(\frac{r}{r_{nd}}-\frac{r_{cc}}{r_{nd}}\right)sin^{2}\theta,$
(31)

while it is slightly different in the alternative at which the field is
supposed to thread the core (without evolving within it though, as explained
in Section 3.2):

$\Psi_{2}(t=0)=\Psi_{0}\left(1.05-\frac{r}{r_{nd}}\right)\frac{r}{r_{nd}}sin^{2}\theta.$
(32)

In both situations, the intensity of the field is controlled by the parameter
$\Psi_{0}$. The poloidal current function is fixed at zero at the beginning of
the simulations.

The temperature profile of the neutron star is of crucial importance for our
calculations too. We choose to adopt a Gaussian-type function for the
temperature distribution:

$T(r,\theta)=T_{0}+\lambda T_{0}\
exp\left(-\frac{(r-r_{0})^{2}}{\sigma_{r}^{2}}\right)exp\left(-\frac{(\mu-\mu_{0})^{2}}{\sigma_{\mu}^{2}}\right),$
(33)

which spreads over both $r$ and $\mu$. In this formula, a hot region centered
at the point $(r_{0},\mu_{0})$ and of $1\sigma$ spread (drop to a factor of
$1/\mathrm{e}$) controlled by the parameters $\sigma_{r}$ and $\sigma_{\mu}$
in radial and meridional directions, respectively, is assumed to be
superimposed to a background of constant temperature $T_{0}$. The magntitude
of the additional temperature at the point $(r_{0},\mu_{0})$ is $\lambda$
times $T_{0}$, where $\lambda$ can be any positive number determining the
steepness of the temperature gradient.

We are of course aware that in most types of stars, their temperature should
rise moving radially inner towards their core. However, since our model
considers only the outermost 5% of the neutron star (its crust), it is rather
safe to assume that a uniform background temperature distribution governs this
region. The additional heat region was added to account for any temperature
anisotropy that may rise at the outer layers of the neutron star (the reasons
for this behavior will be explained more thoroughly in Section 5), which is
subsequently supposed to spread throughout the star. As for the choice of the
temperature function (Gaussian profile), lacking any essential evidence of its
mathematical formulation, we argue that exponential decays appear rather
naturally on such environments. Besides, in any tests we performed with other
functions with similar profiles (e.g. quadratic) we find no signs of
significant changes in our results.

We are mostly interested in the influence of the temperature gradient to the
magnetic field architecture, hence the point of interest is mainly limited to
a region near the pole, in which the temperature shows its maximum value. For
most of our simulations we set the central point equal to
($r_{0}=0.99,\mu_{0}=0.99$) and the standard deviations $\sigma_{r}=0.1$ and
$\sigma_{\mu}=0.1$, while unless otherwise stated we keep $\lambda=1$ (i.e.
the maximum temperature is twice the background temperature). Sections 4.1 and
4.2 are dedicated to exploring the field evolution with different combinations
of $T_{0}$ (high and low, accordingly) and $\Psi_{0}$. A representative
temperature profile, for $T_{0}=10^{9}$ K, is presented in Fig. 1 (the star’s
crust has been scaled up by a factor of eight for visualisation purposes in
all such figures in this paper). On the other hand, alternate values for the
rest of the parameters in equation 33 ($\lambda$, $r_{0}$, $\mu_{0}$,
$\sigma_{r}$, $\sigma_{\mu}$) are examined in Section 4.3.

Lastly, it should be remarked that throughout the simulation time-span, the
temperature is kept constant in every run. In other words, we do not account
for the evolution of temperature with time and assume that some form of heat
source maintains it for some reason at the same levels. While this is
apparently an over-optimistic assumption (especially for the timescales we are
considering), we think that is a necessary approximation we need to make at
the current time, in order to examine battery effects in neutron stars and
keep our calculations manageable. Future studies should include the thermal
and magnetic interactions in a more realistic context.

Figure 1: Temperature profile for $T_{0}=10^{9}$ K, $\lambda=1$, $r_{0}=0.99$,
$\theta_{0}=0.99$, $\sigma_{r}=0.1$, $\sigma_{\mu}=0.1$.

##  4 Results

In our calculations we do not account for the rotation of the star, hence in
all simulations the poloidal lines remain in the same meridional plane. Below
we present the outputs from different models of the temperature profile. All
cases presented here are simulated for a maximum time of 1 million years, in
order to be able to observe the saturation of the system each time.

###  4.1 Model A. High temperature

As stated in Section 3.2 our basic layout has the following parameter values:
$r_{0}=0.99,\mu_{0}=0.99$, $\sigma_{r}=0.1$, $\sigma_{\mu}=0.1$ and
$\lambda=1$. In this subsection, we examine a rather high main crust
temperature, that of $T_{0}=10^{9}$ K.

We begin with initial poloidal strength of $\Psi_{0}=10.0$ and the case of
$\Psi_{1}$ (eq. 31). Three snapshots of the evolution of the magnetic field,
at 10, 100 and 1000 kyr are given in Fig. 2. The colorbar on the left refers
to the value of the toroidal magnetic field, whereas the density of the
contour lines gives an estimation of the poloidal magnetic field.

Soon after the start of the simulation, a region of strong toroidal field is
formed near the point $(r_{0},\theta_{0})$, which gradually overpasses in
strength the structures formed by the Hall effect (yellow and light blue
regions covering the two hemispheres in Fig. 2a). As time passes, the region
of strong field due to the battery shifts deeper and towards the axis of the
star, due to velocity and density differences between different parts of the
crust. The addition of the thermoelectric current, and its resulting toroidal
field creation, have as an effect that poloidal lines have their initial shape
changed, and start to twist around the battery region. These lines eventually
tend to encircle the region of strong toroidal field, creating interesting
arcades and multipoles near the pole, but are kept almost intact at distant
regions.

a |  b |  c
---|---|---
Figure 2: The evolution of Model A through time in the case of $\Psi_{1}$ and
$\Psi_{0}=10.0$. In each snapshot, the same number of poloidal lines are drawn
so that direct comparisons may be made. Figure 3: The evolution of toroidal
and poloidal energies through time of Model A in the case of $\Psi_{1}$,
$\Psi_{0}=10.0$.

We now turn on our attention on the evolution of the poloidal and toroidal
energies through time, which are computed as follows:

$E_{pol}=\frac{1}{8\pi}\int_{V}(B_{r}^{2}+B_{\theta}^{2})\,dV,$ (34)
$E_{tor}=\frac{1}{8\pi}\int_{V}B_{\phi}^{2}\,dV,$ (35)

where $B_{r}=-r^{-2}\partial\Psi/\partial\mu$,
$B_{\theta}=-1/(rsin\theta)\partial\Psi/\partial r$ and
$B_{\phi}=I/(rsin\theta)$ are the components of the magnetic field in
spherical coordinates. Their resulting values are plotted in Fig. 3 as a
function of time. We notice that although the poloidal energy is monotonically
descreasing until reaching zero, the toroidal energy starts from zero (since
no initial toroidal field is assumed) and increases up the value of about
$2.5\times 10^{43}$ erg, at which it stabilizes.

The transition to constant values for the toroidal and poloidal energies is
located around 400 kyr (consistent with the analytical calculation of eq. 23,
with the appropriate scale), at which we assume there is the saturation point.
In other words, saturation occurs when the energy and field provided by the
thermoelectric battery are compensated by the Ohmic losses. After saturation
is achieved, the structure and values of the field components remain constant,
with only a few minor oscillations of the poloidal lines observed.

We have explored the evolution of the above configuration for other values of
$\Psi_{0}$ as well, ranging from 0.01 up to 100.0.The final snapshots (at
simulated time of 1 Myr) are presented in Fig. 4 for the cases of
$\Psi_{0}=0.1,1.0,10.0$ and $100.0$. Simulations with even larger $\Psi_{0}$
were not successful. Interestingly, the larger the initial poloidal field, the
sparser the final poloidal lines appear, but typically the final architecture
of the lines is rather similar, especially for the low-$\Psi_{0}$ cases.
Additionally, all poloidal energies gradually drop to zero, despite starting
at different values depending on $\Psi_{0}$.

Regarding the toroidal field, it appears that the battery term dominates
immediately over the Hall term for lower $\Psi_{0}$, whereas it takes more
time to achieve that for larger values. Still, we note that saturation is
observed at roughly the same time for all of them, while the final value of
the maximum toroidal field in the star is also the same, $3.16\times 10^{14}$
G, as is evident in Fig. 4. This is a value representative for magnetars,
while the rest of the star is described by values around $5\times 10^{12}$ G,
two orders of magnitude lower. Apart from that, for all choices of $\Psi_{0}$
the toroidal energy also reaches the same steady value, but for
$\Psi_{0}\gtrsim 10$, the toroidal energy reaches a peak before dropping at
this value.

We have also experimented with even higher temperatures. For instance, setting
$T_{0}=5\times 10^{9}$ K resulted in a maximum of $\bm{B}_{T}=2.35\times
10^{15}$, noticing even sparser structure of poloidal lines. However, since
even higher temperatures are not physically acceptable for the crusts of the
even highly magnetized magnetars, we do not consider this scenario further.

a |  b
---|---
c |  d
Figure 4: The final snapshots for Model A in the case of $\Psi_{1}$ for different values of $\Psi_{0}$ (0.1 - a, 1.0 - b, 10.0 - c, 100.0 - d). In each snapshot, the same number of poloidal lines are drawn so that direct comparisons may be made. a |  b |  c
---|---|---
Figure 5: The evolution of Model A through time in the case of $\Psi_{2}$ and
$\Psi_{0}=0.01$. In each snapshot, the same number of poloidal lines are drawn
so that direct comparisons may be made. Figure 6: The evolution of toroidal
and poloidal energies through time of Model A in the case of $\Psi_{2}$,
$\Psi_{0}=0.01$.

Things are different for the case of $\Psi_{2}$ (eq. 32). Simulations are only
successful for very low initial values $\Psi_{0}$, as in the opposite case
features of extremely large field are created, rendering the code unstable.
Hence, we simulate temperatures of $T_{0}=10^{9}$ K and $\Psi_{0}=0.01$ (given
in Fig. 5), or even lower. The respective evolution of the energies is given
in Fig. 6. The high-$\bm{B}_{T}$ region is observed to shift to inner depths
and smaller polar angles. At roughly the same time with the previous
configurations, saturation is again achieved (though at a slightly larger
toroidal energy value, while the poloidal energy remains practically zero
during the whole timespan), apparently stabilizing the geometry of poloidal
lines after $\sim$ 700 kyr. The maximum value of toroidal field is now a bit
larger ($3.27\times 10^{14}$ G) and the minimum ones around $2.50\times
10^{13}$ G.

###  4.2 Model B. Low temperature

Keeping the same set of parameters as in Section 4.1
($r_{0}=0.99,\mu_{0}=0.99$, $\sigma_{r}=0.1$, $\sigma_{\mu}=0.1$,
$\lambda=1$), we now consider lower crust temperatures, namely $T_{0}$ being
between $10^{7}$ and $10^{8}$ K. Starting with the case of $\Psi_{1}$, Fig. 7a
presents the final snapshot (1 Myr) when $T_{0}=10^{8}$ K and $\Psi_{0}=1.0$.
Having examined initial poloidal fields $\Psi_{0}$ in the range 0.01 to 10.0
(for lower temperatures the battery term in the resulting toroidal field can
prevail only over weak poloidal fields, the final outputs are insensitive to
this choice, all reaching a steady state of $2.92\times 10^{12}$ G at maximum
for the toroidal field. At this low value, the toroidal field is not strong
enough to twist the existing poloidal lines, hence no such feature is observed
now, apart from minor oscillations after about 700 kyr. As a result, the
toroidal energy rises to values of $\sim$$10^{39}$ erg, while the poloidal
energy follows a similar behavior as in Model A, diminishing to zero as time
passes. The rest of the star (far from the battery region) is maintained at
fields $\sim$$10^{10}$ G.

As expected from eq. 7 lower temperatures lead to even lower thermally induced
magnetic fields. We simulated temperature $T_{0}=5\times 10^{7}$ K and the
final result in the case of $\Psi_{1}$ and $\Psi_{0}=1.0$ is shown Fig. 7b.
$\bm{B}_{T}$ has the value of $7.2\times 10^{11}$ G near the pole, which is
almost an order of magnitude weaker than the field in the case of Fig. 7a, as
is the toroidal energy after saturation.

We notice that for the temperatures examined here ($10^{7}$ \- $10^{8}$ K) the
toroidal fields achieved are in the range of $\sim$$10^{11}$ – $10^{12}$ G,
representative for typical pulsars. At even lower temperatures, the mechanism
becomes rather ineffective. Particularly, for $T_{0}\lesssim 10^{7}$ K no
battery feature is apparent and the resulting magnetic field evolution is
mainly described by the combination of Hall and Ohm terms.

As far as the field $\Psi_{2}$ scenario, even at $10^{8}$ K the Hall term
outweighs battery, as no signs of the latter are observed. We considered
$\Psi_{0}$ as low as 0.01, since as explained in Section 4.1, the code cannot
sustain larger values. Hence, the case of radial field assumed to penetrate
the core is no further considered.

a |  b
---|---
Figure 7: The final snapshots for Model B in the case of $\Psi_{1}$ and
$\Psi_{0}=1.0$ for $T_{0}=10^{8}$ K (a) and $T_{0}=5\cdot 10^{7}$ K (b).

###  4.3 Experiments with the temperature function

This subsection is dedicated in modifying the parameters we were keeping
constant ($r_{0}$, $\mu_{0}$, $\sigma_{r}$, $\sigma_{\mu}$, $\lambda$) and
exploring their effect in the results. Unless otherwise stated, we will adopt
as our main simulation setup the one presented in Figs. 2, 3 (Model A,
$\Psi_{1}$, $\Psi_{0}=10.0$, $T_{0}=10^{9}$ K) and will make any changes (one
at a time) regarding the temperature function parameters to this one.

We first experiment with the position of the point having the maximum
temperature throughout the crust (referred as heat point). We assumed so far
that the maximization of temperature occurs in $r_{0}=0.99$. The reason behind
this option lies in the fact that magnetospheric currents (which contribute to
the augmentation of temperature near the polar cap) are attenuated within the
first $\sim$100 m of the star (at depth corresponding to approximately 1% of
the star’s radius). While this seems a realistic choice, there may be
accumulation of heat at other points inner to that. Consequently, we ran two
tests with $r_{0}=0.95$ and $r_{0}=0.90$. In the former case, results appear
almost identical compared to when $r_{0}=0.99$ in terms of final values and
the architecture of the toroidal and poloidal field, their energies, and
saturation times. Setting $r_{0}=0.90$ (the heat point is located at the base
of the crust) results in some changes though. The final snapshot, at 1 Myr, is
given in Fig. 8, in which the same number of poloidal lines as in Fig. 2 have
been drawn for comparison. It seems that the maximum $\bm{B}_{T}$ values have
dropped to almost half ($\sim$$1.68\times 10^{14}$ G), as does the toroidal
energy reached (about $1.0\times 10^{43}$ erg). By examining the time
evolution, we additionally conclude that saturation occurs at approximately
100 kyr later.

Figure 8: The final snapshot for Model A in the case of $\Psi_{1}$ and
$\Psi_{0}=10.0$ for $r_{0}=0.90$.

Concerning the position of the heat point in the meridional, we simulated the
symmetric (with respect to the equator) configuration, that is
$\mu_{0}=-0.99$. Fig. 9 shows the final output (again the same number of
poloidal lines as in Fig. 2 have been drawn). Clearly, the heat point located
at the south pole yields the opposite values (i.e. the battery region is
surrounded by strong negative magnetic fields), but the absolute values remain
the same. Moreover, the same happens with the other results (saturation time
and magnetic energies). Despite that, moving the heat point at other latitudes
(e.g. at the equator or at mid-latitudes) was not possible. The development of
strong toroidal field ($10^{14}$ G) in regions far from the poles twists
excessively the existing poloidal lines in narrow areas preventing the code
from working properly.

Figure 9: The final snapshot for Model A in the case of $\Psi_{1}$ and
$\Psi_{0}=10.0$ for $\mu_{0}=-0.99$.

Next, we test the spread of Gaussian function in $r$ and $\mu$. Simulating an
idealized neutron star, we have the freedom to vary the values $\sigma_{r}$
and $\sigma_{\mu}$, as the actual temperature profile of a neutron star is
unknown. We choose to examine 3 additional scenarios. The first one is the
most conservative, at which $\sigma_{r}$ and $\sigma_{\mu}$ are reduced 10
times, at 0.01. We then double them, at 0.2, while an even more optimistic
scenario, at which these parameters are as large as 0.5 each, is examined as
well.

These 3 cases, along with the default one, are presented in Fig. 10a-d at
simulated time of 1 Myr (all of the figures have the same number of poloidal
lines). In principle the small yellow region of strong toroidal field is
shifted, as expected, since different spreads also change the point of maximum
temperature gradient. Low standard deviations restrict the battery region near
the pole only, though with multipoles evident farther from it. The maximum
$\bm{B}_{T}$ values look similar for $\sigma_{r}$, $\sigma_{\mu}$ equal to
0.01 and 0.2. Observing the latter case’s evolution reveals that around 400
kyr the poloidal lines start to increase in density forming wavy structures,
whereas the toroidal field reaches a maximum value ($1.83\times 10^{14}$ G)
having more than one regions of strong field, which are later unified as it
drops and saturates at values $1.58\times 10^{14}$ G at roughly 600 kyr. The
large spread of the last case, Fig. 10d, along with the fact that the point of
maximum temperature gradient now belongs in the core, have as a consequence
that the toroidal field is even lower, $7.0\times 10^{13}$ G, and more time
(about 900 kyr) is required for saturation. Lastly, we notice that for all of
these configurations, because of twisting in narrow region near the initial
heat point, field lines become almost radial and parallel around it.

a |  b
---|---
c |  d
Figure 10: The final snapshots for Model A in the case of $\Psi_{1}$ and
$\Psi_{0}=10.0$ for 4 different values of $\sigma_{r}$ and $\sigma_{\mu}$:
(0.01 - a, 0.1 - b, 0.2 - c, 0.5 - d).

Moving on, here is a short discussion of the parameter $\lambda$. We consider
the following situation: the typical temperature is set at $T_{0}=10^{8}$ K
(Model B), which is a more realistic temperature for a neutron star, but we
assume that the temperature rose 11 times at the heat point, that is
$\lambda=10$. The final figure, with the rest model choices being $\Psi_{1}$
and $\Psi_{0}=10.0$, is given in Fig. 11. It is noticed that the field values
reached are almost two orders of magnitude larger ($1.26\times 10^{14}$ G)
than in the case of $\lambda=1$ (Fig. 7). Hence, the battery is now capable of
twisting the poloidal field. We argue that the scenario in this simulation is
quite realistic, as it does not require a high temperature along the crust,
except from only a region that may be greatly enhanced.

Figure 11: The final snapshot for Model B in the case of $\Psi_{1}$,
$\Psi_{0}=10.0$ and $T_{0}=10^{8}$ K for $\lambda=10$.

While it may be tempting to experiment more with variations of these
parameters, or combinations of them, and how they may influence the
architecture and evolution of the thermoelectrically-induced magnetic field,
we restrict ourselves to the experiments performed above. This is because the
temperature distribution of neutron star has not conclusively been determined,
and hence, any deep investigation of the chosen model parameters might be too
theoretical, not resembling any particular relevance to similar future
studies. The practicability and robustness of our model and results are
discussed in the next Section.

##  5 Discussion

It is possible that temperature anisotropies are formed near the surface of a
neutron star. Some possible energy sources that create a hotter region at a
specific part of the crust are the following. Firstly, it can be a result of
the thermo-magnetic evolution of the star. For instance, the magnetic field
itself might form hotspots through phenomena like reconnection. Secondly, it
can be generated by backflowing magnetospheric currents. The current returning
from the light cylinder forms a closed loop and it may heat the region near
the rim of the polar cap as it enters through the outer layers of crust.
Thirdly, another possible mechanism only relevant to non-isolated neutron
stars is accretion. Any accretion (e.g. through a disk) of different material
(e.g. hydrogen) other than the iron-nickel onto the crust could lead to
thermonuclear reactions, which increases the local temperature as well.

Whatever the mechanism is, it is believed that temperature differences may
appear rather naturally on neutron stars. Numerous sources of strong magnetic
fields and ages in the range of $10$ kyr have hotspots whose temperatures
exceed $5\times 10^{6}$ K (Viganò et al. 2013). For instance, in the
magnificent seven, thermal pulses from kilometer-sized areas are observed
(Haberl 2007; Turolla 2009; Potekhin et al. 2015; De Grandis et al. 2022).
These pulses correspond to temperatures of $10^{6}$ K while the characteristic
ages of the sources are in the order of Myrs. Moreover, J0030+0451 has been
thermally mapped in X-rays (Riley et al. 2019; Miller et al. 2019), and both
circular and linear thermal regions of $\sim$$10^{6}$ K were identified near
its south pole. While these temperatures are much smaller than the ones
considered in this work, this source is rather old. Thus, surface temperature
anisotropies can be maintained for rather long timescales.

Thermal blanketing (Gudmundsson et al. 1983; Beznogov et al. 2021) is another
important point. The temperature at the interior of the crust can be much
higher than the one on the surface, as also assumed by thermal evolution
calculations (Viganò et al. 2013; De Grandis et al. 2020). There, the surface
temperature is scaled to the crustal one by application of the following
formula (Igoshev et al. 2021):

$\displaystyle\frac{T_{c}}{10^{8}\ \mathrm{K}}=\left(\frac{T_{s}}{10^{6}\
\mathrm{K}}\right)^{2}\,,$ (36)

where $T_{c}$ is the temperature at the top of the crust and $T_{s}$ is the
observed temperature at the surface. Thus, an observed temperature $5\times
10^{6}$ K, which is the case in several sources, corresponds to a crustal
temperature of $2.5\times 10^{9}$ K, which lies well within the range of our
calculations.

In our model, we assumed that temperature anisotropies appear in an idealized
neutron star crust, ignoring their origin. Nevertheless, the main limitation
of our considerations lies in the timescale at which this temperature gradient
is maintained. We implicitly assumed that the thermal energy of the star
somehow rises and remains constant over a period of hundreds of thousand
years. This implies that the temperatures we consider are somehow conserved by
an external mechanism, although no such term exists in our equations
(including the heat conduction in our set of equations is beyond the scopes of
this paper). Realistically, temperature anisotropies drop, e.g. by neutrino
emission (Yakovlev et al. 2001), over the course of time (e.g. outburts at any
point of the crust typically last for 5-10 years or accretion may last for
some thousand years). While we acknowledge this limitation, we think our study
will be crucial first step in integrating battery terms in equations regarding
magnetic field evolution in neutron stars. Having calculated the maximum
possible impact of a battery term, we show that at least for an idealized
environment thermoelectric effects are rather important and can potentially
provide solutions to open problems regarding neutron star physics, hence
thermo-magnetic interactions should not be omitted.

The suitability of the temperature function used was thoroughly reviewed in
Section 3.2. We now discuss the repercussions of the choices of its relevant
parameters. We have assumed that the hot region of the star is merely
connected to the polar cap. Its position in the meridional is given as a first
approximation (for dipolar field):

$\theta_{0}=sin^{-1}\left(\frac{r_{ns}\Omega}{c}\right),$ (37)

implying that it is related to the rotation of star ($r_{ns}$ stands for the
neutron star’s radius and $\Omega$ is the rotational angular velocity). As we
do not model a specific type of neutron star and restrict ourselves to an
idealized situation, we did not underwent into much detail into the exact
position of the central point, and placing it near the pole(s) will suffice.
Other than that, the temperature gradient may be formed by outbursts as well,
as explained, giving us the freedom to experiment with the geometry of the
temperature profile. In that sense, that is why we examined a wide range for
the extent of the temperature anisotropy.

Lastly, other minor limitations of our considerations are related to the
models adopted for density and conductivity. Particularly, conductivity is a
function of density, which in turn is dependent on the ocean temperature;
thus, conductivity is only indirectly modulated by temperature. Later work may
adopt more accurate microphysics, including both electrical and thermal
conductivities regulated by the respective temperature at each point.

##  6 Conclusions

Under appropriate conditions, magnetic field may be generated and enhanced
through battery mechanisms, one of which is by temperature anisotropies
causing thermal flux. High temperatures gradients are required to achieve
that, and neutron star crusts are a potential candidate for such a phenomenon.
In this paper, we considered an idealized neutron star and showed that
magnetic field enhancement might be possible by thermoelectric effects.

In our model we made the unrealistic approximation that a temperature profile
is kept constant throughout the simulated time, in order to assess the maximum
effect the battery term has in the magnetic field evolution. Regions of
localized strong toroidal field are formed (ranging from $10^{11}$ – $10^{15}$
G, depending on the temperature assumed), around of which poloidal field
twists creating multipoles. For high temperatures ($\gtrsim 10^{8}$ K), the
thermoelectric effect dominates over Hall drift, and the system requires some
thousand years so that Ohmic losses are in equilibrium with the battery, with
only minor oscillations observed afterwards.

Later studies should shed more light into the actual role of thermoelectric
effects in neutron stars. More realistic models need to be adopted, that is
considering temperature decay and spread of heat through time or taking into
account the thermo-magnetic interactions. Follow-up experiments with
microphysical parameters and the temperature profile can be performed in the
future as well, using the most up-to-date data deduced by observations.
Lastly, we expect that 3-D simulations instead of axisymmetry might be worth
investigating.

###### Acknowledgements.

We thank Dimitris Ntotsikas for his assistance in the optimization of the code
used for this work. We thank Davide De Grandis for his insightful comments.
This work was supported by computational time granted from the National
Infrastructures for Research and Technology S.A. (GRNET S.A.) in the National
HPC facility - ARIS - under project ID pr015026/simnstar”.

## References

  * Ascenzi et al. (2024) Ascenzi, S., Viganò, D., Dehman, C., et al. 2024, arXiv preprint arXiv:2401.15711
  * Baym et al. (2018) Baym, G., Hatsuda, T., Kojo, T., et al. 2018, Reports on Progress in Physics, 81, 056902
  * Beznogov et al. (2021) Beznogov, M., Potekhin, A., & Yakovlev, D. 2021, Physics Reports, 919, 1
  * Biermann (1950) Biermann, L. 1950, Zeitschrift Naturforschung Teil A, 5, 65
  * Blandford et al. (1983) Blandford, R., Applegate, J., & Hernquist, L. 1983, Monthly Notices of the Royal Astronomical Society, 204, 1025
  * Chamel & Haensel (2008) Chamel, N. & Haensel, P. 2008, Living Reviews in relativity, 11, 1
  * Ciolfi & Rezzolla (2013) Ciolfi, R. & Rezzolla, L. 2013, MNRAS, 435, L43
  * Coti Zelati et al. (2018) Coti Zelati, F., Rea, N., Pons, J. A., Campana, S., & Esposito, P. 2018, MNRAS, 474, 961
  * Courant et al. (1952) Courant, R., Isaacson, E., & Rees, M. 1952, Communications on pure and applied mathematics, 5, 243
  * Cumming et al. (2004) Cumming, A., Arras, P., & Zweibel, E. 2004, The Astrophysical Journal, 609, 999
  * De Grandis et al. (2022) De Grandis, D., Rigoselli, M., Mereghetti, S., et al. 2022, Monthly Notices of the Royal Astronomical Society, 516, 4932
  * De Grandis et al. (2020) De Grandis, D., Turolla, R., Wood, T. S., et al. 2020, The Astrophysical Journal, 903, 40
  * Dehman et al. (2023) Dehman, C., Viganò, D., Pons, J. A., & Rea, N. 2023, Monthly Notices of the Royal Astronomical Society, 518, 1222
  * Dolginov & Urpin (1980a) Dolginov, A. & Urpin, V. 1980a, Soviet Astronomy, vol. 24, Mar.-Apr. 1980, p. 177-182. Translation Astronomicheskii Zhurnal, vol. 57, Mar.-Apr. 1980, p. 302-310., 24, 177
  * Dolginov & Urpin (1980b) Dolginov, A. & Urpin, V. 1980b, Astrophysics and Space Science, 69, 259
  * Douchin & Haensel (2001) Douchin, F. & Haensel, P. 2001, A&A, 380, 151
  * Geppert (2017) Geppert, U. 2017, Journal of Astrophysics and Astronomy, 38, 1
  * Geppert et al. (2021) Geppert, U., Basu, R., Mitra, D., Melikidze, G., & Szkudlarek, M. 2021, Monthly Notices of the Royal Astronomical Society, 504, 5741
  * Geppert & Wiebicke (1991) Geppert, U. & Wiebicke, H.-J. 1991, Astronomy and Astrophysics Supplement Series (ISSN 0365-0138), vol. 87, no. 2, Feb. 1991, p. 217-228., 87, 217
  * Geppert & Wiebicke (1995) Geppert, U. & Wiebicke, H.-J. 1995, Astronomy and Astrophysics, v. 300, p. 429, 300, 429
  * Gnedin et al. (2000) Gnedin, N. Y., Ferrara, A., & Zweibel, E. G. 2000, The Astrophysical Journal, 539, 505
  * Goldreich & Reisenegger (1992) Goldreich, P. & Reisenegger, A. 1992, Astrophysical Journal, 395, 250
  * Gourgouliatos & Cumming (2014) Gourgouliatos, K. & Cumming, A. 2014, Monthly Notices of the Royal Astronomical Society, 438, 1618
  * Gourgouliatos et al. (2022) Gourgouliatos, K. N., De Grandis, D., & Igoshev, A. 2022, Symmetry, 14, 130
  * Gourgouliatos & Lander (2021) Gourgouliatos, K. N. & Lander, S. K. 2021, Monthly Notices of the Royal Astronomical Society, 506, 3578
  * Gourgouliatos et al. (2016) Gourgouliatos, K. N., Wood, T. S., & Hollerbach, R. 2016, Proceedings of the National Academy of Sciences, 113, 3944
  * Gudmundsson et al. (1983) Gudmundsson, E. H., Pethick, C., & Epstein, R. I. 1983, Astrophysical Journal, Part 1 (ISSN 0004-637X), vol. 272, Sept. 1, 1983, p. 286-300., 272, 286
  * Haberl (2007) Haberl, F. 2007, Ap&SS, 308, 181
  * Hollerbach & Rüdiger (2002) Hollerbach, R. & Rüdiger, G. 2002, Monthly Notices of the Royal Astronomical Society, 337, 216
  * Igoshev et al. (2021) Igoshev, A. P., Hollerbach, R., Wood, T., & Gourgouliatos, K. N. 2021, Nature Astronomy, 5, 145
  * Jones (1988) Jones, P. 1988, Monthly Notices of the Royal Astronomical Society, 233, 875
  * Karageorgopoulos et al. (2019) Karageorgopoulos, V., Gourgouliatos, K., & Contopoulos, I. 2019, Monthly Notices of the Royal Astronomical Society, 487, 3333
  * Kaspi & Beloborodov (2017) Kaspi, V. M. & Beloborodov, A. M. 2017, ARA&A, 55, 261
  * Kojima & Kisaka (2012) Kojima, Y. & Kisaka, S. 2012, Monthly Notices of the Royal Astronomical Society, 421, 2722
  * Konenkov & Geppert (2001) Konenkov, D. & Geppert, U. 2001, Monthly Notices of the Royal Astronomical Society, 325, 426
  * Lander & Gourgouliatos (2019) Lander, S. & Gourgouliatos, K. 2019, Monthly Notices of the Royal Astronomical Society, 486, 4130
  * Lander (2014) Lander, S. K. 2014, MNRAS, 437, 424
  * Makarenko et al. (2021) Makarenko, E. I., Igoshev, A. P., & Kholtygin, A. 2021, Monthly Notices of the Royal Astronomical Society, 504, 5813
  * Marchant et al. (2011) Marchant, P., Reisenegger, A., & Akgün, T. 2011, Monthly Notices of the Royal Astronomical Society, 415, 2426
  * Mazets et al. (1999) Mazets, E. P., Aptekar, R. L., Butterworth, P. S., et al. 1999, ApJ, 519, L151
  * Miller et al. (2019) Miller, M., Lamb, F. K., Dittmann, A., et al. 2019, The Astrophysical Journal Letters, 887, L24
  * Moraga et al. (2023) Moraga, N. A., Castillo, F., Reisenegger, A., Valdivia, J. A., & Gusakov, M. E. 2023, arXiv preprint arXiv:2309.14182
  * Nakagawa et al. (2009) Nakagawa, Y. E., Mihara, T., Yoshida, A., et al. 2009, Publications of the Astronomical Society of Japan, 61, S387
  * Passamonti et al. (2017) Passamonti, A., Akgün, T., Pons, J. A., & Miralles, J. A. 2017, Monthly Notices of the Royal Astronomical Society, 469, 4979
  * Pons & Geppert (2007) Pons, J. A. & Geppert, U. 2007, Astronomy & Astrophysics, 470, 303
  * Possenti et al. (2003) Possenti, A., D’Amico, N., Manchester, R., et al. 2003, The Astrophysical Journal, 599, 475
  * Potekhin & Chabrier (2000) Potekhin, A. Y. & Chabrier, G. 2000, Physical Review E, 62, 8554
  * Potekhin et al. (2015) Potekhin, A. Y., De Luca, A., & Pons, J. A. 2015, Space Sci. Rev., 191, 171
  * Reisenegger et al. (2007) Reisenegger, A., Benguria, R., Prieto, J. P., Araya, P. A., & Lai, D. 2007, Astronomy & Astrophysics, 472, 233
  * Riley et al. (2019) Riley, T. E., Watts, A. L., Bogdanov, S., et al. 2019, The Astrophysical Journal Letters, 887, L21
  * Roxburgh (1966) Roxburgh, I. 1966, Monthly Notices of the Royal Astronomical Society, 132, 201
  * Shalybkov & Urpin (1995) Shalybkov, D. & Urpin, V. 1995, Monthly Notices of the Royal Astronomical Society, 273, 643
  * Skiathas & Gourgouliatos (2024) Skiathas, D. & Gourgouliatos, K. N. 2024, MNRAS[arXiv:2401.08979]
  * Spruit (2008) Spruit, H. C. 2008, in American Institute of Physics Conference Series, Vol. 983, 40 Years of Pulsars: Millisecond Pulsars, Magnetars and More, ed. C. Bassa, Z. Wang, A. Cumming, & V. M. Kaspi, 391–398
  * Stamper et al. (1971) Stamper, J., Papadopoulos, K., Sudan, R., et al. 1971, Physical Review Letters, 26, 1012
  * Tidman & Shanny (1974) Tidman, D. A. & Shanny, R. A. 1974, The Physics of Fluids, 17, 1207
  * Turolla (2009) Turolla, R. 2009, in Astrophysics and Space Science Library, Vol. 357, Astrophysics and Space Science Library, ed. W. Becker, 141
  * Turolla et al. (2015) Turolla, R., Zane, S., & Watts, A. L. 2015, Reports on Progress in Physics, 78, 116901
  * Urpin et al. (1986) Urpin, V., Levshakov, S., & Yakovlev, D. 1986, Monthly Notices of the Royal Astronomical Society, 219, 703
  * Urpin & Yakovlev (1980) Urpin, V. & Yakovlev, D. 1980, Soviet Astronomy, 24, 425
  * Viganò et al. (2013) Viganò, D., Rea, N., Pons, J. A., et al. 2013, Monthly Notices of the Royal Astronomical Society, 434, 123
  * Wiebicke & Geppert (1991) Wiebicke, H.-J. & Geppert, U. 1991, Astronomy and Astrophysics (ISSN 0004-6361), vol. 245, no. 1, May 1991, p. 331-340., 245, 331
  * Wiebicke & Geppert (1992) Wiebicke, H.-J. & Geppert, U. 1992, Astronomy and Astrophysics (ISSN 0004-6361), vol. 262, no. 1, p. 125-130., 262, 125
  * Wiebicke & Geppert (1995) Wiebicke, H.-J. & Geppert, U. 1995, Astronomy and Astrophysics, Vol. 294, p. 303-312 (1995), 294, 303
  * Wiebicke & Geppert (1996) Wiebicke, H.-J. & Geppert, U. 1996, Astronomy and Astrophysics, v. 309, p. 203-212, 309, 203
  * Yakovlev et al. (2001) Yakovlev, D., Kaminker, A., Gnedin, O. Y., & Haensel, P. 2001, Physics Reports, 354, 1