iint iiint iiiint idotsint

# Particle acceleration at ultrarelativistic, perpendicular shock fronts

John G. Kirk, Brian Reville, Zhi-Qiu Huang
Max-Planck-Institut für Kernphysik, Postfach 10 39 80, 69029 Heidelberg,
Germany

###### Abstract

Using an eigenfunction expansion to solve the transport equation, complemented
by Monte-Carlo simulations, we show that ultrarelativistic shocks can be
effective particle accelerators even when they fail to produce large amplitude
turbulence in the downstream plasma. This finding contradicts the widely held
belief that a uniform downstream magnetic field perpendicular to the shock
normal inhibits acceleration by the first order Fermi process. In the
ultrarelativistic limit, we find a stationary power-law particle spectrum of
index $s=4.17$ for these shocks, close to that predicted for a strictly
parallel shock.

###### keywords:

radiation mechanisms: non-thermal — acceleration of particles — gamma rays:
bursts

##  1 Introduction

Discovered in the late 1970s, the theory of diffusive shock acceleration has
established itself as the primary mechanism discussed in connection with the
acceleration of cosmic rays, and has also found many other applications. The
generalisation of this mechanism to mildly relativistic shock fronts followed
roughly a decade later, and, in the early 2000’s, it was found that particles
repeatedly crossing ultrarelativistic shocks can be accelerated into a power
law spectrum of index $s=4.23$ (for recent reviews, see Bell, 2014; Sironi et
al., 2015).

However, in contrast to the nonrelativistic case, an ultrarelativistic shock
that overruns a region containing a uniform magnetic field is generically
superluminal, in the sense that its speed, when projected onto a magnetic
field line, exceeds that of light. At first sight, this poses a problem, since
acceleration requires repeated crossings of the shock front, and a particle
downstream of the shock cannot catch up with it by simply diffusing along a
magnetic field line. Since cross-field diffusion is generally strongly
suppressed, this suggests that the relativistic extension of the diffusive
shock acceleration mechanism might be ineffective (Begelman & Kirk, 1990).

This problem does not arise if the downstream field is effectively scrambled
by strong fluctuations on the scale of a gyroradius (Achterberg et al., 2001).
However, a mechanism for producing such fluctuations has so far not been
identified. For the case of shock propagation into a weakly magnetized plasma,
where the Weibel instability is thought responsible for the formation of the
shock front, particle-in-cell simulations have shown that acceleration is
facilitated by non-resonant scattering on the Weibel-induced filaments (e.g.
Sironi & Spitkovsky, 2009; Plotnikov et al., 2018; Vanthieghem et al., 2020).
However, analytical considerations suggest that scattering exclusively
mediated by such non-resonant interactions is not sufficient to provide the
required cross-field transport in the downstream plasma above a critical
particle energy (Lemoine & Pelletier, 2010; Reville & Bell, 2014; Huang et
al., 2022). Furthermore, such short length-scale fluctuations are susceptible
to damping in the hot downstream plasma, which further reduces the critical
energy (Chang et al., 2008; Keshet et al., 2009; Sironi et al., 2013; Lemoine,
2015)

In this paper, we re-assess these arguments by solving a simple model of a
relativistic, perpendicular shock, thereby demonstrating quantitatively that
particles are accelerated into a power-law distribution whose index lies very
close to that predicted for the idealized, parallel shock case (Kirk et al.,
2000). In our model, energetic particles are assumed to diffuse in angle,
whilst being deflected by a uniform magnetic field that is perpendicular to
the shock normal. We solve this model for the case in which transport upstream
is scattering dominated, i. e., the particles there are unmagnetized, whereas
particles are magnetized when downstream and follow essentially unscattered
trajectories, gyrating about the regular, uniform magnetic field. Two
complementary techniques are employed. First, we use a generalisation of the
eigenfunction approach (Kirk et al., 2000) that takes into account the full,
two-dimensional anisotropy imposed on the particle distribution, including the
drift induced by the downstream magnetic field. This technique is applied in
the ultrarelativistic limit of large shock Lorentz factor, $\Gamma_{\rm
s}\rightarrow\infty$, which enables the eigenfunctions to be found in closed
form. The second method of solution is a Monte-Carlo simulation of many
individual trajectories, using a code described more fully in a companion
paper (Huang et al., 2023, henceforth “ZH”). This is applied here to shocks
ranging from mildly to highly relativistic and shows that, as $\Gamma_{\rm s}$
rises from $2$ to $50$, the power-law index approaches the asymptotic value of
$s=4.17$ found by the eigenfunction method. This value is close to, but
slightly harder than the result obtained previously for the case where
scattering dominates both up and downstream.

In section 2 we formulate the equations describing particle transport, and
present details of the first method of solution. Results of both methods,
consisting of the power law index of accelerated particles and the angular
dependence of their distribution function at the shock front, are presented
and compared in section 3. In section 4 we discuss the physics underlying our
assumptions and their range of applicability and speculate on the implications
of the results for more realistic cases.

##  2 Method

###  2.1 Transport equation

In the presence of isotropic scattering in angle and gyration about a uniform
field, the transport equation governing the phase-space density $f$ of
relativistic particles is given by eq (1) of Takamoto & Kirk (2015),
(henceforth “TK”). Mixed coordinates are used in this equation, with Cartesian
coordinates in configuration space measured in a frame in which the shock is
at rest in the plane $x=0$. The upstream plasma occupies the region $x>0$, the
downstream $x<0$. Momentum space coordinates, on the other hand, are expressed
in spherical polar coordinates, measured in the local (upstream or downstream)
rest frame of the plasma. In this paper, we depart from the notation of TK and
use the shock normal as the axis for these coordinates. Then, the momentum of
an ultrarelativistic particle, in units of $c\,\times\,$the rest mass, is
$\vec{p}=\left(\gamma,\mu,\phi\right)$, with $\gamma$ the Lorentz factor,
$\arccos\mu$ the polar angle to the shock normal and $\phi$ the azimuthal
angle about this axis. The direction of motion of the shock front in this
reference frame is $\mu=1$, and $\mu=0$, $\phi=0$, is the direction of the
uniform magnetic field $\vec{B}$. Solutions are sought that are stationary as
seen in the shock rest frame, and have no dependence on the spatial
coordinates $y$ and $z$. eq (12) of TK is then

$\frac{2\Gamma^{2}\eta c}{\omega_{\rm g}}\left(\mu-\beta\right)\frac{\partial
f}{\partial
x}=\frac{\partial}{\partial\mu}\left(1-\mu^{2}\right)\frac{\partial
f}{\partial\mu}+\frac{1}{1-\mu^{2}}\frac{\partial^{2}f}{\partial\phi^{2}}\\\
+\frac{2\Gamma\eta\mu\cos\phi}{\sqrt{1-\mu^{2}}}\frac{\partial
f}{\partial\phi}-2\Gamma\eta\sqrt{1-\mu^{2}}\sin\phi\frac{\partial
f}{\partial\mu}\,,$ (1)

where $c\beta$ is the shock velocity measured in the local rest frame of the
plasma, $\Gamma=\left(1-\beta^{2}\right)^{-1/2}$ is the Lorentz factor of the
shock front and we have assumed $\gamma\gg\Gamma$. The quantity $\eta$ in this
equation is the ratio of the gyro frequency of an accelerated particle in the
uniform field, $\omega_{\rm g}=\left|eB/\gamma mc\right|$, to $\Gamma$ times
the scattering frequency. Following Reville & Bell (2014), it can be written
in terms of the magnetization parameter $\sigma_{\rm turb}$ associated with
the strength $\delta B$ of the magnetic fluctuations responsible for
scattering, their characteristic length scale $\lambda$, the magnetization
parameter $\sigma_{\rm reg}$ of the uniform, or regular upstream field and the
ion plasma frequency $\omega_{\rm i}$:

$\displaystyle\eta$ $\displaystyle=\left(\frac{mc}{m_{\rm
p}\left<\gamma\right>\omega_{\rm i}\lambda}\right)\left(\frac{\sigma_{\rm
reg}^{1/2}}{\sigma_{\rm turb}}\right)\left(\frac{\gamma}{\Gamma}\right)\,,$
(2)

where $m$ and $m_{\rm p}$ are the rest masses of the accelerating particles
and that of the species dominating the plasma inertia respectively, and
$\left<\gamma\right>m_{\rm p}c^{2}$ is the mean energy per plasma particle.
When $\eta\ll 1$, scattering dominates the transport process and deflections
by the regular magnetic field are unimportant. We then refer to the particles
as being unmagnetized. On the other hand, when $\eta\gg 1$, the particles are
magnetized and follow essentially unscattered trajectories in the uniform,
regular magnetic field. Equation (1) applies in both the upstream and the
downstream regions, where the speeds of the relativistic shock are denoted by
$\beta_{\rm s}$ and $\beta_{\rm d}$, respectively, and the corresponding
Lorentz factors are $\Gamma_{\rm s}\gg 1$ and $\Gamma_{\rm d}\sim 1$.
Analogous notation, $\eta_{\rm s,d}$, is used for the magnetization parameter.

###  2.2 Transport upstream

As is the case for particles accelerated at a parallel, ultrarelativistic
shock front, we expect the stationary particle distribution in the upstream
medium to be concentrated in a narrow cone of opening angle $\sim
1/\Gamma_{\rm s}$ about the shock normal. Consequently, it is convenient to
replace $\mu$ by the stretched variable

$\displaystyle\xi$ $\displaystyle=\left(1-\mu\right)/\left(1-\beta_{\rm
s}\right)$ (3) $\displaystyle\approx 2\Gamma_{\rm s}^{2}\left(1-\mu\right)$

(Kirk & Schneider, 1989) which is zero for particles moving directly along the
shock normal and unity for those moving parallel to the shock front. The
phase-space density is now to be regarded as a function of $x$ and the
momentum space variables $(\gamma,\xi,\phi)$ that lie in the domain
$\gamma>1$, $0\leq\xi<\infty$, $0\leq\phi<2\pi$ which we denote by $A$. Then,
inserting the definition (3) into (1), expanding in powers of the small
parameter $1/\Gamma_{\rm s}$, and assuming the particles are unmagnetized:

$\displaystyle\eta_{\rm s}\ll 1\,,$ (4)

the transport equation upstream becomes

$\displaystyle\frac{\partial^{2}f}{\partial\xi^{2}}+\frac{1}{\xi}\frac{\partial
f}{\partial\xi}+\frac{1}{4\xi^{2}}\frac{\partial^{2}f}{\partial\phi^{2}}-\left(\frac{1}{\xi}-1\right)\frac{\partial
f}{\partial\hat{x}}$ $\displaystyle=0\,,$ (5)

where we have introduced the dimensionless coordinate
$\hat{x}=x\left(4\Gamma_{\rm s}^{2}\omega_{\rm g}\right)/\left(\eta c\right)$.

Equation (5) is separable, resulting in an exponential dependence of $f$ on
$\hat{x}$ and two eigenvalue problems, one for the $\xi$-dependence and one
for the $\phi$-dependence. Taken together with the boundary conditions (i)
$f\rightarrow 0$ as $x\rightarrow\infty$, (ii) $f$ bounded at $\xi=0$ and
$\infty$, (iii) $f$ periodic in $\phi$ with period $2\pi$ and (iv) $f$
invariant under a change of sign in the component of momentum parallel to the
downstream magnetic field: $f(\gamma,\xi,\phi,x)=f(\gamma,\xi,\pi-\phi,x)$,
each eigenvalue problem is self-adjoint. Therefore, the solution can be
expanded in the two-dimensional eigenfunctions $Q_{i}$:

$\displaystyle f\left(\gamma,\xi,\phi,x\right)$
$\displaystyle=F(\gamma)\sum_{i=1}^{\infty}\textrm{e}^{\Lambda_{i}\hat{x}}a_{i}Q_{i}\left(\xi,\phi\right)\,,$
(6)

where the $a_{i}$ are constants and the eigenfunctions are orthogonal over the
weighting function $\xi-1$:

$\displaystyle\iint_{\vec{p}\in
A}\\!\\!\\!\\!\textrm{d}\xi\textrm{d}\phi\,Q_{i}\left(\xi,\phi\right)\left(\xi-1\right)Q_{j}\left(\xi,\phi\right)$
$\displaystyle=0,\quad i\neq j.$ (7)

In Appendix A we give explicit expressions for the eigenvalues $\Lambda_{i}$
and the eigenfunctions $Q_{i}$.

Figure 1: Sample trajectory for a positively charged particle with $\vec{p}\in
R$ as seen in the downstream frame. The shock velocity and magnetic field are
directed along the positive $x$ and $z$ axes respectively. The magnitude of
the particle’s velocity $\beta_{\perp}$ in the $x-y$ plane is a constant of
motion, and its direction along the positive $x$ axis when the phase $\Phi$ is
an integer multiple of $2\pi$. The filled magenta circles denote example
positions where the particle enters and exits the downstream, at phases
$\Phi_{0}$ and $\Phi_{+}$ respectively. To overtake the shock, the
$x$-component of the particle velocity must exceed that of the shock,
requiring $-\textrm{arccos}\zeta<\Phi_{+}<\textrm{arccos}\zeta$, where
$\zeta=\beta_{\rm d}/\beta_{\bot}$. This range is delimited by the two green
circles. For particles with $\vec{p}\in U$, we make the replacements
$\Phi_{0}\rightarrow\Phi_{-}$ and $\Phi_{+}\rightarrow\Phi_{0}$. Further
details are provided in Appendix B.  Figure 2:  Convergence properties of the
approximation scheme: (a) the power-law index $s$ as a function of $i_{\rm
max}$, the number of eigenfunctions used in eq (12), (b) and (c) the phase
space distribution $\bar{f}$ at the shock front ($\hat{x}=0$) as seen in the
frame in which the shock front is at rest, i. e., at constant particle Lorentz
factor in this frame ($\bar{\gamma}$). $\bar{\mu}$ is the cosine of the angle
between the particle momentum and the shock normal, $\bar{\phi}$ the azimuthal
angle about this axis. Particles entering the upstream have $\bar{\mu}>0$. (b)
shows $\bar{f}$ averaged over phase, (c) $\bar{f}$ averaged over $\bar{\mu}$.
(d) shows a slice of $\bar{f}$ for particles that graze the shock,
$\bar{\mu}=0$. Results are plotted for $i_{\rm max}=1,2$ and the fully
converged $i_{\rm max}=30$, as well as for an intermediate value that
indicates the rapidity of convergence. In all cases the speed of the
downstream plasma speed is $\beta_{\rm d}=1/3$. In (b), (c) and (d) only the
relative values of the functions plotted is physically significant.

###  2.3 Transport downstream

Transport in the downstream plasma is potentially more complex, since the
accelerated particles are not concentrated in a narrow beam when viewed from
the frame in which the plasma is at rest. Here, we follow the arguments
presented by Reville & Bell (2014), who conclude that deflection in the
uniform magnetic field perpendicular to the shock normal dominates the
transport process for particles with sufficiently large Lorentz factor, i. e.,

$\displaystyle\eta_{\rm d}\gg 1\,.$ (8)

We concentrate on these high energy particles, because those of lower energy
can be scattered in the turbulence generated at the shock front, and there is
general agreement that the dominance of scattering on both sides of a
relativistic shock results in a power law spectrum $f\propto\gamma^{-s}$ with
index $s\approx 4.2$.

Figure 3: The full angular distribution of the phase space density
$\bar{f}\left(\bar{\gamma},\bar{\mu},\bar{\phi}\right)$ of accelerated
particles at the shock front ($\hat{x}=0$), as seen in the frame in which the
shock front is at rest. The downstream magnetic field is in the direction
$\bar{\mu}=0$, $\bar{\phi}=0$, particles crossing into the upstream region
have $\bar{\mu}>0$.

For $\eta=\eta_{\rm d}\gg 1$ and $\Gamma=\Gamma_{\rm d}\sim 1$, equation (1)
reverts to Liouville’s equation, albeit written in our unconventional mixed
coordinate system. The distribution of particles that are too energetic to
experience scattering downstream is, therefore, controlled by Liouville’s
theorem, which dictates that the phase-space density remains constant along
particle trajectories. When viewed from a frame of reference in which the
downstream plasma is at rest these trajectories are simply helices along which
the Lorentz factor $\gamma$ remains constant.

At the shock front, $x=0$, the domain of momentum space $A$ can be divided
into three non-overlapping sub-domains: (i) those trajectories crossing into
downstream that subsequently return to the shock front, denoted by $R$ (ii)
those crossing into downstream that subsequently escape without re-
encountering the shock, denoted by $E$, and (iii) those crossing into
upstream, (i. e., all trajectories with $0<\xi<1$) denoted by $U$. Given the
coordinates in momentum space $\vec{p}=(\gamma,\xi,\phi)$ of a trajectory in
$R$, it is straightforward to compute the mapping ${\cal M}:R\rightarrow U$
that relates the point $\vec{p}$ on a particle trajectory to the point
$\vec{p}_{+}=\left(\gamma_{+},\xi_{+},\phi_{+}\right)={\cal M}\vec{p}$ at
which this trajectory returns to the shock front. In Figure 1, we provide a
schematic sketch of a returning trajectory (for details see Appendix B). Then,
Liouville’s theorem requires the distribution at the shock front, $x=0$ to
satisfy

$\displaystyle f\left(\vec{p},0\right)$ $\displaystyle=f\left({\cal
M}\vec{p},0\right)\qquad\vec{p}\in R.$ (9)

Conversely, given the coordinates of a trajectory that crosses the shock from
downstream to upstream, one can simply invert this mapping to find the
coordinates $\vec{p}_{-}=\left(\gamma_{-},\xi_{-},\phi_{-}\right)={\cal
M}^{-1}\vec{p}$ with which it previously entered the downstream region, which
implies

$\displaystyle f\left(\vec{p},0\right)$ $\displaystyle=f\left({\cal
M}^{-1}\vec{p},0\right)\qquad\vec{p}\in U\,.$ (10)

Liouville’s theorem does not provide a constraint on the points $\vec{p}\in E$
on trajectories that escape from the shock into the downstream plasma, but it
is convenient to define the mapping ${\cal M}$ to be unity when operating on
these: ${\cal M}\vec{p}=\vec{p}$, $\vec{p}\in E$.

###  2.4 Approximation scheme

Since the problem, as formulated here, does not contain a characteristic scale
for the particle Lorentz factor $\gamma$, we look for solutions that are a
power law in this quantity, $F\left(\gamma\right)=\gamma^{-s}$. The index $s$
is then determined by matching the phase-space density across the shock front,
i. e., by imposing the conditions (9) and (10) on the expression (6) for $f$
in the upstream region, evaluated at the shock front $\hat{x}=0$.

However, to find a numerical value for $s$, it is necessary to implement an
approximation scheme. Here, we use a variant of the Galerkin method, similar
to that used by Kirk & Schneider (1987), which essentially truncates the
expansion (6) after the first $i_{\rm max}$ terms. Writing

$\displaystyle f\left(\vec{p},0\right)$
$\displaystyle=\gamma^{-s}\left[g_{i_{\rm max}}\left(\xi,\phi\right)+{\cal
R}_{i_{\rm max}}\left(\xi,\phi\right)\right]\,,$ (11) where $\displaystyle
g_{i_{\rm max}}\left(\xi,\phi\right)$ $\displaystyle=\sum_{i=1}^{i_{\rm
max}}a_{i}Q_{i}\left(\xi,\phi\right)\,$ (12)

is the desired approximation to the angular part of the distribution function
at the shock front, we demand that the residual ${\cal R}_{i_{\rm max}}$ be
orthogonal to the first $i_{\rm max}$ eigenfunctions:

$\displaystyle\iint_{\vec{p}\in A}\textrm{d}\xi\textrm{d}\phi
Q_{i}\left(\xi-1\right){\cal R}_{i_{\rm max}}\left(\xi,\phi\right)$
$\displaystyle=0\ \textrm{ for }i=1,\dots i_{\rm max}$ (13)

and, additionally, that the constraints (9) and (10) are satisfied not only by
$f$, but also by its approximation, $\gamma^{-s}g_{i_{\rm max}}$. This scheme
can be motivated physically if the terms in the summation in eq (6) are
ordered by the eigenvalue $\Lambda_{i}$, such that

$\displaystyle\Lambda_{i}$ $\displaystyle\geq\Lambda_{i+1},\qquad i\geq 1.$
(14)

Since $\Lambda_{1}<0$, the higher order terms in the expansion decay ever more
rapidly with increasing distance from the shock in the upstream region. The
particles described by these terms move almost in the plane of the shock, and,
as a consequence, receive only a small boost in energy in a cycle of crossing
and re-crossing. It is, therefore, to be expected that the index $s$ is
determined primarily by the first few terms in the summation.

Multiplying equation (11) by $\gamma^{s}Q_{j}$ and the weighting function, and
integrating over $\vec{p}\in A$ leads to

$\displaystyle\sum_{j=1}^{i_{\rm max}}\left[D_{ij}-M_{ij}(s)\right]a_{j}$
$\displaystyle=0\qquad i=1,\dots i_{\rm max}\,,$ (15)

where the diagonal matrix $D_{ij}$ results from substituting the expansion (6)
into the left-hand side of (11):

$\displaystyle D_{ij}$ $\displaystyle=\iint_{\vec{p}\in
A}\\!\\!\textrm{d}\xi\textrm{d}\phi\,Q_{i}\left(\xi,\phi\right)\left(\xi-1\right)Q_{j}\left(\xi,\phi\right)\,.$
(16)

Similarly, the matrix $M_{ij}(s)$ is obtained from the right-hand side by
applying the constraints (9) and (10) to the function $g_{i_{\rm max}}$ and
the constraint (13) to the function ${\cal R}_{i_{\rm max}}$:

$\displaystyle M_{ij}(s)$ $\displaystyle=\gamma^{s}\iint_{\vec{p}\notin
U}\\!\\!\textrm{d}\xi\textrm{d}\phi\,\gamma_{+}^{-s}Q_{i}\left(\xi,\phi\right)\left(\xi-1\right)Q_{j}\left(\xi_{+},\phi_{+}\right)$
$\displaystyle+\gamma^{s}\iint_{\vec{p}\in
U}\\!\\!\textrm{d}\xi\textrm{d}\phi\,\gamma_{-}^{-s}Q_{i}\left(\xi,\phi\right)\left(\xi-1\right)Q_{j}\left(\xi_{-},\phi_{-}\right)\,.$
(17)

Thus, eq (15) represents a system of homogeneous, linear algebraic equations
for the $a_{j}$, which have a non-trivial solution only if

$\displaystyle\textrm{Det}\left[D_{ij}-M_{ij}(s)\right]$ $\displaystyle=0\,.$
(18)

Given a guess for $s$, it is straightforward to evaluate these matrices by
numerical quadrature. Then, a root-finder algorithm applied to equation (18)
yields $s$ and the coefficients $a_{i}$ are found from the null-space of the
corresponding matrix.

##  3 Results

The usefulness of the approximation scheme described in 2.4 is confirmed by
the rapid convergence as the number of eigenfunctions $i_{\rm max}$ is
increased, as shown in fig 2. The converged value of $s$ for an
ultrarelativistic shock front in an ideal fluid, for which $\beta_{\rm
d}=1/3$, is $4.17$, close to, but slightly harder than the result found when
particle transport is dominated by scattering both upstream and downstream.
This value depends only on $\beta_{\rm d}$, and is, in particular, independent
of both the strength of the upstream scattering and the strength of the
downstream magnetic field. As shown in panel (a) of Fig 2, convergence to this
result requires only a few ($\sim 4$) eigenfunctions. Because of the ordering
of the eigenfunctions, equation (14), convergence of the angular dependence of
the phase space density is slowest at the shock front itself. Panels (b), (c)
and (d) show this dependence as seen in the frame in which the shock is at
rest, and the flow is directed along the shock normal. In this reference
frame, the transformed spherical polar coordinates are denoted by
$\bar{\gamma},\bar{\mu},\bar{\phi}$ where

$\displaystyle\bar{\gamma}$ $\displaystyle=\gamma\Gamma_{\rm
s}\left(1-\beta_{\rm s}\mu\right)$ (19) $\displaystyle\bar{\mu}$
$\displaystyle=\left(\mu-\beta_{\rm s}\right)/\left(1-\mu\beta_{\rm s}\right)$
(20) $\displaystyle\bar{\phi}$ $\displaystyle=\phi$ (21)

and, since the phase-space density, denoted in this frame by $\bar{f}$, is a
Lorentz invariant quantity

$\displaystyle\bar{f}\left(\bar{\gamma},\bar{\mu},\bar{\phi},x\right)$
$\displaystyle=f\left(\gamma,\mu,\phi,x\right)\,.$ (22)

Note that although the average over either $\bar{\phi}$ or $\bar{\mu}$ of the
distribution converges with only $\sim 3$ and $\sim 9$ eigenfunctions,
respectively, the phase distribution of those particles that move precisely
along the shock front converges much more slowly, needing $\sim 15$
eigenfunctions. This illustrates the fact that grazing particles have
essentially no impact on the power law index.

Figure 4:  A comparison of the angular distributions just upstream of the
shock front found by the eigenfunction method in the ultrarelativistic limit,
and by Monte-Carlo simulation for $\Gamma_{\rm s}=50$. In (a) the
distributions averaged over phase (about the shock normal), and, in (b), over
the angle to the shock normal are compared. In each case the result for 30
eigenfunctions is shown as a yellow curve and that of the simulation as a
purple histogram. In (c) and (d) the full distributions are shown for the
eigenfunction and Monte Carlo methods, respectively.  Figure 5:  The angle-
averaged distribution just upstream at $\hat{x}=0.1$, as a function of the
particle Lorentz factor $\bar{\gamma}$, measured in the rest frame of the
shock. As the shock Lorentz factor $\Gamma_{\rm s}$ increases, the power-law
section of the distribution converges towards the result found by the
eigenfunction method in the ultrarelativistic limit, shown here as a
horizontal line. For each $\Gamma_{\rm s}$, the power-law index $s$, found
from a least-squares fit to the data in the region
$0.5<\bar{\gamma}/\Gamma_{\rm s}<3$, is listed in the legend.

The angular distribution at the shock front is, as expected, anisotropic. In
addition to the anisotropy with respect to the shock normal, which is well-
known from earlier studies of the scattering-dominated case and which arises
because the relative velocity of the upstream and downstream plasmas is only
slightly less than the particle velocity, a strong anisotropy in the azimuthal
angle $\bar{\phi}$ is present because the downstream magnetic field imposes a
drift along the shock front. This is clearly seen in the full 2-D angular
distribution at the shock front, as shown Fig. 3. In terms of a right-handed
system of coordinates with the magnetic field along the positive $z$-axis, the
drift for a positively charged particle is in the positive $y$-direction
($\sin\bar{\phi}<0$), i. e., in the direction opposite to that of the electric
field seen in the shock rest frame. Thus, as seen in the rest frame of the
shock, particles gain energy during an excursion upstream, but lose some on
their downstream loop.

Monte-Carlo simulations, using the code described in ZH , provide an
independent cross-check of these results, as well as extending them by lifting
the restriction to the ultrarelativistic limit. In Fig 4, we compare the
angular distributions found by both methods. Shock grazing particles are
problematic in the Monte-Carlo approach when sampling the distribution
precisely on the shock surface. Making the comparison a short distance
upstream, here at the surface $\hat{x}=0.1$, mitigates the problem. Results
for the eigenfunction method are shown in the ultrarelativistic limit, since
using a finite value of $\Gamma_{\rm s}$ would necessitate a numerical
evaluation of the eigenfunctions. The Monte-Carlo simulations were performed
with $\Gamma_{\rm s}=50$. Agreement is generally excellent, with small
fluctuations visible only in the Monte-Carlo results for the full 2-D
distribution.

The distribution of accelerated particles in energy for shock speeds ranging
from mildly to highly relativistic is shown in Fig 5. These are found using
Monte-Carlo simulations that inject $10^{6}$ particles into the upstream with
$\gamma=2\Gamma_{\rm s}^{2}$, and a uniform angular distribution within a cone
about the shock normal of opening angle $1/\Gamma_{\rm s}$. The trajectories
are then followed until they escape downstream, whilst registering the value
of $\bar{\gamma}$ at each crossing of the shock. After several crossings, the
distribution settles into a power law, that extends up to a point at which the
statistical noise becomes significant. For simplicity, we choose in each case
the jump conditions for a relativistic gas: $\beta_{\rm d}\beta_{\rm s}=1/3$.
The figure shows the distribution weighted by the factor $\bar{\gamma}^{4.17}$
in order to highlight the departure of the finite $\Gamma_{\rm s}$ results
from the ultrarelativistic result found by the eigenfunction method. It can be
seen that the power-law index $s$ is within a few percent of its asymptotic
value for $\Gamma_{\rm s}>5$.

##  4 Discussion

Whether or not the mechanism of diffusive shock acceleration operates
effectively at a perpendicular shock is a question that is still the subject
of controversy, over four decades after the publication of the discovery
papers, which implicitly addressed parallel shocks. In this context
“effectively” means either that the acceleration rate is comparable to or
faster than that at a parallel shock, or that the power-law index of the
stationary phase-space distribution of accelerated particles is close to or
harder than that produced at a parallel shock. The main result of this paper
concerns the power-law index produced by highly relativistic, and, therefore,
generically perpendicular shocks. Assuming that particles can be treated as
unmagnetized when upstream of the shock (4), but as magnetized when downstream
(8), we demonstrate quantitatively that these shocks are just as effective
accelerators as the possibly less realistic parallel shocks addressed by
previous analytical work. This result has major implications also for the
expected acceleration rate, and the related maximum energy to which particles
can be accelerated at a shock in a given physical situation, questions which
are addressed in ZH .

The persistence of turbulence in the downstream medium at a level sufficient
to demagnetize particles is generally perceived to be a major problem for the
theory of diffusive shock acceleration at relativistic shocks (e. g., Bykov et
al., 2012). Our assumption that turbulence is completely negligible downstream
is specifically designed to address this point. It is, therefore, not
restrictive, in the sense that diffusive acceleration can be expected to
proceed as previously predicted if this assumption is not justified in a
particular application.

On the other hand, the assumption that particles diffuse in angle when
upstream of the shock is important. Unless a degree of randomness enters into
the trajectories that return to the shock from upstream, diffusive
acceleration will cease and particles will receive only a single, finite boost
in energy before being swept away downstream (Begelman & Kirk, 1990; Pelletier
et al., 2009). We have adopted the simple prescription of isotropic diffusion
to describe this randomness, since it renders the analytic approximations
tractable; at least in the parallel shock case, this assumption does not
appear to be crucial (Kirk et al., 2000). However, the validity of our
approach does depend on the presence of turbulence of sufficient amplitude in
the upstream medium. In the case of a weakly magnetized plasma $\sigma_{\rm
reg}<10^{-3}$, such as might be encountered by the blast wave of a gamma-ray
burst (GRB), a lower limit on the level of upstream turbulence is provided by
shock-generated Weibel filaments, which, according to particle-in-cell
simulations and analytical considerations, are amplified to $\sigma_{\rm
turb}\approx 0.1$. According to (2), this is already sufficient to demagnetize
particles with $\gamma<0.1\times\Gamma_{\rm s}\left({m_{\rm
p}/m}\right)\sigma_{\rm reg}^{-1/2}$, and, therefore, enable diffusive
acceleration even if the filaments are strongly damped downstream. However,
this limit is unduly restrictive if turbulence exists in the upstream that is
not directly excited by the processes that form the shock. In the case of a
GRB, for example, the ionization, heating and pair-loading of the surrounding
medium by the prompt emission (Beloborodov, 2002; Grošelj et al., 2022) seems
unlikely to leave behind a quiescent environment. Furthermore, turbulence is
known to be present in the winds of the progenitors of some supernovae of the
type Ibc (Wellons et al., 2012) that is associated with long duration GRBs. On
the other hand, if the upstream medium is strongly magnetized, such as in the
case of the termination shock of a pulsar wind, turbulence is embedded in the
outflow by the pulsar and is thought to facilitate acceleration (Sironi &
Spitkovsky, 2011; Giacchè & Kirk, 2017). Whereas previous discussions assumed
this process to be confined to the equatorial region of the wind, where the
regular field component vanishes, the results presented above suggest that it
may persist to higher latitudes, with interesting implications for modelling
the emission from pulsar wind nebulae (Olmi et al., 2015).

In summary, if sufficient turbulence is present in the upstream medium, our
results demonstrate that the perpendicular magnetic field in the plasma
downstream of a relativistic shock front does not inhibit acceleration. Rather
than being swept away without returning to the shock front, particles indeed
return, and populate a power-law distribution whose index is insensitive to
the nature of the dominant transport mechanism downstream.

## Acknowledgements

We thank Nils Schween, Gwenael Giacinti, Tony Bell, Martin Lemoine and Daniel
Grošelj for helpful discussions.

## Data availability

No data relevant to this paper are available.

## References

  * Achterberg et al. (2001) Achterberg A., Gallant Y. A., Kirk J. G., Guthmann A. W., 2001, MNRAS, 328, 393
  * Begelman & Kirk (1990) Begelman M. C., Kirk J. G., 1990, ApJ, 353, 66
  * Bell (2014) Bell A. R., 2014, Brazilian Journal of Physics, 44, 415
  * Beloborodov (2002) Beloborodov A. M., 2002, ApJ, 565, 808
  * Bykov et al. (2012) Bykov A., Gehrels N., Krawczynski H., Lemoine M., Pelletier G., Pohl M., 2012, Space Sci. Rev., 173, 309
  * Chang et al. (2008) Chang P., Spitkovsky A., Arons J., 2008, ApJ, 674, 378
  * Giacchè & Kirk (2017) Giacchè S., Kirk J. G., 2017, ApJ, 835, 235
  * Grošelj et al. (2022) Grošelj D., Sironi L., Beloborodov A. M., 2022, ApJ, 933, 74
  * Huang et al. (2022) Huang Z.-Q., Kirk J. G., Giacinti G., Reville B., 2022, ApJ, 925, 182
  * Huang et al. (2023) Huang Z.-Q., Reville B., Kirk J. G., Giacinti G., 2023, to be submitted to MNRAS
  * Ince (1956) Ince E. L., 1956, Ordinary Differential Equations. Dover Publications, Inc, New York
  * Keshet et al. (2009) Keshet U., Katz B., Spitkovsky A., Waxman E., 2009, ApJ, 693, L127
  * Kirk & Schneider (1987) Kirk J. G., Schneider P., 1987, ApJ, 315, 425
  * Kirk & Schneider (1989) Kirk J. G., Schneider P., 1989, A&A, 225, 559
  * Kirk et al. (2000) Kirk J. G., Guthmann A. W., Gallant Y. A., Achterberg A., 2000, ApJ, 542, 235
  * Lemoine (2015) Lemoine M., 2015, Journal of Plasma Physics, 81, 455810101
  * Lemoine & Pelletier (2010) Lemoine M., Pelletier G., 2010, MNRAS, 402, 321
  * Olmi et al. (2015) Olmi B., Del Zanna L., Amato E., Bucciantini N., 2015, MNRAS, 449, 3149
  * Pelletier et al. (2009) Pelletier G., Lemoine M., Marcowith A., 2009, MNRAS, 393, 587
  * Plotnikov et al. (2018) Plotnikov I., Grassi A., Grech M., 2018, MNRAS, 477, 5238
  * Reville & Bell (2014) Reville B., Bell A. R., 2014, MNRAS, 439, 2050
  * Sironi & Spitkovsky (2009) Sironi L., Spitkovsky A., 2009, ApJ, 698, 1523
  * Sironi & Spitkovsky (2011) Sironi L., Spitkovsky A., 2011, ApJ, 741, 39
  * Sironi et al. (2013) Sironi L., Spitkovsky A., Arons J., 2013, ApJ, 771, 54
  * Sironi et al. (2015) Sironi L., Keshet U., Lemoine M., 2015, Space Sci. Rev., 191, 519
  * Takamoto & Kirk (2015) Takamoto M., Kirk J. G., 2015, ApJ, 809, 29
  * Vanthieghem et al. (2020) Vanthieghem A., Lemoine M., Plotnikov I., Grassi A., Grech M., Gremillet L., Pelletier G., 2020, Galaxies, 8, 33
  * Wellons et al. (2012) Wellons S., Soderberg A. M., Chevalier R. A., 2012, ApJ, 752, 17

##  Appendix A 2D eigenfunctions

Inserting the expansion (6) into (5) and separating the variables according to

$\displaystyle Q_{i}\left(\xi,\phi\right)$
$\displaystyle=T_{i}\left(\xi\right)S_{i}\left(\phi\right)$ (23)

one obtains for the $\phi$-dependent function

$\displaystyle S^{\prime\prime}_{i}$ $\displaystyle=-j^{2}S_{i}\,,$ (24)

where $j$ is a constant, which, since $S$ is periodic with period $2\pi$, is
an integer. Without loss of generality, we can choose $j\geq 0$ and identify
two families of solutions that satisfy the additional symmetry
$S_{i}(\phi)=S_{i}(\pi-\phi)$:

$\displaystyle S_{i}\propto\left\\{\begin{array}[]{ll}\cos(j\phi)&j\textrm{
even or zero}\\\ \sin(j\phi)&j\textrm{ odd}\,.\end{array}\right.$ (27)

The equation determining the $\xi$-dependent eigenfunction is

$\displaystyle\left(\xi
T_{i}^{\prime}\right)^{\prime}-\left[\frac{j^{2}}{4\xi}+\Lambda_{i}\left(1-\xi\right)\right]T_{i}$
$\displaystyle=0\,.$ (28)

Following Kirk & Schneider (1989, appendix A) (see also Ince (1956, section
7.31)), one looks for a solution of the form

$\displaystyle T_{i}$
$\displaystyle=\textrm{e}^{\lambda\xi}\xi^{\alpha}\sum_{n=0}^{\infty}c_{n}\xi^{n}\,.$
(29)

Inserting this into eq (28), shows that the choice
$\lambda=-\sqrt{-\Lambda_{i}}$ is convenient, since it removes the highest
power of $\xi$ in the term proportional to $T_{i}$, leading to a two-term
recurrence relation for the $c_{n}$. The indicial equation is
$\alpha^{2}=j^{2}/4$, so that there are two possible solutions

$\displaystyle T^{+}_{i}$
$\displaystyle=\textrm{e}^{-\sqrt{-\Lambda}_{i}\xi}\sum_{n=0}^{\infty}\xi^{n}c^{+}_{n}$
$\displaystyle j\textrm{ even or zero}$ (30) $\displaystyle T^{-}_{i}$
$\displaystyle=\textrm{e}^{-\sqrt{-\Lambda}_{i}\xi}\sum_{n=0}^{\infty}\xi^{n+\frac{1}{2}}c^{-}_{n}$
$\displaystyle j\textrm{ odd.}$

Inserting (30) into (28) and equating coefficients gives $j^{2}c^{+}_{0}=0$
and $\left(j^{2}-1\right)c^{-}_{0}=0$ and, for $n\geq 0$, the recurrence
relations

$\displaystyle c^{+}_{n+1}$
$\displaystyle=\left(\frac{4\sqrt{-\Lambda_{i}}\left(2n+1-\sqrt{-\Lambda_{i}}\right)}{(2n+2)^{2}-j^{2}}\right)c^{+}_{n}$
(31) $\displaystyle c^{-}_{n+1}$
$\displaystyle=\left(\frac{4\sqrt{-\Lambda_{i}}\left(2n+2-\sqrt{-\Lambda_{i}}\right)}{(2n+3)^{2}-j^{2}}\right)c^{-}_{n}\,.$

Since $c^{\pm}_{n+1}/c^{\pm}_{n}\rightarrow 2\sqrt{-\Lambda_{i}}/n$, as
$n\rightarrow\infty$ it follows that
$T^{\pm}_{i}\rightarrow\textrm{e}^{\sqrt{-\Lambda_{i}}\xi}$ as
$\xi\rightarrow\infty$, unless the relevant series truncates at finite $n$.
Therefore, the eigenvalues satisfying the boundary conditions of boundedness
at $\xi=0,\infty$ are found by requiring truncation for, say, $n>k\geq 0$.
Then,

$\displaystyle\Lambda_{i}$
$\displaystyle=\left\\{\begin{array}[]{ll}-(2k+1)^{2}&j\textrm{ even or
zero},0\leq j\leq 2k\\\ &\\\ -4(k+1)^{2}&j\textrm{ odd},1\leq j\leq
2k+1\,,\end{array}\right.$ (35)

where $k$ is a positive integer or zero. Thus, each index $i$ corresponds to a
pair of indices $k$ and $j$. The eigenvalues $\Lambda_{i}$ depend on $k$ and
only the parity of $j$. The corresponding eigenfunctions can be evaluated
using the recurrence relations (31). The first few are listed in Table 1.

Table 1: 2D-eigenfunctions: $i$ uniquely identifies the eigenfunction $Q_{i}$, $j$ defines its $\phi$-dependence according to eq (27), $k$ is the largest integer such that $c_{k}\neq 0$ in eq (30), $\Lambda_{i}$ is the eigenvalue determining the $x$-dependence associated with this eigenfunction in eq (6). The (arbitrary) normalisation of the $Q_{i}$ is $c_{n}=1$, where $n={j/2}$ ($j$ even or zero), or $n={(j-1)/2}$ ($j$ odd) is the smallest integer for which $c_{n}\neq 0$. $i$ | $j$ | $k$ | $\Lambda_{i}$ | $Q_{i}$
---|---|---|---|---
$1$ | $0$ | $0$ | $-1$ | $\textrm{e}^{-{\xi}}$
$2$ | $1$ | $0$ | $-4$ | $\textrm{e}^{-2\xi}\sqrt{\xi}\sin(\phi)$
$3$ | $0$ | $1$ | $-9$ | $\textrm{e}^{-3{\xi}}\left(1-6{\xi}\right)$
$4$ | $2$ | $1$ | $-9$ | $\textrm{e}^{-3{\xi}}{\xi}\cos(2\phi)$
$5$ | $1$ | $1$ | $-16$ | $\textrm{e}^{-4{\xi}}\sqrt{{\xi}}\left(1-4{\xi}\right)\sin\phi$
$6$ | $3$ | $1$ | $-16$ | $\textrm{e}^{-4{\xi}}{\xi}^{3/2}\sin(3\phi)$

##  Appendix B Mapping

For an ultrarelativistic particle ($\gamma\gg\Gamma_{\rm s}$), membership of
the three subdomains of momentum space at the shock front — upstream ($U$),
returning ($R$) and escaping ($E$) — of the particle phase space ($A$) is
determined solely by the direction of motion, labelled by the stretched
variable $\xi$, defined in (3), and the azimuthal phase (with respect to the
shock normal) $\phi$. Membership of $U$ requires only $\xi<1$, but the
requirements for $R$ and $E$ are more complicated. First, define the auxiliary
parameter $\zeta$ to be the ratio of the shock speed $c\beta_{\rm d}$ to the
component $c\beta_{\bot}$ of the particle speed that is perpendicular to the
magnetic field, both seen from the frame in which the downstream plasma is at
rest. The computation of $\zeta$ from $(\xi,\phi)$ involves a Lorentz boost
from the rest frame of the upstream plasma to that of the downstream plasma,
followed by a change of axis of the polar coordinates from the shock normal to
the magnetic field. To lowest order in $1/\Gamma_{\rm s}$ one finds:

$\displaystyle\zeta$ $\displaystyle=\frac{\beta_{\rm
d}\left[\left(1+\beta_{\rm d}\right)+\left(1-\beta_{\rm
d}\right)\xi\right]}{\sqrt{\left[\left(1+\beta_{\rm
d}\right)+\left(1-\beta_{\rm d}\right)\xi\right]^{2}-4\xi\left(1-\beta_{\rm
d}^{2}\right)\cos^{2}\phi}}\,\geq\,\beta_{\rm d}\,.$ (36)

This quantity remains constant during the particle’s residence downstream.
Clearly, $\zeta>1$ implies that the particle cannot recross the shock front
($\vec{p}\in E$) and the mapping ${\cal M}$ is unity:

$\displaystyle{\cal M}\vec{p}=\vec{p_{+}},\ \vec{p}\in E$ $\displaystyle:$
$\displaystyle\gamma_{+}=\gamma\qquad\xi_{+}=\xi\qquad\phi_{+}=\phi\,.$ (37)

However, although $\zeta>1$ is sufficient for $\vec{p}\in E$, it is not a
necessary condition. The distance $d$, in units of the particle’s gyroradius,
between the shock and a point on the trajectory depends on $\zeta$, the
particle’s phase $\Phi$ (measured about the downstream magnetic field, such
that it is an increasing function of time for a positively charged particle)
and the phase $\Phi_{0}$ at which the trajectory intersects the shock front:

$\displaystyle d\left(\zeta,\Phi,\Phi_{0}\right)$
$\displaystyle=\zeta\left(\Phi-\Phi_{0}\right)-\left(\sin\Phi-\sin\Phi_{0}\right)\,.$
(38)

From the Lorentz boost and coordinate transformation, one finds, to lowest
order in $1/\Gamma_{\rm s}$,

$\displaystyle\Phi_{0}$
$\displaystyle=\textrm{atan2}\left[2\sin\phi\sqrt{\left(1-\beta_{\rm
d}^{2}\right)\xi},\,\left(1+\beta_{\rm d}\right)-\left(1-\beta_{\rm
d}\right)\xi\right]+2n\pi\,,$ (39)

where, to simplify the discussion, the integer $n$ is to be chosen such that
$-2\pi+\textrm{arccos}\zeta\leq\Phi_{0}<\textrm{arccos}\zeta$. For
$\vec{p}\notin U$, the distance $d$ grows initially, implying
$-2\pi+\textrm{arccos}\zeta<\Phi_{0}<-\textrm{arccos}\zeta$, and subsequently
goes through an alternating series of maxima and minima as the trajectory
gyrates about the magnetic field. The first maximum is reached when
$\Phi=-\textrm{arccos}\zeta$ and the subsequent minimum when
$\Phi=+\textrm{arccos}\zeta$. Therefore, if $\zeta<1$, a sufficient condition
for the particle to escape is that $d>0$ at this minimum, i. e.,

$\displaystyle d\left(\zeta,\textrm{arccos}\zeta,\Phi_{0}\right)$
$\displaystyle>0\,,$ (40)

in which case $\cal M$ is again the unit mapping (37). If, on the other hand,
$d<0$ at this point, then $\vec{p}\in R$, and the phase $\Phi_{+}$, with which
the trajectory returns to the shock front, is given by the single solution of
the equation

$\displaystyle d\left(\zeta,\Phi_{+},\Phi_{0}\right)$ $\displaystyle=0\,,$
(41)

with $-\textrm{arccos}\zeta<\Phi_{+}<\textrm{arccos}\zeta$. Similarly, for
$\vec{p}\in U$, which implies
$-\textrm{arccos}\zeta<\Phi_{0}<\textrm{arccos}\zeta$, the phase $\Phi_{-}$,
at which the trajectory previously entered the downstream region, is the
single solution of the equation

$\displaystyle d\left(\zeta,\Phi_{0},\Phi_{-}\right)$ $\displaystyle=0\,,$
(42)

with $-2\pi+\textrm{arccos}\zeta<\Phi_{-}<-\textrm{arccos}\zeta$. Thus, the
boundaries between the different subdomains for particles in $R$, $E$ and $U$
are uniquely defined in the $\Phi_{0}-\zeta$ plane, as summarised in Figure 6.

To find ${\cal M}$, it remains to transform the coordinates and apply a
Lorentz boost back to the upstream frame, which gives, to lowest order in
$1/\Gamma_{\rm s}$,

$\displaystyle{\cal M}\vec{p}=\vec{p_{+}},\ \vec{p}\in R$
$\displaystyle\textrm{ and }{\cal M}^{-1}\vec{p}=\vec{p_{-}},\ \vec{p}\in U:$
$\displaystyle\gamma_{\pm}=\gamma\left(\zeta+\beta_{\rm
d}\cos\Phi_{\pm}\right)\frac{\left(1-\beta_{\rm
d}\right)\xi+\left(1+\beta_{\rm d}\right)}{2\left(1-\beta_{\rm
d}\right)\zeta}$ $\displaystyle\xi_{\pm}=\frac{\left(1+\beta_{\rm
d}\right)\left(\zeta-\beta_{\rm d}\cos\Phi_{\pm}\right)}{\left(1-\beta_{\rm
d}\right)\left(\zeta+\beta_{\rm d}\cos\Phi_{\pm}\right)}$
$\displaystyle\phi_{\pm}=\textrm{atan2}\left(\beta_{\rm
d}\sin\Phi_{\pm},\sqrt{\zeta^{2}-\beta_{\rm d}^{2}}\right)\,.$ (43)

Figure 6: The three sub-domains of particle phase space at the shock front in
the $\Phi_{0}-\zeta$ plane: purple for $\vec{p}\in R$, green for $\vec{p}\in
U$ and the remaining unshaded area for $\vec{p}\in E$. Note that
$\zeta>\beta_{\rm d}$ .


# Effect of atomic anti-site disorder on the AMR in FeCo alloys

Mingsong Zhang  School of Integrated Circuits Science and Engineering,
University of Electronic Science and Technology of China, 610081, Chengdu,
P.R. China  Bin Peng  School of Integrated Circuits Science and Engineering,
University of Electronic Science and Technology of China, 610081, Chengdu,
P.R. China  Wanli Zhang  School of Integrated Circuits Science and
Engineering, University of Electronic Science and Technology of China, 610081,
Chengdu, P.R. China  Wenxu Zhang Corresponding Author<EMAIL_ADDRESS>School of Integrated Circuits Science and Engineering, University of
Electronic Science and Technology of China, 610081, Chengdu, P.R. China

###### Abstract

In order to understand the anti-site disorder effect on the anisotropic
magnetoresistance (AMR) effect in alloys, $\rm{Fe}_{50}Co_{50}$ alloys were
studied in this work using the fully relativistic spin-polarized screened
(KKR) method. The anti-site effect was modeled by interchanging Fe and Co
atoms and treated by the coherent potential approximation (CPA). We find that
the anti-site disorder broadens the spectral function and decreases the
conductivity. Our work emphasizes that the absolute variations of resistivity
under magnetic moment rotation are less affected by atomic disorders. The
annealing procedure improves the AMR by reduction of the total resistivity. At
the same time, we also find that the fourth-order term in the angular
dependent resistivity becomes weaker when the disorder increases, resulting
from increased scattering of the states around the band-crossing.

##  1 Introduction

The anisotropic magnetoresistance (AMR) effect, where the longitudinal
resistivity in ferromagnetic alloy is dependent on the magnetization
direction, was discovered by Thomson more than one hundred years ago [1]. Its
applications include advanced information storage and modern sensors [2, 3].
However, unrevealing its mechanism has been pondered [1]. To explain the AMR,
Fert and Campbell proposed the two-current model [4], where spin-dependent
scattering plays the central role. Residual resistivity and AMR of bulk NiFe
alloy were calculated using the Kubo-Greenwood formula with the fully
relativistic spin-polarized screened (KKR) method and the coherent potential
approximation (CPA) [5, 6]. Although the calculated residual resistivities are
30~40% smaller than the measured values, the AMR ratios were in excellent
agreement with the experiments in the whole composition percentage. Very
recently, the topological variations of the bands with the magnetization
directions were used to account for the AMR. The spin-orbital coupling (SOC)
dependent gap between the anticrossing bands in the vicinity of the Fermi
energy was the cause of the changing resistivity with respect to the
magnetization direction [7]. At the same time, the fourfold symmetric term of
AMR was attributed to the relaxation time anisotropy due to the variation of
the density of states near the Fermi energy under the rotation of the
magnetization [8]. The fourfold AMR was also proposed to originate from the
splitting of the 3d states of the impurities owing to the SOI [9].

Nowadays, more details are included in the computation of AMR. The AMR at
finite temperatures in hexagonal late transition metals and their alloys was
calculated [10, 11], where phonon and spin fluctuation scattering were taken
into account. The combination of both reproduced the experimental resistivity.
Disorder can not only be introduced by phonon or magnon at finite temperatures
but also by the atom occupations at the lattice. Atomic disorder is widely
observed in alloys and compounds, especially materials with a B2 lattice or
sublattice [12]. In these materials, the measured AMR decreases with the
annealing temperatures. Our previous work shows that the change in AMR is
mainly due to the reduction of the total resistivity [13]. It was usually
attributed to the reduction of the grain boundary and the growth of the grain
size [14, 15]. The structure disorders can also be reflected in anisotropic
magnetoresistance [16]. There are always some complexities in explaining the
AMR influenced by disorder. Disorder may increase the resistivity by enhancing
the electron scattering and smearing the band structure. At the same time, the
gaps formed in the anticrossing bands mentioned above will also change. In
this work, we find that anti-site disorder induced mainly the broadening of
the spectral function and the decrease in the conductivity, which is the main
reason for the decrement of the AMR. In contrast, the difference in
conductivity introduced by the rotation of the magnetization direction is not
influenced by the disorder. Namely, the intrinsic origins of AMR are robust
against disorder. We also find that the disorder increases the symmetry and
leads to the decrement of the fourth-order term in AMR attached to the cubic
symmetry.

The paper is organized as follows: In Sec. II, we present the details related
to the calculation method and the parameters we used. In Sec. III, we analyze
the factors affecting electronic transport based on the electronic structure
and the Bloch spectral function. The paper is summarized in Sec. IV.

##  2 Computational details

In order to model the anti-site disorder in BCC FeCo alloy, the concentration
of Fe at (0, 0, 0) was set to be x, with the concentrations of
Fe(0.5,0.5,0.5), Co(0,0,0), and Co(0.5,0.5,0.5) were set to be 1-x, 1-x, and
x, respectively. We defined the configuration entropy as
$S=-k_{B}[x\ln{x}+(1-x)\ln{(1-x)}]$ to characterize the degree of anti-site
disorder where $k_{B}$ is the Boltzmann constant.

The electronic structure calculation of Fex(Co)1-xCox(Fe)1-x has been
implemented within the SPR-KKR (spin-polarized relativistic Korringa-Kohn-
Rostoker) program package [5]. The anti-site disorder was treated within the
coherent potential approximation (CPA) method [17]. We adopted the fully
relativistic mode (REL) to include the relativistic effect, which introduces
the spin-orbit coupling. Self-consistent potentials were obtained employing
the Vosko-Wilk-Nussair (VWN) parametrization [18] for the exchange-correlation
functional in the local density approximation (LDA). The energy integration
was performed on a semicircular path in the complex plane (GRID = 5) with 32
points of the E-mesh. The atomic sphere approximation (ASA) was used for the
potentials. $l_{max}$, the angular momentum expansion cutoff parameter, plays
an essential role in the convergence of KKR calculations [19]. According to
the previous works [20, 21] and our validation, $l_{max}=3$ is sufficient for
$d$-systems. To obtain the converged calculation, we used 7.2 $\times$
$10^{7}$ k points in the full Brillouin zone. The Bloch spectral function
(BSF) is applied to visualize the electronic structure of alloys, which
represents the spectral weight distribution of electrons in energy-momentum
space [22].

The Kubo linear response formalism was adopted for dealing with electronic
transport in metallic systems. Based on the formalism, using the independent
electron approximation and the thermal limit (T = 0 K), one arrives at the
Kubo-Středa equation for the conductivity tensor [23]

$\begin{split}\sigma_{\mu v}&=\frac{\hbar}{4\pi
V}\operatorname{Tr}\langle\hat{J}_{\mu}(G^{+}-G^{-})\hat{J}_{v}G^{-}-\hat{J}_{\mu}G^{+}\hat{J}_{v}(G^{+}-G^{-})\rangle_{c}\\\
&+\frac{e}{4\pi
iV}\operatorname{Tr}\langle(G^{+}-G^{-})(\hat{r}_{\mu}\hat{J}_{v}-\hat{r}_{v}\hat{J}_{\mu})\rangle_{c}.\end{split}$
(1)

Here $\hat{r}$ and $\hat{J}$ are the position and current density operators,
respectively. $\mu$ and $v$ denote the Cartesian coordinates. $V$ is the unit
cell volume. The retarded $\left(G^{+}\right)$ and advanced
$\left(G^{-}\right)$ Green functions represent the electronic structure of the
system. $\langle\cdots\rangle_{c}$ indicates a configurational average
performed within the CPA method.

##  3 Results and discussions

###  3.1 The lattice constant and the density of states

Figure 1: (a) The lattice constant decreases with the increase of $S$. The
inset is the relationship between $S$ and the concentration of Fe(0,0,0). (b)
The DOS projected on $\mathrm{Fe}$ and $\mathrm{Co}$ in
Fe0.5(Co)0.5Co0.5(Fe)0.5. (c) The densities of states of Fex(Co)1-xCox(Fe)1-x
alloys with different degrees of anti-site disorder.

The entropy (S) grows monotonically from 0 $k_{B}$ to 0.693 $k_{B}$ with the
increasing concentration of Fe(0,0,0). It reaches the maximum when the
concentrations of Fe and Co at the same site become identical. For this
reason, only the name of Fex(Co)1-xCox(Fe)1-x alloy is used to indicate the
degree of anti-site disorder in the subsequent part.

The lattice constant decreases with the increasing degree of disorder, the
same as the results of Kašpar et al. [24] and Jacob et al. [25]. Moreover,
Kašpar et al. suggested that the change in the lattice constant is due to some
occupied anti-bonding majority-spin molecular orbitals becoming partially
unoccupied in the transition from order to disorder [24]. Such a change in the
lattice constant will induce variations in the DOS as well as the magnetic
moment.

The average moment slightly depends on the anti-site disorder, rising from
2.11 $\mu_{B}$ to 2.13 $\mu_{B}$, in agreement with the conclusion of Spooner
et al. [26]. Because the Co atom is more electronegative and has a smaller
exchange splitting than Fe, the minority-spin states of Co lie at lower energy
than the Fe states [27]. Consequently, the DOS at the Fermi energy comes
mainly from the spin-down electrons of Co, as shown in FIG. 1 (b), so the
magnetic moment should primarily originate from cobalt.

Regarding the DOSs with different degrees of anti-site disorder, which are
exhibited in FIG. 1 (c), the ferromagnetically ordered alloys have only slight
changes in the DOS of the spin-up electrons. While for the spin-down parts,
substantial decreases in the peaks above the $E_{\rm{F}}$ smoothen the shape
of DOS, meaning a weaker localization. The sharp contrast between the DOS of
majority and minority electrons is consistent with the conclusions of Kulkova
et al. [27]. The flip of the majority-occupied states to the minority-
unoccupied ones will decrease the spin moment of the alloy, compensating for
the effect of localization.

For ferromagnet, the essential details lie in the minority spin system [28],
so it can be expected that the reduction of the DOS at the $E_{\rm{F}}$ due to
anti-site disorder will increase the resistivity of the alloys.

###  3.2 Resistivity and its angular dependence

Figure 2: (a) The definition of angle $\alpha$ in Cartesian coordinates. The
calculated AMR (b), diagonal resistivity $\rho_{xx}$ (c), and average
resistivity $\rho$ (d) of Fex(Co)1-xCox(Fe)1-x alloys. The curves in (c) and
(d) are fitted to Eq. (2) and (3), respectively.

Since the resistivity varies with the magnetization direction,
$\mathrm{AMR}=\frac{\left(\rho_{\|}-\rho_{\perp}\right)}{\left(\frac{1}{3}\rho_{\|}+\frac{2}{3}\rho_{\perp}\right)}\times
100\%$ is defined to characterize the anisotropic magnetoresistance, where
$\rho_{\|}\left(\rho_{\perp}\right)$ is the resistivity when the magnetization
is parallel (perpendicular) to the current direction. According to the
transformation law for the second rank tensors, we obtained
$\rho_{\|}=\cos^{2}\alpha\rho_{xx}+\cos\alpha\sin\alpha\rho_{xy}+\cos\alpha\sin\alpha\rho_{yx}+\sin^{2}\alpha\rho_{yy}$
and $\rho_{\perp}=$
$\sin^{2}\alpha\rho_{xx}-\cos\alpha\sin\alpha\rho_{xy}-\cos\alpha\sin\alpha\rho_{yx}+\cos^{2}\alpha\rho_{yy}$
from the resistivity tensors, where $\rho_{ij}$ is the component of the
resistivity tensor and $\alpha$ is the angle between the magnetization
direction and the $x$-axis. The AMR values for different anti-site disordered
Fex(Co)1-xCox(Fe)1-x alloys are shown in FIG. 2 (b), and the solid curves
represent the fitting results. As anti-site disorder increases, the AMR
decreases from 14.1% to 5.3% dramatically, meaning the alloy becomes more
isotropic than before.

The longitudinal and average resistivity $\rho_{xx}$ and $\rho$ can be written
as series expansions with respect to $\alpha$ based on the phenomenological
model [1]. Since we consider the situation when the magnetic moment rotates in
the $xy$-plane, the expansions can be simplified as

$\rho_{xx}=A+B\cos 2\alpha-D\cos 4\alpha,$ (2)

$\rho=\frac{(\rho_{xx}+\rho_{yy}+\rho_{zz})}{3}=(A-\frac{B}{3}-\frac{D}{3}+\frac{\mathrm{E}}{12})-(\frac{2D}{3}+\frac{\mathrm{E}}{12})\cos
4\alpha.$ (3)

$\rho_{xx}$ and $\rho$ are fitted to the magnetic moment angle $\alpha$ and
shown in FIG. 2 (c) and (d). The angular dependence of $\rho_{xx}$ is in good
agreement with the results of Zeng et al. [7]. In addition to the change in
magnetization direction, which affects the resistivity, the anti-site disorder
also significantly affects $\rho_{xx}$ and $\rho$. The comparison of the five
Fex(Co)1-xCox(Fe)1-x alloys in FIG. 2 (c) and (d) can lead to the conclusion
that anti-site disorder causes enhanced scattering and a substantial increase
in resistivity. According to the experiments of Freitas et al., the FeCo
resistivity increases from 2.09 $\mu\Omega\rm{cm}$ to 3.24 $\mu\Omega\rm{cm}$
with the increasing disorder by varying the quenching temperatures [29], which
clearly connects the atomic disorder effect on the electronic transport.

Figure 3: (a) The $\Delta\rho$ and $\Delta\rho_{xx}$ of different disordered
Fex(Co)1-xCox(Fe)1-x alloys. (b) The coefficients of the phenomenological
expression of $\rho_{xx}$ as a function of the concentration of Fe(0,0,0).

To find out the reason for the decrease in AMR with the anti-site disorder, we
define $\Delta\rho=\rho_{\max}-\rho_{\min}$ and
$\Delta\rho_{xx}=\rho_{xx,\max}-\rho_{xx,\min}$. It can be observed in FIG. 3
(a) that $\Delta\rho_{xx}$ fluctuates around 0.22 $\mu\Omega\rm{cm}$, in
agreement with the data of Zeng et al. [7] and Wimmer et al. [30]. The values
of $\Delta\rho_{xx}$ and $\Delta\rho$ are much smaller than $\rho_{xx}$ and
$\rho$, however, indicating that the difference in AMR between different
disordered systems is mainly a result of the resistivity increase caused by
disorder rather than the change in $\rho_{xx}$ and $\rho$ due to the rotation
of the magnetization direction. The enhancement of anti-site disorder also
dramatically raises the value of the constant term $A$. The $\cos 4\alpha$
term’s coefficient $D$, which characterizes the single crystal anisotropic
magnetoresistance, reduces remarkably with the degree of disorder, as shown in
FIG. 3 (b).

The vertex correction is crucial for discussing the scattering of impurities
in electronic transport [31, 32], which corresponds to the scattering-in term
in the Boltzmann equation [33, 34]. However, the effect of the vertex
correction is much weaker if the states at the Fermi energy have $d$
properties primarily [35, 36, 37]. In KKR-CPA calculation, the vertex
correction of conductivity obtained based on the Kubo linear response theory
is the difference between correlated and uncorrelated configurational averages

$\sigma_{\mu\nu}^{\rm{VC}}=\left\langle\hat{O}_{\mu}\hat{G}^{+}\hat{J}_{\nu}\hat{G}^{-}\right\rangle_{c}-\left\langle\hat{O}_{\mu}\hat{G}^{+}\right\rangle_{c}\left\langle\hat{J}_{\nu}\hat{G}^{-}\right\rangle_{c}.$
(4)

FIG. 4 (a) shows the vertex correction of the longitudinal conductivity,
$\sigma_{xx}^{\mathrm{VC}}$, which strongly depends on the direction of
magnetization. From the numerical point of view, however,
$\sigma_{xx}^{\mathrm{VC}}$ accounts for only 5% to 7% of $\sigma_{xx}$ in the
Fex(Co)1-xCox(Fe)1-x alloy since its Fermi level is mainly occupied by the
3$d$ electrons. In a semiclassical picture, when $\mathbf{M}$ gradually
approaches the $x$-axis, the scattering enhancement causes the scattering
events experienced by different electrons to tend to be the same. Therefore,
the difference between correlated and uncorrelated configurational averages,
i.e., the vertex correction $\sigma_{xx}^{\mathrm{VC}}$, decreases. The higher
the degree of anti-site disorder, the stronger the scattering, which leads to
less sensitivity to the rotation of $\mathbf{M}$ and a smaller value of
$\sigma_{xx}^{\mathrm{VC}}$, as shown in FIG. 4 (a) and (b). Moreover, the
variation of $\sigma_{xx}^{\mathrm{VC}}$ due to the rotation of $\mathbf{M}$
is much smaller than that introduced by the anti-site disorder, confirming
that the disorder-induced scattering is more substantial in the
Fex(Co)1-xCox(Fe)1-x alloy.

Figure 4: (a) The vertex correction of the longitudinal conductivity
$\sigma_{xx}^{\mathrm{VC}}$ and $\sigma_{xx}^{\mathrm{VC}}$ / $\sigma_{xx}$ as
a function of $\mathbf{M}$. (b) Dependencies of $\rho$, $\rho_{xx}$,
$\rm{AMR}$, and $\sigma_{xx}^{\mathrm{VC}}$ on the anti-site disorder.

Interestingly, the maximal values of $\rho$ and $\rho_{xx}$ do not occur at
$x=0.5$, but where the concentration of $\mathrm{Fe}_{(0,0,0)}$ is around 45%,
as shown in FIG. 4 (b). Therefore, the differences of $\rho$ and $\rho_{xx}$
between $x=0.5$ and $x=0.4$ in FIG. 2 are negligible. The non-monotonic
increment of the conductivity with the disorder depends on the detailed
electronic structure, which will be seen in the following discussions.

###  3.3 The Bloch spectral function

The AMR comes from the relative variation of longitudinal resistivity
$\rho_{xx}$ and total resistivity $\rho$ when the magnetic moment rotates. The
two values may have various dependencies on the magnetic moment for different
alloys. The electronic reason for the difference can be observed from the
Bloch spectral function (BSF), with the color bar indicating the magnitude of
$\log_{10}A(\vec{k},E)$. The absence of translational symmetry due to random
occupation of the lattice sites makes the Bloch state no longer an eigenstate
of the system. As a result, the system does not have the well-defined
dispersion found in ordered systems. The smeared bands in the BSF caused by
scattering imply a reduction of the electron mean free path [38]. At the same
time, when the relativistic effect is taken into account, the spin is no
longer a good quantum number. The spin-up and down states are mixed when the
magnetization directions vary in space. [6]. We show the total Bloch spectral
function in FIG. 5, where $\mathrm{X}_{1}$ is the high symmetry point X in the
[010] direction.

Figure 5: The total BSF of Fe0.1(Co)0.9Co0.1(Fe)0.9 when $\alpha=0^{\circ}$
(left), $\alpha=45^{\circ}$ (middle), and $90^{\circ}$ (right). The Fermi
level is set at 0 eV.

The disorder-induced smearing of the state at the Fermi level remains under
different magnetization directions. In addition, we find no difference in the
DOS at the $E_{\mathrm{F}}$ for the different magnetization directions.
However, the band-crossing in the $\mathrm{X}-\Gamma$ path about 0.5 eV below
the Fermi level gradually opens. The band-crossing responds to the AMR in the
ordered compounds, as recently proposed by Zeng et al. [7]. Simultaneously, a
change of degeneracy occurs near the high symmetry point $\mathrm{X}_{1}$ 0.4
eV below the $E_{\mathrm{F}}$ in the $\mathrm{X}_{1}-\Gamma$ path, as shown in
FIG. 5. Although they are slightly away from the $E_{\mathrm{F}}$, because of
the smeared band, it can be qualitatively assumed that the alteration of the
magnetization direction affects the electronic transport mainly in the form of
intrinsic mechanisms [7]. There is no difference between the BSFs when
$\alpha=0^{\circ}$ and $90^{\circ}$, except that the band
crossing/anticrossing positions are interchanged due to the symmetry of the
crystal. For $\alpha=45^{\circ}$, however, the crossing opens on both
$\mathrm{X}_{1}-\Gamma$ and $\mathrm{X}-\Gamma$ paths. In other words, the
band-crossing caused by the change of $\mathbf{M}$ has the same fourfold
symmetry as what contributes to the $\cos 4\alpha$ term of AMR.

To examine the specific effect of changes in the band structure on the
electronic transport, we manually alter the positions of $E_{\mathrm{F}}$ of
Fex(Co)1-xCox(Fe)1-x with $x=0.1$ and $x=0.5$, from which the resistivity
contributed from the corresponding positions can be obtained. Thereby, one can
further probe the relationship between the band crossing/anticrossing caused
by the variation of $\mathbf{M}$ and the angular dependence of resistivity. In
FIG. 6, the band-anticrossing of Fe0.1(Co)0.9Co0.1(Fe)0.9 appears around -0.6
eV, and $\rho_{xx}$ is notably enhanced around there. While for
Fe0.5(Co)0.5Co0.5(Fe)0.5, the band-anticrossing occurs around -0.5 eV due to
the lowered $E_{\mathrm{F}}$ caused by the higher anti-site disorder. The
energy corresponding to the peak of $\rho_{xx}$ is accordingly lifted to -0.5
eV. On the other hand, since Fe0.5(Co)0.5Co0.5(Fe)0.5 is much more disordered
than Fe0.1(Co)0.9Co0.1(Fe)0.9, its $\rho_{xx}$ curve with energy is much
smoother and numerically higher than that of Fe0.1(Co)0.9Co0.1(Fe)0.9.

Figure 6: The total BSF of Fe0.1(Co)0.9Co0.1(Fe)0.9 (a) and
Fe0.5(Co)0.5Co0.5(Fe)0.5 (b) at $\alpha=0^{\circ}$. The Fermi level is set at
0 eV in the BSF. $\rho_{xx,\alpha=0^{\circ}}$ is plotted as a function of
energy.

From the comparison of Fex(Co)1-xCox(Fe)1-x with $x=0.1$ and $x=0.5$ in FIG.
6, in addition to the change in the band structure, the smearing is
significantly strengthened in the BSF of Fe0.5(Co)0.5Co0.5(Fe)0.5, indicating
the pronounced scattering enhancement. With the results of vertex correction,
$\sigma_{xx}^{\mathrm{VC}}$ decreases with the augmentation of $S$, confirming
that anti-site disorder drastically enhances the scattering, leading to a
larger resistivity.

In Sec. III. B, we have pointed out that the maximum values of resistivity and
AMR under different anti-site disorders occur at Fe0.45(Co)0.55Co0.45(Fe)0.55
instead of Fe0.5(Co)0.5Co0.5(Fe)0.5. To find out the reason for this
phenomenon, we analyzed the BSFs of the two Fex(Co)1-xCox(Fe)1-x alloys at
different energy in FIG. 7 and found that they overlap in most cases. However,
at -0.20 eV, the BSFs of the two diverge distinctly at the $\mathrm{X}$ and
$\mathrm{M}$ points. It can be inferred that the increment of the resistivity
when $x=0.45$ originates from the states 0.20 eV below the Fermi level.

Figure 7: The BSFs of Fex(Co)1-xCox(Fe)1-x with $x=0.5$ and $0.45$ with the
energy levels at $E_{\mathrm{F}}$ (a) and 0.20 eV below that (b).

##  4 Conclusions

In summary, we have investigated the influences of anti-site disorder in FeCo
alloy on the AMR effect. Our work confirmed that the band-crossing, which
directly determines the variations of the resistivity with respect to the
magnetization direction, gives rise to the AMR. The anti-site disorder
significantly enhances the scattering of the electrons, consequently
increasing the resistivity, while the difference between the resistivity with
magnetic moment parallel and perpendicular to the current only decreases
slightly. As a consequence, the anti-site disorder affects AMR mainly by
reducing conductivity. This explains the fact that suitable annealing
generally enhances the AMR. The result is explained by examining the
variations of the Bloch spectral function with different degrees of disorder
and magnetic moment directions.

Moreover, we fit the calculated angular dependent resistivity to the
phenomenological model. We find that an alloy with a lower anti-site disorder
has $\rho_{xx}$ deviating further from the $\cos 2\alpha$ relationship. The
$\cos 4\alpha$ term rises remarkably with decreasing disorder. According to
the analysis of the Bloch spectral function, the variation of the $\cos
4\alpha$ term can be attributed to the smearing of the states near the
anticrossing points in the Bloch spectral function.

## References

  * McGuire and Potter [1975] T. McGuire and R. Potter, IEEE Transactions on Magnetics 11, 1018 (1975).
  * Popovic  _et al._ [2006] R. S. Popovic, P. M. Drljaca, and P. Kejik, Sensors and Actuators A: Physical 129, 94 (2006), eMSA 2004.
  * Wang  _et al._ [2018] C. Wang, W. Su, Z. Hu, J. Pu, M. Guan, B. Peng, L. Li, W. Ren, Z. Zhou, Z. Jiang, and M. Liu, IEEE Transactions on Magnetics 54, 1 (2018).
  * Fert and Campbell [1968] A. Fert and I. A. Campbell, Phys. Rev. Lett. 21, 1190 (1968).
  * Ebert  _et al._ [2011] H. Ebert, D. Kodderitzsch, and J. Minar, Reports on Progress in Physics 74, 096501 (2011).
  * Ebert  _et al._ [1997] H. Ebert, A. Vernes, and J. Banhart, Solid State Communications 104, 243 (1997).
  * Zeng  _et al._ [2020] F. L. Zeng, Z. Y. Ren, Y. Li, J. Y. Zeng, M. W. Jia, J. Miao, A. Hoffmann, W. Zhang, Y. Z. Wu, and Z. Yuan, Phys Rev Lett 125, 097201 (2020).
  * Dai  _et al._ [2022] Y. Dai, Y. W. Zhao, L. Ma, M. Tang, X. P. Qiu, Y. Liu, Z. Yuan, and S. M. Zhou, Phys. Rev. Lett. 128, 247202 (2022).
  * Yahagi  _et al._ [2020] Y. Yahagi, D. Miura, and A. Sakuma, Journal of the Physical Society of Japan 89, 044714 (2020), https://doi.org/10.7566/JPSJ.89.044714 .
  * Wagenknecht  _et al._ [2020] D. Wagenknecht, K. Výborný, K. Carva, and I. Turek, Journal of Magnetism and Magnetic Materials 513, 167078 (2020).
  * Wagenknecht and Turek [2020] D. Wagenknecht and I. Turek, Journal of Physics: Condensed Matter 32, 185802 (2020).
  * Seeger [2006] A. Seeger, International Journal of Materials Research 97, 861 (2006).
  * Zhang  _et al._ [2017] W. Zhang, T. Wu, B. Peng, and W. Zhang, Journal of Alloys and Compounds 696, 234 (2017).
  * Toth  _et al._ [2010] B. G. Toth, L. Peter, A. Revesz, J. Padar, and I. Bakonyi, European Physical Journal B 75, 167 (2010).
  * Wang and Zhang [2011] Z. Wang and H. Zhang,  _PROCEEDINGS OF THE 4TH INTERNATIONAL CONFERENCE ON MECHANICAL ENGINEERING AND MECHANICS_ , 343 (2011), 4th International Conference on Mechanical Engineering and Mechanics, Suzhou, PEOPLES R CHINA, AUG 11-12, 2011.
  * Bruski  _et al._ [2012] P. Bruski, K.-J. Friedland, R. Farshchi, J. Herfort, and M. Ramsteiner, Solid State Communications 152, 1131 (2012).
  * Soven [1967] P. Soven, Physical Review 156, 809 (1967).
  * Vosko  _et al._ [1980] S. H. Vosko, L. Wilk, and M. Nusair, Canadian Journal of Physics 58, 1200 (1980).
  * Alam  _et al._ [2014] A. Alam, S. N. Khan, A. V. Smirnov, D. M. Nicholson, and D. D. Johnson, Physical Review B 90, 205102 (2014).
  * Mankovsky  _et al._ [2013] S. Mankovsky, D. Kodderitzsch, G. Woltersdorf, and H. Ebert, Physical Review B 87, 014430 (2013).
  * Ayaz Khan  _et al._ [2016] S. Ayaz Khan, P. Blaha, H. Ebert, J. Minar, and O. Sipr, Physical Review B 94, 144436 (2016).
  * Chadova  _et al._ [2016] K. Chadova, D. Kodderitzsch, J. Minar, H. Ebert, J. Kiss, S. W. D’Souza, L. Wollmann, C. Felser, and S. Chadov, Physical Review B 93, 195102 (2016).
  * Streda [1982] P. Streda, Journal of Physics C: Solid State Physics 15, L717 (1982).
  * Kaspar and Salahub [1983] J. Kaspar and D. R. Salahub, Journal of Physics F: Metal Physics 13, 311 (1983).
  * Jacob  _et al._ [2020] G. A. Jacob, S. Sellaiyan, A. Uedono, and R. J. Joseyphus, Applied Physics A 126, 120 (2020).
  * Spooner  _et al._ [1972] S. Spooner, J. W. Lynn, and J. W. Cable, AIP Conference Proceedings 5, 1415 (1972) .
  * Kulkova  _et al._ [2002] S. E. Kulkova, D. V. Valujsky, J. Sam Kim, G. Lee, and Y. M. Koo, Physica B: Condensed Matter 322, 236 (2002).
  * Schwarz  _et al._ [1984] K. Schwarz, P. Mohn, P. Blaha, and J. Kubler, Journal of Physics F: Metal Physics 14, 2659 (1984).
  * Freitas and Berger [1988] P. P. Freitas and L. Berger, Physical Review B 37, 6079 (1988).
  * Wimmer [2018] S. Wimmer,  _Spin caloric transport and related phenomena from first principles_ , Ph.D. thesis (2018).
  * Tulip  _et al._ [2008] P. R. Tulip, J. B. Staunton, S. Lowitzer, D. Kodderitzsch, and H. Ebert, Phys. Rev. B 77, 165116 (2008).
  * Sipr  _et al._ [2020] O. Sipr, S. Wimmer, S. Mankovsky, and H. Ebert, Physical Review B 101, 085109 (2020).
  * Butler [1985] W. H. Butler, Physical Review B 31, 3260 (1985).
  * Turek  _et al._ [2019] I. Turek, J. Kudrnovsky, and V. Drchal, Physical Review B 100, 134435 (2019).
  * Banhart  _et al._ [1994] J. Banhart, H. Ebert, P. Weinberger, and J. Voitländer, Phys. Rev. B 50, 2104 (1994).
  * Turek  _et al._ [2004] I. Turek, J. Kudrnovsky, V. Drchal, and P. Weinberger, Journal of Physics: Condensed Matter 16, S5607 (2004).
  * Ebert  _et al._ [2015] H. Ebert, S. Mankovsky, K. Chadova, S. Polesya, J. Minar, and D. Kodderitzsch, Physical Review B 91, 165132 (2015).
  * Zhang  _et al._ [2015] Y. Zhang, G. M. Stocks, K. Jin, C. Lu, H. Bei, B. C. Sales, L. Wang, L. K. Béland, R. E. Stoller, G. D. Samolyuk, M. Caro, A. Caro, and W. J. Weber, Nature Communications 6, 8736 (2015).


where $\kappa$ denotes the tension of the domain wall connecting the de Sitter
region and the AdS region. This is distinct from the tension associated with
the excitation used to glue two copies of a Euclidean AdS black hole, and
$\kappa$ in Eq. (A.2) is in general different from the mass of the excitation
$\kappa_{0}$ in Eq. (A.1).

Denoting the metric of the Euclidean de Sitter (sphere) by
$ds^{2}=d\theta^{2}+\sin^{2}\\!\theta\,d\varphi^{2}$, the dilaton profile of
the de Sitter side is $\Phi_{\rm dS}=B_{0}+B\cos\theta$.111111 In this
appendix, we include the constant part to the dilaton profiles.  It has two
horizons, the cosmological horizon at $\theta=0$ and the black hole horizon at
$\theta=\pi$. Let $\theta=\theta(\varphi)$ be the brane trajectory in the de
Sitter side. Then, the junction condition (A.2) is again recast into the
equation of motion of particle in one dimension

$\dot{\theta}^{2}+V(\theta)=0$ (A.3)

with

$V(\theta)=\left(\frac{(B\cos\theta+B_{0})^{2}-\kappa^{2}-B^{2}\sin^{2}\\!\theta}{2\kappa
B\sin\theta}\right)^{2}-1,\quad A\ll B,B_{0}.$ (A.4)

This is the two-dimensional analog of the Euclidean potential problem for
higher dimensions, Eq. (3.2). Again, the trajectory oscillates in the bounded
region $\theta_{{\rm min}}\leq\theta\leq\theta_{{\rm max}}$.

It remains to specify the location of the branching $\rho=\rho_{1}$. Let
$\xi_{R}$ be the outward point normal of the right part of the domain wall at
the branching point and $\xi_{L}$ be the similar normal vector for the left
part, as in Fig. 13.

Figure 13:  The branching point $\rho=\rho_{1}$ of the excitation in AdS2
(black line) with the domain wall (two blue lines). We denote the angle
between the blue and the black line by $\alpha$.

Denoting by $\dot{\rho}_{+}$ the velocity of the domain wall in $\rho$
direction at the branching point, its normal vector is

$\xi_{R}=\left(\sqrt{1-\dot{\rho}_{+}^{2}},\;\frac{\dot{\rho}_{+}}{\sinh\rho_{1}}\right).$
(A.5)

The normal vector for the excitation trajectory is

$\xi_{{\rm
e}}=\left(\frac{\sinh\rho_{m}}{\sinh\rho_{1}},\;\sqrt{1-\left(\frac{\sinh\rho_{m}}{\sinh\rho_{1}}\right)}\right),\quad\sinh\rho_{m}=\frac{\kappa_{0}}{2A},$
(A.6)

where $\rho=\rho_{m}$ is the closest approach to the horizon, located at
$\rho=0$. If we denote the angle between these two normal vectors by
$\pi-\alpha$, then the angle between two normal vectors $\xi_{R},\xi_{L}$ in
the de Sitter side is $2\alpha$. This is because there is no conical
singularity in the geometry at the branching point, a conclusion coming from
the equations of motion for the dilaton profile.

This yields the following equation:

$\frac{\sinh\rho_{m}}{\sinh\rho_{1}}\sqrt{1-\dot{\rho}_{+}^{2}}+\left(1-\frac{B\sin\theta_{1}}{B\cos\theta_{1}+B_{0}}\right)\sqrt{1-\left(\frac{\sinh\rho_{m}}{\sinh\rho_{1}}\right)^{2}}=0,$
(A.7)

which is solved to obtain $\theta=\theta_{1}$ (or equivalently $\rho=\rho_{1}$
in the coordinate in the AdS side). When $\rho_{1}\gg\rho_{m}$, one can easily
see that $\theta_{1}\sim\theta_{{\rm min}}$.

## References

  * (1) G. Gibbons and S. Hawking, _Cosmological event horizons, thermodynamics, and particle creation_ , _Phys. Rev. D_ 15 (1977) 2738.
  * (2) Y. Chen, V. Gorbenko and J. Maldacena, _Bra-ket wormholes in gravitationally prepared states_ , _JHEP_ 02 (2021) 009 [2007.16091].
  * (3) T. Hartman, Y. Jiang and E. Shaghoulian, _Islands in cosmology_ , _JHEP_ 11 (2020) 111 [2008.01022].
  * (4) V. Balasubramanian, A. Kar and T. Ugajin, _Islands in de Sitter space_ , _JHEP_ 02 (2021) 072 [2008.05275].
  * (5) W. Sybesma, _Pure de Sitter space and the island moving back in time_ , _Class. Quant. Grav._ 38 (2021) 145012 [2008.07994].
  * (6) H. Geng, Y. Nomura and H.-Y. Sun, _Information paradox and its resolution in de Sitter holography_ , _Phys. Rev. D_ 103 (2021) 126004 [2103.07477].
  * (7) S. Fallows and S. F. Ross, _Islands and mixed states in closed universes_ , _JHEP_ 07 (2021) 022 [2103.14364].
  * (8) L. Aalsma and W. Sybesma, _The price of curiosity: information recovery in de Sitter space_ , _JHEP_ 05 (2021) 291 [2104.00006].
  * (9) K. Langhoff, C. Murdia and Y. Nomura, _Multiverse in an inverted island_ , _Phys. Rev. D_ 104 (2021) 086007 [2106.05271].
  * (10) S. E. Aguilar-Gutierrez, A. Chatwin-Davies, T. Hertog, N. Pinzani-Fokeeva and B. Robinson, _Islands in multiverse models_ , _JHEP_ 11 (2021) 212 [2108.01278].
  * (11) J. Kames-King, E. M. H. Verheijden and E. P. Verlinde, _No Page curves for the de Sitter horizon_ , _JHEP_ 03 (2022) 040 [2108.09318].
  * (12) L. Aalsma, S. E. Aguilar-Gutierrez and W. Sybesma, _An outsider’s perspective on information recovery in de Sitter space_ , _JHEP_ 01 (2023) 129 [2210.12176].
  * (13) J.-H. Baek and K.-S. Choi, _Islands in proliferating de Sitter spaces_ , _JHEP_ 05 (2023) 098 [2212.14753].
  * (14) G. Penington, _Entanglement wedge reconstruction and the information paradox_ , _JHEP_ 09 (2020) 002 [1905.08255].
  * (15) A. Almheiri, N. Engelhardt, D. Marolf and H. Maxfield, _The entropy of bulk quantum fields and the entanglement wedge of an evaporating black hole_ , _JHEP_ 12 (2019) 063 [1905.08762].
  * (16) A. Almheiri, R. Mahajan, J. Maldacena and Y. Zhao, _The Page curve of Hawking radiation from semiclassical geometry_ , _JHEP_ 03 (2020) 149 [1908.10996].
  * (17) D. Marolf and H. Maxfield, _Transcending the ensemble: baby universes, spacetime wormholes, and the order and disorder of black hole information_ , _JHEP_ 08 (2020) 044 [2002.08950].
  * (18) J. McNamara and C. Vafa, _Baby universes, holography, and the swampland_ , 2004.06738.
  * (19) Y. Nomura, N. Salzetta, F. Sanches and S. J. Weinberg, _Toward a holographic theory for general spacetimes_ , _Phys. Rev. D_ 95 (2017) 086002 [1611.02702].
  * (20) Y. Nomura, P. Rath and N. Salzetta, _Spacetime from unentanglement_ , _Phys. Rev. D_ 97 (2018) 106010 [1711.05263].
  * (21) Y. Nomura, _Spacetime and universal soft modes — black holes and beyond_ , _Phys. Rev. D_ 101 (2020) 066024 [1908.05728].
  * (22) C. Murdia, Y. Nomura and K. Ritchie, _Black hole and de Sitter microstructures from a semiclassical perspective_ , _Phys. Rev. D_ 107 (2023) 026016 [2207.01625].
  * (23) L. Susskind, _De Sitter holography: fluctuations, anomalous symmetry, and wormholes_ , _Universe_ 7 (2021) 464 [2106.03964].
  * (24) L. Susskind, _Entanglement and chaos in de Sitter space holography: an SYK example_ , _JHAP_ 1 (2021) 1 [2109.14104].
  * (25) E. Shaghoulian, _The central dogma and cosmological horizons_ , _JHEP_ 01 (2022) 132 [2110.13210].
  * (26) E. Shaghoulian and L. Susskind, _Entanglement in de Sitter space_ , _JHEP_ 08 (2022) 198 [2201.03603].
  * (27) A. Strominger, _The dS / CFT correspondence_ , _JHEP_ 10 (2001) 034 [hep-th/0106113].
  * (28) V. Balasubramanian, J. de Boer and D. Minic, _Mass, entropy and holography in asymptotically de Sitter spaces_ , _Phys. Rev. D_ 65 (2002) 123508 [hep-th/0110108].
  * (29) V. Balasubramanian, J. de Boer and D. Minic, _Notes on de Sitter space and holography_ , _Class. Quant. Grav._ 19 (2002) 5655 [hep-th/0207245].
  * (30) M. Alishahiha, A. Karch, E. Silverstein and D. Tong, _The dS/dS correspondence_ , _AIP Conf. Proc._ 743 (2004) 393 [hep-th/0407125].
  * (31) M. Alishahiha, A. Karch and E. Silverstein, _Hologravity_ , _JHEP_ 06 (2005) 028 [hep-th/0504056].
  * (32) X. Dong, E. Silverstein and G. Torroba, _De Sitter holography and entanglement entropy_ , _JHEP_ 07 (2018) 050 [1804.08623].
  * (33) V. Gorbenko, E. Silverstein and G. Torroba, _dS/dS and $T\overline{T}$_, _JHEP_ 03 (2019) 085 [1811.07965].
  * (34) L. Aalsma and G. Shiu, _Chaos and complementarity in de Sitter space_ , _JHEP_ 05 (2020) 152 [2002.01326].
  * (35) H. Geng, _Non-local entanglement and fast scrambling in de Sitter holography_ , _Annals Phys._ 426 (2021) 168402 [2005.00021].
  * (36) L. Aalsma, A. Cole, E. Morvan, J. P. van der Schaar and G. Shiu, _Shocks and information exchange in de Sitter space_ , _JHEP_ 10 (2021) 104 [2105.12737].
  * (37) A. Svesko, E. Verheijden, E. P. Verlinde and M. R. Visser, _Quasi-local energy and microcanonical entropy in two-dimensional nearly de Sitter gravity_ , _JHEP_ 08 (2022) 075 [2203.00700].
  * (38) V. Franken, H. Partouche, F. Rondeau and N. Toumbas, _Bridging the static patches: de Sitter holography and entanglement_ , 2305.12861.
  * (39) V. Balasubramanian, A. Kar and T. Ugajin, _Entanglement between two gravitating universes_ , _Class. Quant. Grav._ 39 (2022) 174001 [2104.13383].
  * (40) G. Penington, S. H. Shenker, D. Stanford and Z. Yang, _Replica wormholes and the black hole interior_ , _JHEP_ 03 (2022) 205 [1911.11977].
  * (41) A. Almheiri, T. Hartman, J. Maldacena, E. Shaghoulian and A. Tajdini, _Replica wormholes and the entropy of Hawking radiation_ , _JHEP_ 05 (2020) 013 [1911.12333].
  * (42) H. Geng, A. Karch, C. Perez-Pardavila, S. Raju, L. Randall, M. Riojas et al., _Information transfer with a gravitating bath_ , _SciPost Phys._ 10 (2021) 103 [2012.04671].
  * (43) L. Anderson, O. Parrikar and R. M. Soni, _Islands with gravitating baths: towards ER = EPR_ , _JHEP_ 21 (2020) 226 [2103.14746].
  * (44) H. Geng, S. Lüst, R. K. Mishra and D. Wakeham, _Holographic BCFTs and communicating black holes_ , _JHEP_ 08 (2021) 003 [2104.07039].
  * (45) V. Balasubramanian, B. Craps, M. Khramtsov and E. Shaghoulian, _Submerging islands through thermalization_ , _JHEP_ 10 (2021) 048 [2107.14746].
  * (46) N. Iizuka, A. Miyata and T. Ugajin, _A comment on a fine-grained description of evaporating black holes with baby universes_ , _JHEP_ 09 (2022) 158 [2111.07107].
  * (47) A. Miyata and T. Ugajin, _Entanglement between two evaporating black holes_ , _JHEP_ 09 (2022) 009 [2111.11688].
  * (48) H. Geng, A. Karch, C. Perez-Pardavila, S. Raju, L. Randall, M. Riojas et al., _Entanglement phase structure of a holographic BCFT in a black hole background_ , _JHEP_ 05 (2022) 153 [2112.09132].
  * (49) Z. Fu and D. Marolf, _Bag-of-gold spacetimes, Euclidean wormholes, and inflation from domain walls in AdS/CFT_ , _JHEP_ 11 (2019) 040 [1909.02505].
  * (50) M. Mirbabayi, _Uptunneling to de Sitter_ , _JHEP_ 09 (2020) 070 [2003.05460].
  * (51) A. Sahu, P. Simidzija and M. Van Raamsdonk, _Bubbles of cosmology in AdS/CFT_ , 2306.13143.
  * (52) S. E. Aguilar-Gutierrez, E. Bahiru and R. Espíndola, _The centaur-algebra of observables_ , 2307.04233.
  * (53) S. Antonini, M. Sasieta and B. Swingle, _Cosmology from random entanglement_ , 2307.14416.
  * (54) B. Freivogel, V. E. Hubeny, A. Maloney, R. C. Myers, M. Rangamani and S. Shenker, _Inflation in AdS/CFT_ , _JHEP_ 03 (2006) 007 [hep-th/0510046].
  * (55) E. Farhi and A. H. Guth, _An obstacle to creating a universe in the laboratory_ , _Phys. Lett. B_ 183 (1987) 149.
  * (56) E. Farhi, A. H. Guth and J. Guven, _Is it possible to create a universe in the laboratory by quantum tunneling?_ , _Nucl. Phys. B_ 339 (1990) 417.
  * (57) K. Sato, M. Sasaki, H. Kodama and K.-i. Maeda, _Creation of wormholes by first order phase transition of a vacuum in the early universe_ , _Prog. Theor. Phys._ 65 (1981) 1443.
  * (58) K.-i. Maeda, K. Sato, M. Sasaki and H. Kodama, _Creation of de Sitter-Schwarzschild wormholes by a cosmological first order phase transition_ , _Phys. Lett. B_ 108 (1982) 98.
  * (59) K. Sato, H. Kodama, M. Sasaki and K.-i. Maeda, _Multiproduction of universes by first order phase transition of a vacuum_ , _Phys. Lett. B_ 108 (1982) 103.
  * (60) H. Kodama, M. Sasaki and K. Sato, _Abundance of primordial holes produced by cosmological first order phase transition_ , _Prog. Theor. Phys._ 68 (1982) 1979.
  * (61) D. Anninos and D. M. Hofman, _Infrared realization of dS 2 in AdS2_, _Class. Quant. Grav._ 35 (2018) 085003 [1703.04622].
  * (62) E. Witten, _Deformations of JT gravity and phase transitions_ , 2006.03494.
  * (63) D. Anninos and D. A. Galante, _Constructing AdS 2 flow geometries_, _JHEP_ 02 (2021) 045 [2011.01944].
  * (64) S. Chapman, D. A. Galante and E. D. Kramer, _Holographic complexity and de Sitter space_ , _JHEP_ 02 (2022) 198 [2110.05522].
  * (65) F. Ecker, D. Grumiller and R. McNees, _dS 2 as excitation of AdS2_, _SciPost Phys._ 13 (2022) 119 [2204.00045].
  * (66) V. Chandrasekaran, R. Longo, G. Penington and E. Witten, _An algebra of observables for de Sitter space_ , _JHEP_ 02 (2023) 082 [2206.10780].
  * (67) M. Van Raamsdonk, _Building up spacetime with quantum entanglement_ , _Gen. Rel. Grav._ 42 (2010) 2323 [1005.3035].
  * (68) J. Maldacena and L. Susskind, _Cool horizons for entangled black holes_ , _Fortsch. Phys._ 61 (2013) 781 [1306.0533].
  * (69) L. Susskind and J. Uglum, _Black hole entropy in canonical quantum gravity and superstring theory_ , _Phys. Rev. D_ 50 (1994) 2700 [hep-th/9401070].
  * (70) J. H. Cooperman and M. A. Luty, _Renormalization of entanglement entropy and the gravitational effective action_ , _JHEP_ 12 (2014) 045 [1302.1878].
  * (71) A. Goel, H. T. Lam, G. J. Turiaci and H. Verlinde, _Expanding the black hole interior: partially entangled thermal states in SYK_ , _JHEP_ 02 (2019) 156 [1807.03916].
  * (72) V. Balasubramanian, A. Lawrence, J. M. Magan and M. Sasieta, _Microscopic origin of the entropy of black holes in general relativity_ , 2212.02447.
  * (73) J. Chandra and T. Hartman, _Coarse graining pure states in AdS/CFT_ , 2206.03414.
  * (74) V. Balasubramanian, A. Lawrence, J. M. Magan and M. Sasieta, _Microscopic origin of the entropy of astrophysical black holes_ , 2212.08623.
  * (75) V. Balasubramanian, A. Kar and T. Ugajin, _Entanglement between two disjoint universes_ , _JHEP_ 02 (2021) 136 [2008.05274].
  * (76) J. Maldacena, G. J. Turiaci and Z. Yang, _Two dimensional nearly de Sitter gravity_ , _JHEP_ 01 (2021) 139 [1904.01911].
  * (77) J. Cotler, K. Jensen and A. Maloney, _Low-dimensional de Sitter quantum gravity_ , _JHEP_ 06 (2020) 048 [1905.03780].
  * (78) A. Almheiri and J. Polchinski, _Models of AdS 2 backreaction and holography_, _JHEP_ 11 (2015) 014 [1402.6334].
  * (79) J. Maldacena, D. Stanford and Z. Yang, _Conformal symmetry and its breaking in two dimensional nearly anti-de Sitter space_ , _PTEP_ 2016 (2016) 12C104 [1606.01857].
  * (80) D. Bak, C. Kim and S.-H. Yi, _Bulk view of teleportation and traversable wormholes_ , _JHEP_ 08 (2018) 140 [1805.12349].
  * (81) Y. Nakata, T. Takayanagi, Y. Taki, K. Tamaoka and Z. Wei, _New holographic generalization of entanglement entropy_ , _Phys. Rev. D_ 103 (2021) 026005 [2005.13801].
  * (82) M. Miyaji, _Island for gravitationally prepared state and pseudo entanglement wedge_ , _JHEP_ 12 (2021) 013 [2109.03830].
  * (83) A. Lewkowycz and J. Maldacena, _Generalized gravitational entropy_ , _JHEP_ 08 (2013) 090 [1304.4926].
  * (84) P. Calabrese and J. L. Cardy, _Entanglement entropy and quantum field theory_ , _J. Stat. Mech._ 0406 (2004) P06002 [hep-th/0405152].
  * (85) V. Balasubramanian, V. Jejjala and J. Simon, _The library of babel_ , _Int. J. Mod. Phys. D_ 14 (2005) 2181 [hep-th/0505123].
  * (86) A. Miyata and T. Ugajin, _Evaporation of black holes in flat space entangled with an auxiliary universe_ , _PTEP_ 2022 (2022) 013B13 [2104.00183].