aainstitutetext: School of Physics, Peking University, Beijing 100871,
Chinabbinstitutetext: Center for High Energy Physics, Peking University,
Beijing 100871, Chinaccinstitutetext: Department of Physics, Nanjing
University, Nanjing 210093, China

# Distinct photon-ALP propagation modes

Qing-Hong Cao c Zuowei Liu a Jun-Chen Wang<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

Measurement of cosmic photons may reveal their propagation in the interstellar
environment, thereby offering a promising way to probe axions and axion-like
particles (ALPs). Numerical methods are usually used to compute the
propagation of the photon-ALP beam due to the complexity of both the
interstellar magnetic field and the evolution equation. However, under certain
conditions, the evolution equation can be greatly simplified so that the
photon-ALP propagation can be analytically solved. By using analytic methods,
we find two distinct photon-ALP propagation modes, determined by the relative
magnitude of the photon-ALP mixing term in comparison to the photon
attenuation term. In one mode, the intensity of photons decreases with the
increasing distance; in the other mode, it also exhibits oscillatory behavior.
To distinguish the two propagation modes, we compute the observable quantities
such as the photon survival probability and the degree of polarization.

##  1 Introduction

The existence of axions and axion-like particles (ALPs) is a topic of great
interest in modern particle physics Ringwald:2012hr ; Irastorza:2018dyq ;
Adams:2022pbo . The QCD axion was initially proposed as a natural solution to
the strong CP problem Peccei:1977hh ; Wilczek:1977pj . Recently, the study of
ALPs has gained widespread attention due to the much wider range of mass and
coupling parameter than the QCD axion Jaeckel:2010ni . ALPs arise naturally in
a plethora of extensions of the standard model, including supersymmetric
models Chang:1999si ; Turok:1995ai  and superstring theories Svrcek:2006yi ;
Ringwald:2008cu ; Arvanitaki:2009fg .

The presence of ALPs can significantly alter photon propagation in the
universe, leading to distinct photon-ALP oscillation signatures Maiani:1986md
; Raffelt:1987im . Because of the complex astrophysical environment through
which photon-ALP beams propagate, most recent studies used numerical
simulations to analyze their propagation properties; see e.g., Li:2020pcn ;
Horns:2012pp ; Gill:2011yp ; Bi:2020ths ; Galanti:2022yxn ; Galanti:2022tow ;
Xia:2019yud ; Liang:2018mqm . Nevertheless, there are also a number of
analytical studies on the photon-ALP propagation, see e.g., Raffelt:1987im ;
Grossman:2002by ; Csaki:2003ef ; Lai:2006af ; DeAngelis:2007dqd ;
Agarwal:2008ac ; Ganguly:2008kh ; Mirizzi:2009aj ; Wang:2015dil ;
Kartavtsev:2016doq . Although numerical simulations lead to more accurate
results than analytical calculations, some physical phenomena that can be
discerned in analytical calculations are often obscure in numerical
simulations. In this paper we use analytic methods to study the propagation of
the photon-ALP beams. We find that under certain conditions the evolution
equation of the photon-ALP beam can be significantly simplified so that
analytic methods can be used. Our analytic methods reveal two distinct photon
propagation modes, determined by the relative magnitude of the photon-ALP
mixing term in comparison to the photon attenuation term. More precisely, the
two propagation modes are determined by the sign of the quantity $D$ given in
Eq. (14): in the $D>0$ case, photons exhibit clear oscillations during the
propagation; conversely, in the $D<0$ case, the oscillations are absent. We
investigate the astrophysical conditions that lead to the emergence of these
two propagation modes, and further study the implications for relevant
physical observables such as the photon survival probability and degree of
polarization. We study the distinguishing features of the two propagation
modes in detail, which can be used to discriminate between the two modes, as
well as against the standard model background.

The rest of paper is organized as follows. In section 2, we introduce the ALPs
and the propagation equation governing the photon-ALP system. In section 3, we
introduce the two distinct photon propagation modes that are characterized by
the sign of a quantity $D$. We compute the density matrix via analytic methods
and further discuss the effects of ALPs on the density matrix. In sections 4
and 5, we discuss the physical observables such as the photon survival
probability and the photon degree of polarization in the two propagation
modes: We focus on photons with energies above 100 TeV in section 4, and
photons with energies in the range of $10^{-3}-10^{4}$ GeV in section 5. In
section 6, we summarize our findings. In the appendix, we compare our analytic
method in which a uniform magnetic field is assumed, with numerical
calculations in which magnetic field models that aim to describe the
astrophysical conditions are used.

##  2 Propagation equation

Consider the following effective interaction Lagrangian between the axion-like
particle (ALP) $a$ and the photon

$\mathcal{L}_{\rm
int}=-\frac{1}{4}g_{a\gamma}F^{\mu\nu}\tilde{F}_{\mu\nu}\,a=g_{a\gamma}\mathbf{E}\cdot\mathbf{B}\,a,$
(1)

where $F^{\mu\nu}$ is the electromagnetic stress tensor and
$\tilde{F}_{\mu\nu}$ is the dual, $\mathbf{E}$ and $\mathbf{B}$ are the
electric and magnetic fields, respectively, and $g_{a\gamma}$ is the coupling
constant.

The propagation of the photon-ALP beam along the $z$-direction with energy
$\omega$ can be described by a three-component vector
$\psi=(A_{x},A_{y},a)^{T}$, where $A_{x}$ and $A_{y}$ are the electromagnetic
potentials linearly polarized along the $x$ and $y$ axes, respectively.
Without loss of generality, we consider the case where the external magnetic
field is along the $y$ direction and the equation that governs the propagation
of $\psi$ is given by Raffelt:1987im ; DeAngelis:2011id ; Galanti:2022ijh ,

$\left(i\frac{d}{dz}+\omega+\mathcal{M}\right)\psi=0,$ (2)

where $\mathcal{M}$ is given by

$\mathcal{M}=\begin{pmatrix}\Delta_{x}&0&0\\\ 0&\Delta_{y}&\Delta_{a\gamma}\\\
0&\Delta_{a\gamma}&\Delta_{aa}\end{pmatrix}+\frac{1}{2}\begin{pmatrix}i\Gamma&0&0\\\
0&i\Gamma&0\\\ 0&0&0\end{pmatrix}.$ (3)

Here, $\Delta_{x}$ and $\Delta_{y}$ describe the medium effects, $\Gamma$ is
the absorption rate accounting for the attenuation of photons,
$\Delta_{a\gamma}=g_{a\gamma}B/2$ is the photon-ALP mixing term with
$B=|\boldsymbol{B}|$, and $\Delta_{aa}=-m_{a}^{2}/(2\omega)$ with $m_{a}$
being the ALP mass. The absorption rate $\Gamma$ is mainly caused by the
reaction $\gamma\gamma\to e^{+}e^{-}$ between the propagating photons and
ambient photons, such as cosmic microwave background, extragalactic background
light, and so on Vernetto:2016alq ; Lipari:2018gzn . As shown in the black
line of Fig. (1), $\Gamma$ becomes significant as energy $\omega\gtrsim 100$
TeV.

The medium effects are given by

$\Delta_{x,y}=N\Delta_{\rm QED}+\Delta_{\rm pl}+\Delta_{\rm dis},$ (4)

where $N=2$ $(7/2)$ for $\Delta_{x}$ ($\Delta_{y}$), $\Delta_{\rm QED}$
represents the QED birefringence, $\Delta_{\rm pl}$ represents the plasma
effect, and $\Delta_{\rm dis}$ accounts for dispersion effects from photon-
photon scattering on environmental radiation field Bi:2020ths ;
Galanti:2022yxn . The plasma effect $\Delta_{\rm pl}$ is given by

$\Delta_{\rm pl}=-\frac{\omega_{\rm
pl}^{2}}{2\omega}=-\frac{e^{2}n_{e}}{2m_{e}\omega},$ (5)

where $\omega_{\rm pl}$ is the plasma frequency, $e$ is the QED coupling, and
$m_{e}$ and $n_{e}$ is the mass and number density of the electron,
respectively. For the case of $\omega\lesssim 100\ \text{TeV}$, the QED
birefringence and dispersion effects can be obtained from the Euler-Heisenberg
Lagrangian Dittrich:2000zu ; Latorre:1994cv ; Cougo-Pinto:1999agw ;
Tarrach:1983cb ; Adler:1971wn

$\displaystyle\Delta_{\rm QED}$
$\displaystyle=\frac{\alpha\omega}{45\pi}\left(\frac{Be}{m_{e}^{2}}\right)^{2},$
(6) $\displaystyle\Delta_{\rm dis}$ $\displaystyle=\frac{44\alpha^{2}\rho_{\rm
RF}\omega}{135m_{e}^{4}},$ (7)

where $\alpha$ is the fine structure constant, and $\rho_{\rm RF}$ is the
ambient photon energy density. For the high energy cosmic gamma-ray with
$\omega\gtrsim\omega_{2}=100$ TeV, Euler-Heisenberg approximation breaks down,
and one can use the a scaling function $g_{2}(\omega/\omega_{2})$ to take into
account the gamma-gamma scattering due to background photon Dobrynina:2014qba
. Thus, at photon energy $\omega\gtrsim\omega_{2}=100$ TeV, the quantities
$\Delta_{\rm QED}$ and $\Delta_{\rm dis}$ are both modified to $\Delta_{\rm
QED(dis)}g_{2}(\omega/\omega_{2})$, as follows Dobrynina:2014qba :

$\begin{split}g_{2}(u)&=\frac{15}{\pi^{4}}\int_{0}^{\infty}dx\frac{x^{3}}{e^{x}-1}g_{1}\left(ux\frac{30\zeta(3)}{\pi^{4}}\right),\\\
g_{1}(u)&=\frac{3}{8}\int^{+1}_{-1}d\mu(1-\mu)^{2}g_{0}\left(u\frac{1-\mu}{2}\right),\\\
g_{0}(u)&=\frac{45}{44}\int_{1}^{\infty}du^{\prime}\frac{f(u^{\prime})}{u^{\prime
2}-u^{2}},\\\
f(u)&=\frac{4u(u+1)-2}{u^{3}}\ln(\sqrt{u}+\sqrt{u-1})-\frac{2(u+1)\sqrt{u-1}}{u^{5/2}}.\end{split}$
(8)

Note that Eq. (2) is similar to the Schrödinger equation if the coordinate $z$
is replaced by the time $t$ and $(\omega+{\cal M})$ is replaced by
$-\mathcal{H}$, where $\mathcal{H}$ denotes the Hamiltonian. Thus, analogous
to solving the Schrödinger equation, one can also construct the transition
matrix $U(z)=e^{i\mathcal{M}z}$ for the photon-ALP progapation. For a more
general case, one can define the density matrix
$\rho\equiv\psi\psi^{\dagger}$, which satisfies an equation similar to the Von
Neumann equation,

$i\frac{d}{dz}\rho=\rho\mathcal{M}^{\dagger}-\mathcal{M}\rho.$ (9)

The solution can then be obtained via $\rho(z)=U(z)\rho_{0}U(z)^{\dagger}$,
where $\rho_{0}$ is the initial condition. In our analysis, we assume an
unpolarized photon beam such that $\rho_{0}=\text{diag}(1/2,1/2,0)$. The
upper-left $2\times 2$ submatrix of $\rho$, denoted as $\rho_{\gamma}$, can be
parameterized as Kosowsky:1998mb

$\rho_{\gamma}=\begin{pmatrix}\rho_{x}&\rho_{xy}\\\
\rho_{xy}^{*}&\rho_{y}\end{pmatrix}=\frac{1}{2}\begin{pmatrix}I+Q&U-iV\\\
U+iV&I-Q\end{pmatrix},$ (10)

where $I$, $Q$, $U$, and $V$ are the Stokes parameters. The photon degree of
polarization $\Pi_{L}$ is then given by defin

$\Pi_{L}=\frac{\sqrt{Q^{2}+U^{2}}}{I}.$ (11)

The photon survival probability after propagation is given by DeAngelis:2011id

$P_{\gamma\to\gamma}=\rho_{x}+\rho_{y}=I.$ (12)

##  3 Two different propagation modes

Although recent studies focused on numerical methods to solve the photon-ALP
propagation, due to the complex medium effects and the magnetic fields in the
astrophysical environments, there are instances where analytic calculations
are good approximations to the photon propagation. In this section we first
identify the conditions in which analytic methods to the photon propagation
can be used. We then discuss two different propagation modes that are found in
our analysis.

We will use the high-energy photons as the example in this section, though
analytic methods can also be used for low-energy photons. We first consider
the very high energy photons with $\omega\gtrsim 100$ TeV. At such high
energy, the matrix $\cal M$ can be greatly simplified, because the various
$\Delta$ terms can be neglected except the new physics term $\Delta_{a\gamma}$
and the absorption term $\Gamma$, as shown in Fig. (1), where
$g_{a\gamma}=3\times 10^{-12}\ \text{GeV}^{-1}$ and $B=2\ \mu$G are assumed.
Thus we have

$\mathcal{M}\simeq\frac{1}{2}\begin{pmatrix}i\Gamma&0&0\\\
0&i\Gamma&g_{a\gamma}B\\\ 0&g_{a\gamma}B&0\end{pmatrix}.$ (13)

We further assume that the variation of external magnetic field is relatively
small so that it can be approximated by a uniform magnetic field. In this
case, one can then analytically solve the propagation. We note that the
analytical solution can facilitate the calculations and can reveal some
physics pictures of the problem that are difficult to be seen in the numerical
calculations.

Figure 1: Matrix elements of $\mathcal{M}$ as a function of the photon energy
$\omega$. Here we take $g_{a\gamma}=3\times 10^{-12}\ \text{GeV}^{-1}$, $B=2\
\mu$G, $m_{a}=1\times 10^{-19}\ \text{GeV}$, $n_{e}=1.1\times 10^{-2}\
\text{cm}^{-3}$ accounting for the electron number density in the Milky Way
JMYAO  and $\rho_{\rm RF}=2.01\times 10^{-51}\ \text{GeV}^{4}$. accounting for
the energy density of CMB Dobrynina:2014qba . We use Eq. (6) and Eq. (7) to
compute $\Delta_{\rm QED}$ and $\Delta_{\rm dis}$; for $\omega>$100 TeV, we
further multiply a factor of $g_{2}(\omega/\omega_{2})$ where $\omega_{2}=100$
TeV given in Ref. Dobrynina:2014qba . The $\Gamma$ curve is adopted from Ref.
Lipari:2018gzn .

We then find that in the analytic solutions for the simplified $\cal M$, there
exist two different propagation modes of photons, characterized by the sign of
$D$, which is defined as

$D\equiv g_{a\gamma}^{2}B^{2}-\frac{\Gamma^{2}}{4}.$ (14)

We discuss these two propagation modes: $D>0$ and $D<0$ below.

###  3.1 The $D>0$ propagation mode

We first discuss the $D>0$ case. To solve the propagation analytically, we
first compute the transition matrix $U(z)=e^{i\mathcal{M}z}$, where
$\mathcal{M}$ is given by Eq. (13). Thus, for the $D>0$ case, we have

$U_{D>0}(z)=\frac{e^{-\frac{\Gamma}{4}z}}{\cos\varphi}\begin{pmatrix}e^{-\frac{\Gamma}{4}z}\cos{\varphi}&0&0\\\
0&\cos\left(\varphi+\frac{\Delta}{2}z\right)&i\sin\frac{\Delta}{2}z\\\
0&i\sin\frac{\Delta}{2}z&\cos\left(\varphi-\frac{\Delta}{2}z\right)\end{pmatrix},$
(15)

where we have defined $\Delta=\sqrt{D}$ and
$\varphi=\arctan(\Gamma/(2\Delta))$. We then compute $\rho(z)$ at distance $z$
from the source via $\rho(z)=U(z)\rho_{0}U(z)^{\dagger}$, where
$\rho_{0}=\text{diag}(1/2,1/2,0)$ is initial condition for an unpolarized
photon beam at the source. The matrix elements of $\rho_{\gamma}(z)$ in this
case are given by

$\displaystyle\rho_{x}(z)$ $\displaystyle=\frac{1}{2}e^{-\Gamma z},$ (16)
$\displaystyle\rho_{y}(z)$
$\displaystyle=\frac{1}{2}\frac{e^{-\frac{\Gamma}{2}z}}{\cos^{2}\varphi}\cos^{2}\left(\varphi+\frac{\Delta}{2}z\right),$
(17) $\displaystyle\rho_{xy}(z)$ $\displaystyle=0.$ (18)

Here $\rho_{x}(z)$ and $\rho_{y}(z)$ describe the intensities of photons
polarized along $x$ and $y$ directions, respectively. Because the off-diagonal
element $\rho_{xy}(z)$ vanishes in this case, the photon degree of
polarization $\Pi_{L}$ becomes

$\Pi_{L}=\frac{|\rho_{x}-\rho_{y}|}{\rho_{x}+\rho_{y}}.$ (19)

We note that the intensities of photons with polarization along the $x$ and
$y$ directions exhibit different dependencies on the propagation distance $z$.
We discuss the distinct behaviors below.

First, the intensity of photons polarized along the $y$ direction depends both
on the attenuation term $\Gamma$ and on the ALP-interaction term
$\Delta_{a\gamma}$ (through $\Delta$ and $\varphi$), but the intensity along
the $x$ direction depends on $\Gamma$ only. This discrepancy arises from the
fact that the external magnetic field is taken to be along the $y$ direction
so that photons polarized along the $x$ direction are not directly affected by
the ALP.

Second, the intensity of photons polarized along the $x$ direction decreases
exponentially as the distance $z$ increases, with a decay length of
$L_{d}=\Gamma^{-1}$, as shown in Eq. (16). On the other hand, the intensity
for photons polarized along the $y$ direction appears to exhibit a longer
decay length of $L_{d}=(\Gamma/2)^{-1}$ (twice of that along the $x$
direction), based solely on the argument of the exponential function in Eq.
(17). This change on the decay length is due to the fact that photons
polarized along the $y$ direction can convert into ALPs that do not experience
the photon attenuation effects $\Gamma$. We emphasize that the actual decay
length for photons polarized along the $y$ direction will deviate somewhat
from the value of $L_{d}=(\Gamma/2)^{-1}$ due to the additional $z$-dependence
in Eq. (17). Specifically, taking the limit $g_{a\gamma}\to 0$ in Eq. (17)
introduces an additional exponential factor $e^{-(\Gamma/2)z}$, which must be
taken into account. 111Note that the $g_{a\gamma}\to 0$ limit is not allowed
in the $D>0$ case. See more discussions in the $D<0$ case.

Third, the intensity of photons polarized along the $y$ direction, while
propagating, exhibits an oscillatory behavior. We define the oscillation
length, denoted by $L_{o}$, as the period of the absolute value of the cosine
function in Eq. (17). The oscillation length is given by

$L_{o}=\frac{2\pi}{\Delta}=\frac{2\pi}{\sqrt{g_{a\gamma}^{2}B^{2}-\left(\frac{\Gamma}{2}\right)^{2}}}.$
(20)

Thus, the increase in the difference between $g_{a\gamma}B$ and $\Gamma$ leads
to an increase in the oscillation length. We note that the oscillation is
attenuated by the exponential factor $e^{-(\Gamma/2)z}$.

###  3.2 The $D<0$ propagation mode

We next discuss the $D<0$ case. For the case of $D<0$, the transition matrix
$U(z)$ can be written as

$U_{D<0}(z)=\frac{e^{-\frac{\Gamma}{4}z}}{\cos
2\alpha}\begin{pmatrix}e^{-\frac{\Gamma}{4}z}\cos 2\alpha&0&0\\\
0&e^{-\frac{\Delta}{2}z}\cos^{2}\alpha-e^{\frac{\Delta}{2}z}\sin^{2}\alpha&i\left(e^{\frac{\Delta}{2}z}-e^{-\frac{\Delta}{2}z}\right)\sin\alpha\cos\alpha\\\
0&i\left(e^{\frac{\Delta}{2}z}-e^{-\frac{\Delta}{2}z}\right)\sin\alpha\cos\alpha&e^{\frac{\Delta}{2}z}\cos^{2}\alpha-e^{-\frac{\Delta}{2}z}\sin^{2}\alpha\end{pmatrix},$
(21)

where $\Delta=\sqrt{-D}$ and
$\alpha=\arcsin{\sqrt{\Gamma/2-\Delta}}/{\sqrt{\Gamma}}$. 222The definition of
$\Delta$ in the $D<0$ case is different from the $D>0$ case.

For an unpolarized photon beam at the source, the diagonal matrix elements of
$\rho_{\gamma}(z)$ are given by

$\displaystyle\rho_{x}(z)$ $\displaystyle=\frac{1}{2}e^{-\Gamma z},$ (22)
$\displaystyle\rho_{y}(z)$
$\displaystyle=\frac{1}{2}\frac{e^{-\frac{\Gamma}{2}z}}{\cos^{2}(2\alpha)}e^{-{\Delta}z}\left|\cos^{2}\alpha-e^{\Delta
z}\sin^{2}\alpha\right|^{2},$ (23) $\displaystyle\rho_{xy}(z)$
$\displaystyle=0.$ (24)

We note that both $\rho_{x}(z)$ and $\rho_{xy}(z)$ are the same as in the
$D>0$ case. Similarly to the $D>0$ case, the vanishing off-diagonal element
$\rho_{xy}(z)$ leads to a simplified expression of the photon as given in Eq.
(19). We next discuss the distinct dependencies on the propagation distance
$z$ of the intensities of photons with polarization along the $x$ and $y$
directions, and also make comparison to the $D>0$ case.

First, photons polarized along the $x$ direction only undergo attenuation
characterized by the photon attenuation coefficient $\Gamma$, whereas photons
polarized along the $y$ direction are influenced both by $\Gamma$ and by the
mixing term with the ALP. This is similar to the $D>0$ case.

Second, the ALP-photon mixing term weakens the photon attenuation. This can be
seen from the argument of the exponential function in Eq. (23), which hints a
decay length of $L_{d}=(\Gamma/2+\Delta)^{-1}$ for $y$-polarized photons. The
decay length of $x$-polarized photons is $L_{d}=\Gamma^{-1}$. Since
$\Delta<\Gamma/2$, photons polarized along the $y$-direction has a larger
decay length than photons polarized along the $x$-direction. The reason behind
this is the same as the $D>0$ case: when propagating, $y$-polarized photons
can convert into ALPs which are unaffected by photons attenuation effects.

Third, in contrast to the $D>0$ scenario, where photons polarized in the $y$
direction exhibit oscillations, photons in the $D<0$ scenario do not exhibit
any oscillatory behavior.

##  4 Propagation modes for high-energy photons with $\omega\geq 100$ TeV

In this section, we discuss the propagation modes for photons with energy
$\omega\geq 100$ TeV. We will adopt the value of the coupling constant
$g_{a\gamma}=3\times 10^{-12}$ GeV-1 Irastorza:2021tdu  as the default value
for all subsequent analyses, unless otherwise stated.

In the Milky Way (MW) galaxy, the strength of the magnetic field is of $O(1)\
\mu$G Jansson:2012pc ; Jansson:2012rt . Taking $B=1\ \mu$G , we find that
$Bg_{a\gamma}\simeq 0.01$ kpc-1; equating $Bg_{a\gamma}$ with $\Gamma/2$ leads
to the photon energy at $\sim 250$ TeV. Thus, the photon propagation mode is
the $D<0$ mode for $\omega\gtrsim 250$ TeV, and the $D>0$ mode for
$\omega\lesssim 250$ TeV.

For active galactic nuclei, galaxy clusters, and extragalactic space, the
typical magnetic field strengths are $O(1-10^{5})\ \mu$G Pudritz:2012xj ;
Tavecchio:2009zb , $O(1-10)\ \mu$G Govoni:2004as , and $O(10^{-7}-1)$ nG
Neronov:2010gir ; Pshirkov:2015tua , respectively. Because the magnetic field
in active galactic nuclei and galaxy clusters is typically much larger than
that in the MW galaxy, the propagation mode is expected to be $D>0$ for
photons with energy $\omega\geq 100$ TeV. 333However, if the variation of
external magnetic field is large, the mean value of $B$ may be very small and
the propagation mode would be more similar to the $D<0$ mode. See Fig. (8) for
details. On the other hand, the magnetic field in extragalactic space is much
smaller than that in the MW galaxy, leading to the $D<0$ mode for photons with
energy $\omega\geq 100$ TeV.

We note that our general discussions here are not applicable in extreme
environments where the magnitude of the magnetic field and/or electron number
density is significantly large such that the terms $\Delta_{\rm QED}$ and/or
$\Delta_{\rm pl}$ in the $\mathcal{M}$ matrix can become substantial.

###  4.1 The $D>0$ case

We first discuss the physics in the $D>0$ case. For the $D>0$ mode, we first
compute the photon degree of polarization $\Pi_{L}$, which is given by Eq.
(19), in the absence of the off-diagonal elements of the $\rho_{\gamma}$
matrix. Due to the interaction with the ALP, photons observed at Earth can be
fully polarized, i.e., $\Pi_{L}=1$, in spite of the unpolarized initial
condition at the source position. This is a remarkable signature for ALP
detection. Thus, it is of great interest to determine under what conditions
photons observed can be fully polarized. The nearly fully-polarized
$\Pi_{L}\sim 1$ can be achieved in two cases: (1) $\rho_{x}\gg\rho_{y}$, (2)
$\rho_{x}\ll\rho_{y}$.

We first discuss the $\rho_{x}\gg\rho_{y}$ case. Since $\rho_{x}(z)>0$, the
full-polarization case $\Pi_{L}=1$ can be achieved if $\rho_{y}(z)=0$; by
using the analytic expression of $\rho_{y}(z)$ given in Eq. (17), we obtain
the $z$ values for $\Pi_{L}=1$:

$z_{n}=\frac{2}{\Delta}\left[\frac{n-1}{2}\pi+\arctan\left(\frac{2\Delta}{\Gamma}\right)\right],$
(25)

where $n$ is a positive odd integer. Thus, the full-polarization $\Pi_{L}=1$
phenomena are evenly distributed in space, and the distance between adjacent
$\Pi_{L}$ points is $2\pi/\Delta$. Because $\rho_{x}(z)>0$ decreases with $z$,
the detection should be performed with the smallest distance $z_{1}$ in Eq.
(25), if possible. Note that $z_{1}=2\arctan(2\Delta/\Gamma)/\Delta$ can be
significantly small if $g_{a\gamma}B$ is sufficiently large.

Figure 2: The photon polarization degree $\Pi_{L}$ (left) and the photon
survival probability $P_{\gamma\to\gamma}$ (right) as a function of the
propagation distance $z$ (blue-solid), in the $D>0$ case. Here $B=10\
\mu\text{G}$, $g_{a\gamma}=3\times 10^{-12}\ \text{GeV}^{-1}$, and
$\Gamma=0.02\ \text{kpc}^{-1}$. The approximated expression of $\Pi_{L}$ given
in Eq. (27) is also shown (red-dotted).

Fig. (2) shows the photon polarization degree $\Pi_{L}$ as a function of the
propagation distance $z$, where a constant external magnetic field $B=10\
\mu\text{G}$, the coupling constant $g_{a\gamma}=3\times 10^{-12}\
\text{GeV}^{-1}$, and the absorption rate $\Gamma=0.02\ \text{kpc}^{-1}$
(corresponding to $\omega\sim 250$ TeV) are used. We note that in Fig. (2),
$z_{1}\simeq 31.8$ kpc, and the distance between adjacent peaks in the
$\rho_{x}\gg\rho_{y}$ case is $2\pi/\Delta\simeq 68.3$ kpc.

We next discuss the $\rho_{x}\ll\rho_{y}$ case. Because the $y$-polarized
photons have a longer decay length than the $x$-polarized photons, eventually
$\rho_{x}$ can become much smaller than $\rho_{y}$, leading to a nearly full-
polarization, $\Pi_{L}\simeq 1$. By using the analytic expressions given in
Eq. (16) and Eq. (17), we obtain the maximum values of $\Pi_{L}$ in the
$\rho_{x}\ll\rho_{y}$ case

$\Pi_{L}^{n}=\tanh\left(\frac{\Gamma}{4}z_{n}\right),\,\,\text{with}\,\,z_{n}=\frac{\pi
n}{\Delta},$ (26)

where $n$ is a positive even integer. Unlike the $\rho_{x}\gg\rho_{y}$ case,
where $\Pi_{L}=1$ can occur at relatively small distance, in the
$\rho_{x}\ll\rho_{y}$ case, large $\Pi_{L}$ only occurs at large distance; for
example, in Fig. (2), the maximum value of $\Pi_{L}$ in the
$\rho_{x}\ll\rho_{y}$ case is rather small ($\Pi_{L}<0.77$) until it reaches
the third maximum point in this case, which is $z_{6}=6\pi/\Delta\simeq 205$
kpc.

We further extrapolate the maximum value of $\Pi_{L}$ on $z_{n}$, as given in
Eq. (26), to all $z$ values,

$\Pi_{L}=\tanh\left(\frac{\Gamma}{4}z\right),$ (27)

which can be interpreted as the theoretical upper bound of $\Pi_{L}$, as shown
in Fig. (2). It is interesting to note that Eq. (27) is independent of the
external magnetic field.

We further discuss the photon survival probability. Fig. (2) shows the photon
survival probability $P_{\gamma\to\gamma}$ as a function of the propagation
distance $z$, for the same parameters as in Fig. (2). Because of the
attenuation term $\Gamma$, the photon survival probability decreases with the
propagation distance. Moreover, the presence of the ALP results in an
oscillatory behavior in the photon survival probability. The photon survival
probability $P_{\gamma\to\gamma}$ is bounded from above by
$P_{\gamma\to\gamma}<(e^{-\Gamma z}+e^{-\Gamma z/2}/\cos\phi)/2$, which is
larger than the case without ALP, $P_{\gamma\to\gamma}=e^{-\Gamma z}$. Thus,
the presence of ALP makes distant galaxies brighter, resulting in a better
detection probability.

As shown in Fig. (2), $\Pi_{L}\sim 1$ can occur at large $z$ values both in
the $\rho_{x}\gg\rho_{y}$ case (odd-numbered peaks) and in the
$\rho_{x}\ll\rho_{y}$ case (even-numbered peaks). However, as shown in Fig.
(2), the photon survival probability at odd-numbered peaks can be several-
order-of-magnitude smaller than the adjacent even-numbered peaks at large $z$
values. Thus, photons emitted from a distant source are likely to be observed
in the $\rho_{x}\ll\rho_{y}$ case than the $\rho_{x}\gg\rho_{y}$ case.

###  4.2 The $D<0$ case

Figure 3: The photon polarization degree $\Pi_{L}$ (left) and the photon
survival probability $P_{\gamma\to\gamma}$ (right) as a function of the photon
propagating distance $z$, in the $D<0$ case. A constant external magnetic
$B=2\ \mu\text{G}$, the coupling constant $g_{a\gamma}=3\times 10^{-12}\
\text{GeV}^{-1}$, and the absorption rate $\Gamma=0.1\ \text{kpc}^{-1}$
(corresponding to $\omega\sim 850$ TeV) are used.

We next discuss the physics in the $D<0$ case. The left panel figure of Fig.
(3) shows the degree of polarization $\Pi_{L}$, in the $D<0$ case, as a
function of the propagation distance $z$, where a constant external magnetic
field $B=2\ \mu\text{G}$, $g_{a\gamma}=3\times 10^{-12}\ \text{GeV}^{-1}$, and
$\Gamma=0.1\ \text{kpc}^{-1}$ are used. Similarly to the $D>0$ propagation
mode, the nearly fully-polarized $\Pi_{L}\sim 1$ in the $D<0$ propagation mode
can be also achieved in two cases: (1) $\rho_{x}\gg\rho_{y}$, (2)
$\rho_{x}\ll\rho_{y}$.

We first discuss the $\rho_{x}\gg\rho_{y}$ case. In this case, the full-
polarization cases $\Pi_{L}=1$ occurs when $\rho_{y}(z)=0$, since
$\rho_{x}(z)>0$. Unlike the infinite $z$ values for $\rho_{y}(z)=0$ in the
$D>0$ case, there is only a single point for $\rho_{y}(z)=0$ in the $D<0$
case. By using the analytic expression of $\rho_{y}(z)$ in Eq. (23), we obtain
the single $z$ value for $\rho_{y}(z)=0$ as

$z_{0}=\frac{1}{\Delta}\ln{\frac{\Gamma+2\Delta}{\Gamma-2\Delta}=\frac{2}{\sqrt{\left(\frac{\Gamma}{2}\right)^{2}-g_{a\gamma}^{2}B^{2}}}\ln{\frac{\frac{\Gamma}{2}+\sqrt{\left(\frac{\Gamma}{2}\right)^{2}-g_{a\gamma}^{2}B^{2}}}{g_{a\gamma}B}}}.$
(28)

In Fig. (3), one has that $z_{0}\simeq 71$ kpc, which is the first point where
the full-polarization cases $\Pi_{L}$ reaches one.

We next discuss the $\rho_{x}\ll\rho_{y}$ case. Because the $y$-polarized
photons have a longer decay length than the $x$-polarized photons, eventually
$\rho_{x}$ can become much smaller than $\rho_{y}$, leading to a nearly full-
polarization, $\Pi_{L}\simeq 1$. This occurs at a relatively large $z$ value.
Note that $\rho_{y}$ has a local maximal point at $z=2z_{0}$. Thus, we take
the $z\gtrsim 2z_{0}$ as the region for the nearly full-polarization,
$\Pi_{L}\simeq 1$. In Fig. (3), the nearly full-polarization, $\Pi_{L}\simeq
1$ occurs at $z\gtrsim 2z_{0}\simeq 140$ kpc.

The right panel figure of Fig. (3) shows the photon survival probability
$P_{\gamma\to\gamma}$, in the $D<0$ case. Unlike the $D>0$ case, here
$P_{\gamma\to\gamma}$ does not oscillate with the propagation distance $z$.
Nevertheless, the presence of ALPs can still facilitate the detection of
remote photons more effectively than in the absence of ALPs.

We next investigate the optimal conditions for observing polarization signals
in the $D<0$ propagation mode. In order to observe a distinct polarization
signature with a strong signal strength, two conditions must be met. Firstly,
the degree of polarization should be close to one, namely $\Pi_{L}\approx 1$.
This condition can be achieved either at $z\simeq z_{0}$ or at $z\geq 2z_{0}$.
Secondly, the total intensity of photons should be significant, which implies
a significant photon survival probability
$P_{\gamma\to\gamma}=\rho_{x}+\rho_{y}$. Since $\rho_{x}(z)$ decreases with
increasing $z$, observation at $z=z_{0}$ should be perused as the first
option.

We next study the $z\geq z_{0}$ region. In this region, while the intensity
$\rho_{x}(z)$ continues to decrease with increasing $z$, the intensity
$\rho_{y}(z)$ starts to grow from zero at $z=z_{0}$ to its local maximum point
at $z=2z_{0}$, and then decreases with increasing $z$. Therefore, the maximum
value of the photon survival probability in the $z\geq z_{0}$ region should
occur in the interval of $z_{0}\leq z\leq 2z_{0}$; within this interval, the
maximum value of $\rho_{x}$ is $e^{-\Gamma z_{0}}/2$ at $z=z_{0}$, which is
the same as the maximum value of $\rho_{y}$ at $z=2z_{0}$. The theoretical
upper bound on the photon survival probability can be obtained by simply
summing these two maximum values, leading to $P_{\gamma\to\gamma}<e^{-\Gamma
z_{0}}$. 444The upper bound on the photon survival probability obtained in
this way is larger than or equal to its maximum value. We next study the
possible maximum value of the above upper bound. By using Eq. (28), we have
$\Gamma z_{0}=2x\ln[(x+1)/(x-1)],$ where $x\equiv{\Gamma}/(2\Delta)\geq 1$.
When $x\to\infty$, the quantity $\Gamma z_{0}$ reaches its minimum value,
which is 4. This can be achieved by setting $\Delta=0$, which arises if
$\Gamma/2=g_{a\gamma}B$. This then leads to a maximum value for the upper
bound on the photon survival probability $P_{\gamma\to\gamma}<e^{-4}\simeq
1.8\%$.

Therefore, for $D<0$ propagation mode, the observations on the degree of
polarization $\Pi_{L}$ should be first conducted at the distance of $z\simeq
z_{0}$ and $z\simeq 2z_{0}$. The theoretical upper bound on the photon
survival probability in the $z\geq z_{0}$ region is $\simeq 1.8\%$. We note
that this upper bound can be somewhat larger than the actual maximum value of
the photon survival probability in this region; for example, the photon
survival probability is $\lesssim 0.04\%$ in the $z\geq z_{0}$ region for the
model parameters on Fig. (3).

Figure 4: The photon survival probability $P_{\gamma\to\gamma}$ (upper) and
the photon degree of polarization $\Pi_{L}$ (lower) as a function of the
photon energy. Blue lines are results adopted from Figure 11 of Ref.
Galanti:2022yxn , where photons originate from clusters at $z=0.03$ (left
column) and at $z=0.4$ (right column). The red lines indicate
$P_{\gamma\to\gamma}=1.8\%$.

In Fig. (4), we further compare the theoretical upper bound on the photon
survival probability in the $z\geq z_{0}$ region, which is
$P_{\gamma\to\gamma}<1.8\%$, to the actual calculations where different medium
effects and varying magnetic fields are taken into consideration. We find that
in the energy range where the photon polarization $\Pi_{L}\approx 1$, which is
$\omega\gtrsim 20$ TeV ($\omega\gtrsim 1$ TeV) on the left (right) panel
figure of Fig. (4), the photon survival probability is indeed below $1.8\%$.

##  5 Propagation modes for photons with energy $\omega=10^{-3}-10^{4}$ GeV

In this section, we discuss the propagation modes for photons with with energy
$\omega=10^{-3}-10^{4}$ GeV. In this energy range, as compared to the ALP-
photon mixing term, the quantity $\Gamma$ can also be considered negligible,
in addition to the $\Delta_{x}$ and $\Delta_{y}$ terms, which are considered
to be negligible in section 3. Thus our analytical analysis carried out in
section 3 is still valid and can be further simplified by taking the
$\Gamma\to 0$ limit. In section 3 we have neglected the ALP mass term. In this
section we also study the effects of the ALP mass term.

We first discuss the case where the ALP mass term is neglected. In this case,
we take the $\Gamma\to 0$ limit in Eqs. (16) and (17), which lead to

$\displaystyle\rho_{x}(z)=\frac{1}{2},$ (29)
$\displaystyle\rho_{y}(z)=\frac{1}{2}\cos^{2}\left(\frac{g_{a\gamma}Bz}{2}\right).$
(30)

Thus, both $\rho_{x}(z)$ and $\rho_{y}(z)$ are independent on energy, which
also result in energy-independence of the photon survival probability
$P_{\gamma\to\gamma}$ and the photon degree of polarization $\Pi_{L}$.

Figure 5: Photon survival probability as a function of the energy. Here we
take the magnetic field as $B=5\ \mu\text{G}$, the coupling constant as
$g_{a\gamma}=3\times 10^{-12}$ GeV, the propagation distance as $z=1\
\text{Mpc}$, and the ALP mass as $m_{a}=10^{-10}\ \text{eV}$.

We next discuss the case where the ALP mass term is important. We consider the
case where $m_{a}=10^{-10}\ \text{eV}$. In this case, we have
$\rho_{x}(z)=1/2$. The analytical form of $\rho_{y}(z)$ can be obtained by
replacing $\Gamma$ in Eq. (17) with $2i\Delta_{aa}=-im_{a}^{2}/\omega$. Thus,
we have

$\rho_{y}(z)=\frac{1}{2}\left|\cos\left(\frac{\Delta}{2}z\right)+i\frac{\Delta_{aa}}{\Delta}\sin\left(\frac{\Delta}{2}z\right)\right|^{2},$
(31)

where $\Delta=\sqrt{g_{a\gamma}^{2}B^{2}+m_{a}^{4}/(4\omega^{2})}$. We further
compute the photon survival probability $P_{\gamma\to\gamma}$ and the
polarization degree $\Pi_{L}$ by $P_{\gamma\to\gamma}=1/2+\rho_{y}$ and

$\Pi_{L}=\frac{1-2\rho_{y}}{1+2\rho_{y}}=(P_{\gamma\to\gamma})^{-1}-1.$ (32)

Thus, the polarization degree is solely determined by the photon survival
probability $P_{\gamma\to\gamma}$.

Fig. (5) shows the relation between the photon survival probability
$P_{\gamma\to\gamma}$ and the photon energy $\omega$. Interestingly, we find
that there exist two energy regions where the photon survival probability
$P_{\gamma\to\gamma}$ is (nearly) independent of the photon energy: (1)
$\omega\lesssim 1$ MeV, (2) $\omega\gtrsim 300$ MeV. This can be understood as
follows. If the photon energy $\omega$ is high such that the term
$\Delta_{aa}=-m_{a}^{2}/(2\omega)$ can be neglected, this is just the case
where the ALP mass can be neglected. Thus one should have energy-independent
$\rho_{x}$ and $\rho_{y}$, as given in Eqs. (29) and (30), resulting in
energy-independence in $P_{\gamma\to\gamma}$ for $\omega\gtrsim 300$ MeV in
Fig. (5). If the photon energy is low such that $\Delta_{aa}/\Delta\to 1$, we
have $\rho_{y}\to 1/2$. This then leads to $P_{\gamma\to\gamma}\approx 1$ for
$\omega\lesssim 1$ MeV in Fig. (5).

The measurement of photon polarization degree can be carried out in the energy
range of $10^{-1}-10^{4}\ \text{MeV}$ through advanced gamma-ray detectors,
such as COSI COSI , e-ASTROGAM Kawata:2017zyo ; e-ASTROGAM:2018jlu , and AMEGO
Kierans:2020otl , in the upcoming future. We note that the polarization
degree, $\Pi_{L}$, can be indirectly determined using the relation
$\Pi_{L}=(1/P_{\gamma\to\gamma}-1)$. Additionally, due to the photon-ALP
interaction, the photon energy spectrum can exhibit a distinct oscillatory
pattern sandwiched between two almost energy-independent regions, as
illustrated in Fig. (5), which may serve as a novel signature for the ALP
detection.

##  6 Conclusions

In this paper we have identified two distinct photon propagation modes in the
presence of ALPs. We classify the two different modes by the sign of
$D=(g_{a\gamma}B)^{2}-\Gamma^{2}/4$. For the $D>0$ propagation mode, the
intensity of photon oscillates as propagating, producing multiple peaks along
its propagation path. For the $D<0$ case, on the contrary, the intensity of
photon does not exhibit any oscillatory behavior.

We use analytic methods to study the two photon propagation modes, because
they are not readily discernible in numerical simulations, which have been
extensively used in the literature to model the photon-ALP propagation. In our
analytic methods, we assume a uniform magnetic field and negligible medium
effects so that the propagation equation of the photon-ALP system can be
solved in a simple analytic form.

We investigate photon propagation in two energy regions where our analytic
methods are appropriate: (1) photon with energy $\gtrsim 100$ TeV; (2) photon
in the energy range of $10^{-3}-10^{4}$ GeV. We identify the two photon
propagation modes by comparing the magnetic field and the photon attenuation
rate $\Gamma$ in different astrophysical environments. For the two propagation
modes, We compute the photon survival probability $P_{\gamma\to\gamma}$ and
the degree of photon polarization $\Pi_{L}$.

In the $D>0$ propagation mode of photons with energy $\gtrsim 100$ TeV, the
fully-polarization can occur either in the $\rho_{x}\ll\rho_{y}$ case or in
the $\rho_{y}\ll\rho_{x}$ case. Because of the oscillatory behavior in the
intensity, the fully-polarization exhibits as various peaks in the propagation
distance. In the $D<0$ propagation mode of photons with energy $\gtrsim 100$
TeV, there is no oscillation in the photon intensity. The detection condition
that yields optimal results should include both a a nearly full-polarization
signal and a considerable photon survival probability. The distances at which
this condition is met in the $D<0$ case are $z\simeq z_{0}$ and $z\simeq
2z_{0}$, where $z_{0}$ is defined in Eq. (28). We further find an upper bound
on the photon survival probability, $\simeq 1.8$ %, for the full-polarization
region, $\Pi_{L}\approx 1$.

In the energy interval of $10^{-3}<\omega<10^{4}$ GeV, both medium and
attenuation effects can become small compared to the ALP-photon mixing term,
leading to even simpler analytic results. In this energy range, the
propagation mode is predominately the $D>0$ mode in most of the parameter
space of interest. We further find some distinguishing signatures associated
with the ALP mass $m_{a}$: If $m_{a}\simeq 0$, both $P_{\gamma\to\gamma}$ and
$\Pi_{L}$ are energy-independent; If $m_{a}$ cannot be neglected, there exist
two energy-independent regions separated by an oscilating pattern in the
energy spectrum of $P_{\gamma\to\gamma}$ and $\Pi_{L}$, which may serve as a
novel signature to detect ALPs in the future experiments.

##  7 Acknowledgements

The work is supported in part by the National Natural Science Foundation of
China under Grant Nos. 11675002, 11725520, 12147103, 12235001, and 12275128.
ZL thank Yonglin Li and Zi-Wei Tang for discussions.

##  Appendix A Comparison between analytic and numerical calculations

In this section we compare analytic calculations based on a uniform magnetic
field assumption, with the numerical calculations using a more realistic
consideration of the magnetic field, both in the MW galaxy and in galaxy
clusters.

###  A.1 MW galaxy

We first discuss the photon propagation in the MW galaxy. We adopt the MW
magnetic field model given in Refs. Jansson:2012pc ; Jansson:2012rt , in which
the magnetic field consists of two components: the regular component and the
random component. Following Ref. Bi:2020ths , we neglect the random fields in
our analysis. For completeness, we present the regular component of the MW
magnetic field from Refs. Jansson:2012pc ; Jansson:2012rt  in Table 1.

$i$ | 0 | 1 | 2 | 3 | 4 | 5 | 6 | 7
---|---|---|---|---|---|---|---|---
$r_{xi}$ (kpc) | 5.1 | 6.3 | 7.1 | 8.3 | 9.8 | 11.4 | 12.7 | 15.5
$b_{i}$ ($\mu$G) |  | 0.1 | 3.0 | -0.9 | -0.8 | -2.0 | -4.2 | 0
Table 1: The regular component of the magnetic field model in the MW disk
Jansson:2012pc . The disk is partitioned into 7 regions with 8 boundaries that
are given by $r_{i}=r_{xi}\exp(\phi\tan(90^{\circ}-\beta))$, where
$\beta=11.5^{\circ}$; the $i$-th region is bounded by the two curves $r_{i-1}$
and $r_{i}$. The coordinates $(r,\phi)$ are defined such that the Galactic
center is (0,0) and the Sun is (-8.5 kpc, 0). The direction of the regular
component is $\hat{b}=\sin(\beta)\hat{r}+\cos(\beta)\hat{\phi}$; the magnitude
of the regular component in the $i$-th region is given by $B=b_{i}{r_{0}}/{r}$
where $r_{0}=5$ kpc.

Figure 6: The polarization degree (left) and survival probability (right) as a
function of the photon energy in the MW galaxy. The blue lines are obtained in
analytic calculations where the medium effects are neglected and the magnetic
field $B_{0}=0.98\ \mu$G is used. The red lines are obtained in numerical
calculations where the medium effects ($\Delta_{x}$ and $\Delta_{y}$) and the
MW magnetic field model are used (see main text).

We consider a source in the direction of $(\ell,b)=(180^{\circ},0^{\circ})$
(the opposite direction of GC) and with a distance of 4 kpc from Earth. In our
analytic calculation we use $B_{0}=0.98\ \mu$G, which is the average value of
the magnitude of the magnetic field along the propagation path. Fig. (6) shows
the the polarization degree and survival probability, computed both in the
analytic calculations and in the numerical calculations where the medium
effects and the MW magnetic field model are used. The agreements between the
analytic and numerical calculations indicate that the analytic calculation
with a uniform magnetic field is a good approximation.

Figure 7: Simulated magnetic fields in the $x$\- (left) and $y$\- (right)
directions in the galaxy cluster as a function of propagation distance $z$.

Figure 8: The polarization degree (left) and survival probability (right) as a
function of the photon energy in the a galaxy cluster. The blue lines are
obtained in analytic calculations where the medium effects are neglected and a
uniform magnetic field of $B_{0}=0.19\ \mu$G is used. The red lines are
obtained in numerical calculations where the medium effects are considered and
the magnetic fields shown in Fig. 7 are used.

###  A.2 Galaxy clusters

We next discuss galaxy clusters. To model the complex magnetic field in galaxy
clusters, we use the parametrization given in Refs. Galanti:2022yxn ;
Bonafede:2010xg

$B_{\rm clu}(z)=B_{\rm clu0}\left(1+\frac{z^{2}}{r_{\rm
core}^{2}}\right)^{-\frac{3}{2}\eta\beta}.$ (33)

For galaxy clusters, we adopt the following parameters: $B_{\rm clu0}=15$
$\mu$G, $r_{\rm core}=100$ kpc, $\beta=2/3$, and $\eta=0.75$ Galanti:2022yxn .
Following Ref. Bonafede:2010xg , we assume that the direction of the magnetic
field is random; hence, we simulate the magnetic field direction by using
Monte Carlo method with a step of 1 kpc along the propagation path. Fig. (7)
shows the simulated magnetic fields in the $x$\- (left) and $y$\- (right)
directions in the galaxy cluster as a function of propagation distance $z$.

We obtain the mean value of the magnetic field in the galaxy cluster by taking
the average of the magnetic field along the propagation path. In our
simulation, we obtain $B_{0}=0.19\ \mu$G, which is the average of the magnetic
field from $z=0$ to $z=500$ kpc. Fig. (8) shows the the polarization degree
and survival probability, computed both in the analytic calculations and in
the numerical calculations where the medium effects are considered and the
magnetic field model given in Eq. (33) is used. We note that although the
analytic calculation deviates from the numerical calculation, it describes the
behavior qualitatively and provide valuable insights for ALP effects.

## References

  * (1) A. Ringwald, “Exploring the Role of Axions and Other WISPs in the Dark Universe,” Phys. Dark Univ. 1 (2012) 116–135 [arXiv:1210.5081].
  * (2) I. G. Irastorza and J. Redondo, “New experimental approaches in the search for axion-like particles,” Prog. Part. Nucl. Phys. 102 (2018) 89–159 [arXiv:1801.08127].
  * (3) C. B. Adams et al. in Snowmass 2021.  2022\.  arXiv:2203.14923.
  * (4) R. D. Peccei and H. R. Quinn, “CP Conservation in the Presence of Instantons,” Phys. Rev. Lett. 38 (1977) 1440–1443.
  * (5) F. Wilczek, “Problem of Strong $P$ and $T$ Invariance in the Presence of Instantons,” Phys. Rev. Lett. 40 (1978) 279–282.
  * (6) J. Jaeckel and A. Ringwald, “The Low-Energy Frontier of Particle Physics,” Ann. Rev. Nucl. Part. Sci. 60 (2010) 405–437 [arXiv:1002.0329].
  * (7) S. Chang, S. Tazawa, and M. Yamaguchi, “Axion model in extra dimensions with TeV scale gravity,” Phys. Rev. D 61 (2000) 084005 [hep-ph/9908515].
  * (8) N. Turok, “Almost Goldstone bosons from extra dimensional gauge theories,” Phys. Rev. Lett. 76 (1996) 1015–1018 [hep-ph/9511238].
  * (9) P. Svrcek and E. Witten, “Axions In String Theory,” JHEP 06 (2006) 051 [hep-th/0605206].
  * (10) A. Ringwald in 4th Patras Workshop on Axions, WIMPs and WISPs, pp. 17–20.  2008\.  arXiv:0810.3106.
  * (11) A. Arvanitaki, S. Dimopoulos, S. Dubovsky, N. Kaloper, and J. March-Russell, “String Axiverse,” Phys. Rev. D 81 (2010) 123530 [arXiv:0905.4720].
  * (12) L. Maiani, R. Petronzio, and E. Zavattini, “Effects of Nearly Massless, Spin Zero Particles on Light Propagation in a Magnetic Field,” Phys. Lett. B 175 (1986) 359–363.
  * (13) G. Raffelt and L. Stodolsky, “Mixing of the Photon with Low Mass Particles,” Phys. Rev. D 37 (1988) 1237.
  * (14) H.-J. Li, J.-G. Guo, X.-J. Bi, S.-J. Lin, and P.-F. Yin, “Limits on axion-like particles from Mrk 421 with 4.5-year period observations by ARGO-YBJ and Fermi-LAT,” Phys. Rev. D 103 (2021) 083003 [arXiv:2008.09464].
  * (15) D. Horns, L. Maccione, A. Mirizzi, and M. Roncadelli, “Probing axion-like particles with the ultraviolet photon polarization from active galactic nuclei in radio galaxies,” Phys. Rev. D 85 (2012) 085021 [arXiv:1203.2184].
  * (16) R. Gill and J. S. Heyl, “Constraining the photon-axion coupling constant with magnetic white dwarfs,” Phys. Rev. D 84 (2011) 085001 [arXiv:1105.2083].
  * (17) X.-J. Bi, et al., “Axion and dark photon limits from Crab Nebula high energy gamma-rays,” Phys. Rev. D 103 (2021) 043018 [arXiv:2002.01796].
  * (18) G. Galanti, “Photon-ALP oscillations inducing modifications to photon polarization,” Phys. Rev. D 107 (2023) 043006 [arXiv:2202.11675].
  * (19) G. Galanti, M. Roncadelli, F. Tavecchio, and E. Costa, “ALP induced polarization effects on photons from galaxy clusters,” Phys. Rev. D 107 (2023) 103007 [arXiv:2202.12286].
  * (20) Z.-Q. Xia, et al., “Searching for the possible signal of the photon-axionlike particle oscillation in the combined GeV and TeV spectra of supernova remnants,” Phys. Rev. D 100 (2019) 123004 [arXiv:1911.08096].
  * (21) Y.-F. Liang, et al., “Constraints on axion-like particle properties with TeV gamma-ray observations of Galactic sources,” JCAP 06 (2019) 042 [arXiv:1804.07186].
  * (22) Y. Grossman, S. Roy, and J. Zupan, “Effects of initial axion production and photon axion oscillation on type Ia supernova dimming,” Phys. Lett. B 543 (2002) 23–28 [hep-ph/0204216].
  * (23) C. Csaki, N. Kaloper, M. Peloso, and J. Terning, “Super GZK photons from photon axion mixing,” JCAP 05 (2003) 005 [hep-ph/0302030].
  * (24) D. Lai and J. Heyl, “Probing Axions with Radiation from Magnetic Stars,” Phys. Rev. D 74 (2006) 123003 [astro-ph/0609775].
  * (25) A. De Angelis, M. Roncadelli, and O. Mansutti, “Evidence for a new light spin-zero boson from cosmological gamma-ray propagation?” Phys. Rev. D 76 (2007) 121301 [arXiv:0707.4312].
  * (26) N. Agarwal, P. Jain, D. W. McKay, and J. P. Ralston, “Signatures of Pseudoscalar Photon Mixing in CMB Radiation,” Phys. Rev. D 78 (2008) 085028 [arXiv:0807.4587].
  * (27) A. K. Ganguly, P. Jain, and S. Mandal, “Photon and axion oscillation in a magnetized medium: A general treatment,” Phys. Rev. D 79 (2009) 115014 [arXiv:0810.4380].
  * (28) A. Mirizzi and D. Montanino, “Stochastic conversions of TeV photons into axion-like particles in extragalactic magnetic fields,” JCAP 12 (2009) 004 [arXiv:0911.0015].
  * (29) C. Wang and D. Lai, “Axion-photon Propagation in Magnetized Universe,” JCAP 06 (2016) 006 [arXiv:1511.03380].
  * (30) A. Kartavtsev, G. Raffelt, and H. Vogel, “Extragalactic photon-ALP conversion at CTA energies,” JCAP 01 (2017) 024 [arXiv:1611.04526].
  * (31) A. De Angelis, G. Galanti, and M. Roncadelli, “Relevance of axion-like particles for very-high-energy astrophysics,” Phys. Rev. D 84 (2011) 105030 [arXiv:1106.1132]. [Erratum: Phys.Rev.D 87, 109903 (2013)].
  * (32) G. Galanti and M. Roncadelli, “Axion-like Particles Implications for High-Energy Astrophysics,” Universe 8 (2022) 253 [arXiv:2205.00940].
  * (33) S. Vernetto and P. Lipari, “Absorption of very high energy gamma rays in the Milky Way,” Phys. Rev. D 94 (2016) 063009 [arXiv:1608.01587].
  * (34) P. Lipari and S. Vernetto, “Diffuse Galactic gamma ray flux at very high energy,” Phys. Rev. D 98 (2018) 043003 [arXiv:1804.10116].
  * (35) W. Dittrich and H. Gies, Probing the quantum vacuum. Perturbative effective action approach in quantum electrodynamics and its application, vol. 166.  2000\.
  * (36) J. I. Latorre, P. Pascual, and R. Tarrach, “Speed of light in nontrivial vacua,” Nucl. Phys. B 437 (1995) 60–82 [hep-th/9408016].
  * (37) M. V. Cougo-Pinto, C. Farina, F. C. Santos, and A. Tort, “The speed of light in confined QED vacuum: Faster or slower than c?” Phys. Lett. B 446 (1999) 170–174.
  * (38) R. Tarrach, “Thermal Effects on the Speed of Light,” Phys. Lett. B 133 (1983) 259–261.
  * (39) S. L. Adler, “Photon splitting and photon dispersion in a strong magnetic field,” Annals Phys. 67 (1971) 599–647.
  * (40) A. Dobrynina, A. Kartavtsev, and G. Raffelt, “Photon-photon dispersion of TeV gamma rays and its role for photon-ALP conversion,” Phys. Rev. D 91 (2015) 083003 [arXiv:1412.4777]. [Erratum: Phys.Rev.D 95, 109905 (2017)].
  * (41) A. Kosowsky, “Introduction to microwave background polarization,” New Astron. Rev. 43 (1999) 157 [astro-ph/9904102].
  * (42) G. B. Rybicki and A. P. Lightman, Radiative Processes in Astrophysics.  1979\.
  * (43) J. M. Yao, R. N. Manchester, and N. Wang, “A new electron density model for estimation of pulsar and FRB distances,” The Astrophysical Journal 835 (2017) 29.
  * (44) I. G. Irastorza, “An introduction to axions and their detection,” SciPost Phys. Lect. Notes 45 (2022) 1 [arXiv:2109.07376].
  * (45) R. Jansson and G. R. Farrar, “A New Model of the Galactic Magnetic Field,” Astrophys. J. 757 (2012) 14 [arXiv:1204.3662].
  * (46) R. Jansson and G. R. Farrar, “The Galactic Magnetic Field,” Astrophys. J. Lett. 761 (2012) L11 [arXiv:1210.7820].
  * (47) R. E. Pudritz, M. J. Hardcastle, and D. C. Gabuzda, “Magnetic Fields in Astrophysical Jets: From Launch to Termination,” Space Sci. Rev. 169 (2012) 27–72 [arXiv:1205.2073].
  * (48) F. Tavecchio, G. Ghisellini, G. Ghirlanda, L. Foschini, and L. Maraschi, “TeV BL Lac objects at the dawn of the Fermi era,” Mon. Not. Roy. Astron. Soc. 401 (2010) 1570 [arXiv:0909.0651].
  * (49) F. Govoni and L. Feretti, “Magnetic field in clusters of galaxies,” Int. J. Mod. Phys. D 13 (2004) 1549–1594 [astro-ph/0410182].
  * (50) A. Neronov and I. Vovk, “Evidence for strong extragalactic magnetic fields from Fermi observations of TeV blazars,” Science 328 (2010) 73–75 [arXiv:1006.3504].
  * (51) M. S. Pshirkov, P. G. Tinyakov, and F. R. Urban, “New limits on extragalactic magnetic fields from rotation measures,” Phys. Rev. Lett. 116 (2016) 191302 [arXiv:1504.06546].
  * (52) C. Y. Yang, et al. in Space Telescopes and Instrumentation 2018: Ultraviolet to Gamma Ray, J.-W. A. den Herder, S. Nikzad, and K. Nakazawa, eds., vol. 10699 of Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, p. 106992K.  2018\.
  * (53) K. Kawata, et al., “Energy determination of gamma-ray induced air showers observed by an extensive air shower array,” Exper. Astron. 44 (2017) 1–9.
  * (54) G. Z. Angeli and P. Dierickx, eds., “The e-ASTROGAM gamma-ray space observatory for the multimessenger astronomy of the 2030s,” Proc. SPIE Int. Soc. Opt. Eng. 10699 (2018) 106992J [arXiv:1805.06435].
  * (55) AMEGO Team Collaboration, “AMEGO: Exploring the Extreme Multimessenger Universe,” Proc. SPIE Int. Soc. Opt. Eng. 11444 (2020) 1144431 [arXiv:2101.03105].
  * (56) A. Bonafede, et al., “The Coma cluster magnetic field from Faraday rotation measures,” Astron. Astrophys. 513 (2010) A30 [arXiv:1002.0594].


# Fermi-Surface Modeling of Light-Rare-Earth Hexaborides
with 2D-ACAR Spectroscopy

###### Abstract

Two dimensional angular correlation of the positron annihilation radiation
(2D-ACAR) spectra are measured for $\mathrm{LaB}_{6}$ along high symmetry
directions and compared with first principle calculations based on density
functional theory (DFT). This allows the modeling of the Fermi surface in
terms of ellipsoid electron pockets centered at $X$-points elongated along the
$\Sigma$ axis (${\Gamma-M}$ direction). The obtained structure is in agreement
with quantum oscillation measurements and previous band structure
calculations. For the isostructural topologically not-trivial
$\mathrm{SmB}_{6}$ the similar ellipsoids are connected through necks that
have significantly smaller radii in the case of $\mathrm{LaB}_{6}$. A
theoretical analysis of the 2D-ACAR spectra is also performed for
$\mathrm{CeB}_{6}$ including the on-site repulsion $U$ correction to the
local-density approximation (LDA+$U$) of the DFT. The similarities of 2D-ACAR
spectra and the Fermi-surface projections of these two compounds allow to
infer that both $\mathrm{LaB}_{6}$ and $\mathrm{CeB}_{6}$ are topologically
trivial correlated metals.

###### keywords:

positron annihilation, electronic structure, rare-earth hexaborides

Josef Ketels* Michael Leitner Peter Böni Christoph Hugenschmidt*
Mikheil Sekania* Alyn D. N. James Jakob A. E. Bonart Nico Unglert Liviu
Chioncel*

J. Ketels, P. Böni
Physik-Department, Technische Universität München, James-Franck-Straße 1,
85748 Garching, Germany
Email Address<EMAIL_ADDRESS>
M. Leitner, C. Hugenschmidt
Physik-Department, Technische Universität München, James-Franck-Straße 1,
85748 Garching, Germany
Heinz Maier-Leibnitz Zentrum (MLZ), Technische Universität München,
Lichtenbergstraße 1, Garching, 85748, Germany
Email Address<EMAIL_ADDRESS>
M. Sekania
Theoretical Physics III, Center for Electronic Correlations and Magnetism,
Institute of Physics, University of Augsburg, Universitätsstraße 1, 86159
Augsburg, Germany
Institut für Physik, Martin-Luther Universität Halle-Wittenberg, 06120
Halle/Saale, Germany
Andronikashvili Institute of Physics, Javakhishvili Tbilisi State University,
Tamarashvili str. 6, 0177 Tbilisi, Georgia
Email Address<EMAIL_ADDRESS>
A. D. N. James
H. H. Wills Physics Laboratory, University of Bristol, Tyndall Avenue,
Bristol, BS8 1TL, United Kingdom

N. Unglert, J. A. E. Bonart, L. Chioncel
Theoretical Physics III, Center for Electronic Correlations and Magnetism,
Institute of Physics, University of Augsburg, Universitätsstraße 1, 86159
Augsburg, Germany
Augsburg Center for Innovative Technologies, University of Augsburg,
Universitätsstraße 1, 86159 Augsburg, Germany
Email Address<EMAIL_ADDRESS>
##  1 Introduction

The rare-earth hexaborides ($R\mathrm{B_{6}}$, $R=\mathrm{La}$, $\mathrm{Ce}$,
$\mathrm{Pr}$, $\mathrm{Nd}$, and $\mathrm{Sm}$) provide an interesting
subject for experimental and theoretical studies because of a variety of
features that can be ascribed to the particular role played by their
$4f$-electrons [1, 2]. The unfilled $4f$-states are a challenging problem for
an accurate theoretical description of the electronic structure of rare-earth
compounds. The $4f$ energy levels often overlap with the non-$4f$ broad bands,
form narrow resonances, but nevertheless are frequently treated as core
states. A proper description, however, requires inclusion of the $4f$-states
in valence bands where these states are also subjected to a strong on-site
Coulomb repulsion. A simple theory that can capture this strong on-site
Coulomb repulsion for $4f$ orbitals at the mean-field level is the LDA+$U$
approach [3]. Another important ingredient in the physics of the rare-earth
compounds is the presence of significant spin-orbit coupling, which
fortunately can be included on an equal footing with the strong on-site
Coulomb interactions (the Hubbard $U$). The recent theoretical models beyond
the local approximations, such as $GW$ [4], or dynamic electronic correlations
such as Dynamical Mean Field Theory (DMFT) [5, 6], however, are difficult to
use in case of the multi-orbitals (like $4f$) and the strong spin-orbit
couplings.

In the $R\mathrm{B}_{6}$ series of compounds, $\mathrm{LaB}_{6}$ is considered
as the reference non $f$-electron system, $\mathrm{CeB}_{6}$ is a Kondo
lattice heavy fermion system, while $\mathrm{PrB}_{6}$ and $\mathrm{NdB}_{6}$
are local-moment magnetically ordered metals. Among these compounds
$\mathrm{SmB}_{6}$ is considered to be a “topological Kondo insulator” [7] due
to the narrow hybridization band gap that opens at the Fermi level. In
general, some of the $4f$-bands of the rare-earth elements are located around
the Fermi level ($E_{\mathrm{F}}$) and are crossed by the more dispersive $5d$
band. The finite spin-orbit coupling splits these crossings and opens a gap.
The band-character inversion due to the spin-orbit coupling between the rare-
earth $4f$/$5d$ states situated around $E_{\mathrm{F}}$ signals the presence
of a topologically non-trivial phase. This type of behavior was observed in
DFT calculations for $\mathrm{SmB}_{6}$ [8] and further searched in other
rare-earth hexaborides [9]. The effects due to a finite Hubbard interaction,
$U$, were also investigated for $\mathrm{Sm}$-$4f$ and it was found that the
topological features in $\mathrm{SmB}_{6}$ are insensitive to the values of
the on-site Coulomb interaction. The latter, however, plays a crucial role in
the other isostructural hexaborides [9].

In the present study we investigate the electronic properties of two members
of the $R\mathrm{B}_{6}$ family, namely $\mathrm{LaB}_{6}$ and
$\mathrm{CeB}_{6}$, with two dimensional angular correlation of annihilation
radiation (2D-ACAR) and first principle Density Functional Theory (DFT) [10,
11, 12, 13]. We analyze the radial anisotropies of the obtained data in the
extended momentum space and the corresponding back-folded spectra in the
crystal-momentum space in order to identify the different Fermi surface
features. We also compare the experimentally measured and computed 2D-ACAR
spectra for $\mathrm{LaB}_{6}$. The crystal structures of the
$R\mathrm{B}_{6}$ belong to the simple cubic $\mathrm{CsCl}$ structure type
($Pm\overline{3}m$), with lattice constant
${a=$4.1569\text{\,}\mathrm{\SIUnitSymbolAngstrom}$}$ and
${a=$4.1391\text{\,}\mathrm{\SIUnitSymbolAngstrom}$}$, for $\mathrm{LaB}_{6}$
and $\mathrm{CeB}_{6}$, respectively. The rare-earth atoms occupy the corner
of the unit cell, corresponding to the ${1a\,(0,0,0)}$ Wyckoff site, while the
B atoms are located at the octahedral sites in the body-centered position at
the ${6f\,(0.5,0.5,z)}$ Wyckoff sites, where $z$ is $0.2011$ for
$\mathrm{CeB}_{6}$ [14] and $0.1996$ for $\mathrm{LaB}_{6}$ [15].

Although the $f$-shell of $\mathrm{La}$ is empty while $\mathrm{Ce}$ has a
single $f$-electron and $\mathrm{CeB}_{6}$ exhibits a far richer phase diagram
than $\mathrm{LaB}_{6}$, undergoing several different magnetic phase
transitions, de Haas-Van Alphen (dHvA) measurements [16] yielded similar Fermi
surface (FS) features for both $\mathrm{CeB}_{6}$ and $\mathrm{LaB}_{6}$. The
FS of $\mathrm{CeB}_{6}$ consists of large “ellipsoid” pockets around the $X$
points of the Brillouin zone (BZ), whereas the states around the zone center
($\Gamma$-point) are shifted away by band-renormalization effects leading to a
hole pocket [17]. The first-principles calculations agree with the
experimental results around the $X$-points but fail to capture the strongly
renormalized electronic states around the $\Gamma$-point [17]. The dHvA
measurements [16] at temperatures ${T<T_{K}}$ can be brought into accordance
with a model in which $f$-electrons do not contribute to the Fermi volume.
Above the Kondo temperature ${T>T_{K}}$, in the paramagnetic phase, dHvA
measurements and 2D-ACAR results were also in agreement in the case of
$\mathrm{CeB}_{6}$ [16]. Results of the 2D-ACAR experiments on
$\mathrm{LaB}_{6}$ were reported earlier in the literature [18, 19], where the
analysis in terms of the back-folded momentum densities were performed. Three-
dimensional reconstruction of the Fermi surface using 2D-ACAR projections were
also reported [20], and contrary to the electronic structure calculations the
strong hybridization between $\mathrm{B}$-$p$ and $\mathrm{La}$-$d$ states
($f$-states are unoccupied) that produced additional FS sheets, were not
detected in the experiment.

In our work we analyze the momentum densities of $\mathrm{LaB}_{6}$
(experiment/theory) both in the extended $p$ and reduced $k$ spaces. Possible
similarities between the theoretical spectra with the isostructural
$\mathrm{CeB}_{6}$ are also discussed. For the $\mathrm{CeB}_{6}$ compound, in
the absence of recent experimental measurements, we theoretically study the
effects of the different $U$ values of the Hubbard interaction and diverse
double counting schemes of the LDA+$U$ method. Based on our 2D-ACAR results,
we conclude that both $\mathrm{LaB}_{6}$ and $\mathrm{CeB}_{6}$ are
topologically trivial but correlated metals. Therefore, we believe that the
future studies with DMFT which takes into account the full local correlations
starting from itinerant (valence) $4f$-states will provide a suitable
description of the electronic structure of the rare-earth hexaborides.

##  2 Methods

###  2.1 2D-ACAR

2D-ACAR is a powerful tool to investigate the bulk electronic structure [21,
22]. It is based on the annihilation of positrons with electrons of a sample
leading to the emission of two $\gamma$-quanta in nearly anti-parallel
directions. The small angular deviation from collinearity is caused by the
transverse component of the electron’s momentum. Hence, the coincident
measurement of the annihilation quanta for many annihilation events yields a
projection of the so called two photon momentum density (TPMD)
$\rho^{2\gamma}(\mathrm{\bf{p}})$. This is usually computed as the Fourier
transform of the product of positron wave function $\Psi^{+}(\mathrm{\bf{x}})$
and electron wave function $\Psi^{-}({\mathrm{\bf{x}}})$:

$\rho^{2\gamma}({\mathrm{\bf{p}}})\propto\sum_{j,\mathrm{\bf{k}}}n_{j}({\mathrm{\bf{k}}})\left|\int\mathrm{d}\mathrm{\bf{x}}\,e^{-i{\mathrm{\bf{xp}}}}\,\Psi^{+}({\mathrm{\bf{x}}})\Psi^{-}({\mathrm{\bf{x}}})\,\sqrt{\gamma({\mathrm{\bf{x}}})}\right|^{2}$
(1)

The sum runs over all $\mathrm{\bf{k}}$ states in all bands, $j$, with the
occupation $n_{j}(\mathrm{\bf{k}})$. The so-called “enhancement factor”
$\gamma({\mathrm{\bf{x}}})$ [23], takes into account the electron positron
correlation. The 2D-ACAR spectrum
$\rho^{\vphantom{\dagger}}_{\mathrm{ACAR}}(p_{x},p_{y})$, the quantity which
is actually accessible by an experiment, is a 2D projection of the 3D
momentum-density distribution $\rho^{2\gamma}(\mathrm{\bf{p}})$ along a chosen
($p_{z}$) axis

$\rho^{\vphantom{\dagger}}_{\mathrm{ACAR}}(p_{x},p_{y})=\int\rho^{2\gamma}({\mathrm{\bf{p}}})\mathrm{d}p_{z}\,.$
(2)

The positron annihilation probes all electrons in the system. Filled bands,
especially bands of core electrons, give a nearly isotropic distribution which
is superimposed by an anisotropic contribution mainly produced by the
electrons near the Fermi surface. This anisotropic
$\rho^{\vphantom{\dagger}}_{\mathrm{aniso}}(p_{x},p_{y})$ contribution is
therefore one of the most interesting feature of an ACAR spectrum
$\rho^{\vphantom{\dagger}}_{\mathrm{ACAR}}(p_{x},p_{y})$. It can be calculated
by subtracting the isotropic contribution
${\overline{\rho}^{\vphantom{\dagger}}_{\mathrm{ACAR}}(p_{x},p_{y})}$:

$\rho^{\vphantom{\dagger}}_{\mathrm{aniso}}(p_{x},p_{y})=\rho^{\vphantom{\dagger}}_{\mathrm{ACAR}}(p_{x},p_{y})-\overline{\rho}^{\vphantom{\dagger}}_{\mathrm{ACAR}}(p_{x},p_{y})\,,$
(3)

where the isotropic contribution is the radial average
${\overline{\rho}^{\vphantom{\dagger}}_{\mathrm{ACAR}}(p_{x},p_{y})\equiv\overline{\rho}^{\vphantom{\dagger}}_{\mathrm{ACAR}}(\sqrt{p_{x}^{2}+p_{y}^{2}})}$
constructed from the original spectrum
$\rho^{\vphantom{\dagger}}_{\mathrm{ACAR}}(p_{x},p_{y})$ by averaging over all
data points in equidistant intervals ${[p_{r},p_{r}+\Delta p_{r})}$ from the
center.

The discontinuities in $\rho^{2\gamma}(\mathrm{\bf{p}})$ correspond to the
density jumps in the Fermi distribution and are used to identify the Fermi
surface sheets. These discontinuities occurring at
${\mathrm{\bf{p}}=\mathrm{\bf{k}}+\mathrm{\bf{G}}}$, where $\mathrm{\bf{k}}$
is the wave vector in the first BZ and $\mathrm{\bf{G}}$ is a reciprocal
lattice vector, can be folded back from the extended zone scheme $p$-space to
the reduced zone scheme $k$-space by applying the so called Lock-Crisp-West
(LCW) theorem [24].

###  2.2 Experiment Details

The ACAR experiments were performed at the 2D-ACAR spectrometer at the
Technische Universität München. The detector-detector distance is
$17.5\text{\,}\mathrm{m}$ with the sample positioned exactly in the middle.
This leads to an angular resolution of
${\sigma_{h}=$0.538(1)\text{\,}\mathrm{m}\mathrm{r}\mathrm{a}\mathrm{d}$}$ and
${\sigma_{v}=$0.655(1)\text{\,}\mathrm{m}\mathrm{r}\mathrm{a}\mathrm{d}$}$ in
horizontal and vertical direction, respectively [25]. All measurements were
conducted at a temperature of $15\text{\,}\mathrm{K}$ to reduce the
contribution of the thermal motion of the positron. The positrons were guided
onto the sample by a magnetic field of $1.2\text{\,}\mathrm{T}$ at the sample
position. For full details on the experimental setup we refer to ref [26].

The $\mathrm{LaB}_{6}$ sample used in the 2D-ACAR experiment was a cuboid-
shaped high-quality single crystal of ${10\times 5\times 2}$
${\mathrm{mm}}^{3}$ size. It was grown by vertical crucible-free inductive
floating zone melting in argon gas atmosphere. The grown crystal was
characterized by Laue back-scattering, optical spectral analysis, X-ray
diffraction and density measurements. The surfaces of the plate were polished
with diamond powder.

We recorded 2D-ACAR projections along the high symmetry $[110]$, $[001]$ and
$[111]$ directions within the $(1\bar{1}0)$ plane. For each spectrum more than
$9\text{\times}{10}^{7}\text{\,}\mathrm{c}\mathrm{o}\mathrm{u}\mathrm{n}\mathrm{t}\mathrm{s}$
were collected. To correct for the varying efficiency and limited aperture of
the Anger camera detectors, the raw data were divided by the so-called
momentum sampling function (MSF). To calculate the MSF the single events on
both detectors were collected in parallel to the coincident events throughout
the experiment. The convolution of those single spectra gives the MSF. The
ACAR spectra were symmetrized according to the crystal point-symmetry to
further enhance the statistics.

###  2.3 Electronic Structure Calculation and Enhancement Models

Electronic structure calculations were performed using the ELK code [27],
which is an all-electron full-potential linearised augmented plane-wave (LAPW)
code for determining the properties of crystalline solids. The valence
electron configuration for the rare-earth atoms ($\mathrm{La/Ce}$) is
${6s6p5d4f}$, while the valence electrons of B are located in the $2s$ and
$2p$ orbitals. The self-consistent calculations were performed using the Local
Density Approximation (LDA) exchange-correlation functional parameterized by
Perdew and Wang [28]. For the LDA+$U$ the rationally invariant formulation of
the Coulomb interaction between the valence electrons was used [3] with
different double counting schemes [29]. Further details of the LDA+$U$
calculations can be found in Section 3.2. In the self-consistent calculations
the $\mathrm{\bf{k}}$-summations are performed using the tetrahedron method,
with a $20\times 20\times 20$ $k$-mesh in the irreducible part of the
Brillouin zone and convergence has been achieved in total energy with an
accuracy of better than $10^{-5}$ Ry. For the detailed Fermi surface
calculations a significantly larger, $80\times 80\times 80$, $k$-mesh was
employed.

The DFT can be generalized to electron-positron systems by including the
positron density, in the form of the two-component DFT [30, 31]. In this
formalism the positron is considered to be thermalized and described by a
state with ${{\mathrm{\bf{p}}}_{p}=0}$ with $s$-type symmetry at the bottom of
the positronic band. The photons resulting from the electron-positron
annihilation carry the momentum of the annihilated pair up to a reciprocal
lattice vector, reflecting the fact that the annihilation takes place in a
crystal. Hence, an electron with a transverse wave vector ${\mathrm{\bf{k}}}$
contributes to spin-resolved two-photon momentum density,
$\rho^{2\gamma}(\textbf{p})$ not only at
${{\mathrm{\bf{p}}}={\mathrm{\bf{k}}}}$ (normal process) but also at
${{\mathrm{\bf{p}}}={\mathrm{\bf{k}}}+{\mathrm{\bf{G}}}}$, with
${\mathrm{\bf{G}}}$ a vector of the reciprocal lattice (Umklapp process). In
the LDA(+$U$) framework the electron-positron momentum density
$\rho({\mathrm{\bf{p}}})$ is computed directly with the spin-resolved versions
of Equations 1 and 2. We calculate the 2D-ACAR spectra according to the method
described in ref [32]. These TPMD were calculated to a maximum momentum of $6$
a.u. and the electron positron correlations are taken into account within the
$\gamma(\mathrm{\bf{x}})$ term, the so-called enhancement factor. Within the
Independent Particle Model (IPM) ${\gamma(\mathrm{\bf{x}})=1}$, while models
of enhancement beyond IPM consider momentum and energy dependencies and a
separate treatment of $f$, $d$ states in comparison with the $s$, $p$ states.
These have been quantified via many different approximations [33]. In the
present studies we employ the so called Drummond parametrization [34, 35].

Figure 1: Radial anisotropy of measured and symmetrized $\mathrm{LaB}_{6}$
ACAR spectra (top) and theoretical spectra calculated in the IPM (bottom left)
and with the Drummond enhancement (bottom right). The spectra from first-
principle calculations are convolved with a Gaussian accounting for the
experimental resolution of $0.619\text{\,}\mathrm{mrad}$. The borders of the
projected first BZ (grey lines) in a repeated zone scheme are shown in all
plots. From left to right $(110)$, $(001)$, and $(111)$ high symmetry
projections.

##  3 Results and Discussion

###  3.1 $\mathrm{LaB}_{6}$: Experiment and Theory

To compare experimental and theoretical data on $\mathrm{LaB}_{6}$ in the
$p$-space we look at the radial anisotropy
$\rho^{\vphantom{\dagger}}_{\mathrm{aniso}}(p_{x},p_{y})$, which emphasizes
the spectral contribution from electrons near the Fermi level. Figure 1 shows
the radial anisotropy for all three measured projections (from left to right
$(110)$, $(001)$ and $(111)$) in the upper half and the corresponding
theoretical calculations in the lower half of the plots. The calculated
results are for the momentum distributions with IPM (left) and Drummond
enhancement (right). To account for the limited experimental resolution the
first-principle spectra are convolved with a Gaussian of width
${\sigma^{\vphantom{\dagger}}_{\mathrm{exp}}=$0.619\text{\,}\mathrm{mrad}$}$.

Overall, we have a very good agreement between theory and experiment. Dominant
features seen in the theoretical calculations are also visible in the
experimental data. Especially the four-fold and sixfold symmetric projections
are well reproduced. Some minor differences can be seen in the $(110)$
projection. Particularly, the weight distribution in the region between $-5$
to $5\text{\,}\mathrm{mrad}$ on both axes shows a significant departure
between the theory and experiment. However, the high intensity positive signal
along ${[1\bar{1}0]}$ direction at about $8\text{\,}\mathrm{mrad}$ is well
reproduced. We attribute the better agreement between the theory and the
experiment for the $(001)$ and $(111)$ projections to the higher point
symmetry of these directions. In all experimental projections the anisotropy
profiles are slightly stretched out within the plane. The similar behaviour is
also observed when the Drummond enhancement is included in the computation as
compared to the IPM results. However, the radial expansion of the profile is
not sufficient to match the observed “stretching” in the experimental
measurement. The differences between the IPM and Drummond model are marginal,
though.

Figure 2: $\mathrm{LaB}_{6}$ LCW back-folded $k$-space density of the ACAR
measurements (top) and the theoretical calculations in the IPM (bottom left)
and the Drummond model (bottom right). The spectra from first-principle
calculations are convolved with a Gaussian accounting for the experimental
resolution of $0.619\text{\,}\mathrm{mrad}$. The borders of the projected
first BZ (white lines) in a repeated zone scheme are shown in all plots. From
left to right $(110)$, $(001)$, and $(111)$ high symmetry projections.

Figure 2 shows the LCW back-folded data, experiments as well as theory, and is
organized according to the structure of Figure 1. The top part of each plot
shows the $k$-space experimental 2D-ACAR spectrum, while the bottom part is
split between the calculated results using IPM (left) and Drummond enhancement
(right), broadened with the experimental resolution.

The measured and calculated LCW back-folded spectra show an overall good
agreement. The $(001)$ and $(111)$ projections show the expected features of
the Fermi surface topology (see Figure 3). In the $(111)$ projection, the
sixfold rose structure around the $\Gamma$-point becomes apparent. As expected
from the calculated Fermi surface, there is no density at the $\Gamma$\- and
$R$-point in the BZ. Note that one can not see the hole-space around the
$\Gamma$-point in the $(001)$ projection, as they are covered by the
“ellipsoids” along the $[001]$ direction. While the theoretical calculations
explicitly reproduce the “ellipsoids” structure in $(110)$ projection, the
experimental $(110)$ projection shows stripes along $[1\bar{1}0]$ direction.
The latter can be also attributed to the observed substantial discrepancies
between the theory and experiment in the radial anisotropy (cf. Figure 1).

The Fermi surface of $\mathrm{LaB}_{6}$ is shown in Figure 3 in the simple
cubic BZ where the high symmetry points are ${\Gamma\equiv[0,0,0]}$,
${X\equiv[1/2,0,0]}$, ${M\equiv[1/2,1/2,0]}$ and ${R\equiv[1/2,1/2,1/2]}$. The
Fermi surface of $\mathrm{LaB}_{6}$ consists of a set of equivalent
“ellipsoids” centered at the $X$-points and connected by necks which intersect
along $\Sigma$ ($\Gamma-M$) direction. In the panels (b) and (c) in Figure 3
we show a cross-section in the ${\Gamma-X-M}$-plane (${k_{[100]}=0}$ plane).
The DFT (LDA) Fermi surface calculation for $\mathrm{SmB}_{6}$ with a downward
shift of the Fermi energy is presented in ref [36] and has a FS similar to
that of $\mathrm{LaB}_{6}$ as shown in Figure 3). The angular variations of
the dHvA frequencies and the disappearance of oscillations in some angular
regions have been recently discussed for both $\mathrm{SmB}_{6}$ and
$\mathrm{LaB}_{6}$ [36]. These quantum oscillation frequencies identify the
so-called $\alpha,\gamma,\varepsilon$-branches associated to FS regions
visible in the ${\Gamma-X-M}$-plane, see Figure 3. The $\alpha$-branches were
associated with the ellipsoids centered at the $X$-point. As one can see the
nearest-neighbor ellipsoid FS pieces touch along the ${\Gamma-M}$ line and
connect through the small distorted circular shape neck building up a multiply
connected FS. Both the $\gamma$ and $\varepsilon$-orbits are hole-like orbits
and are centered around $M$-point and $\Gamma$-point, respectively. According
to the dHvA frequencies the angular region for the hole-like $\gamma$\- and
$\varepsilon$-orbits is significantly larger in $\mathrm{SmB}_{6}$ than in
$\mathrm{LaB}_{6}$ which may be the cause of the different physical properties
of these two compounds. From our calculations we can compare linear dimensions
of the Fermi surface features to both previous computations and the present
experimental data. The corresponding values are $0.848$, $0.668$, and $0.624$
for the ${X-\Gamma}$, ${X-M}$, and ${X-R}$ directions, respectively, as a
fraction of the BZ size. Our results agree with the previous results reported
in the literature [20].

To determine the FS linear dimension from the experimental data [20, 37],
taking into account a finite experimental resolution, we model the electron-
momentum density with prolate ellipsoids with equatorial radius
$R_{\mathrm{eq}}$ and polar radius $R_{\mathrm{pl}}$, centered at the $X$
points in reciprocal space and pointing towards the neighbouring $\Gamma$
points

$F(\mathrm{\bf{k}};R_{\mathrm{eq}},R_{\mathrm{pl}})=\sum\limits_{\mathrm{\bf{G}}}\sum\limits_{\mathrm{\bf{b}}_{i}}\Theta\left(1-\frac{|(\mathrm{\bf{k}}-\mathrm{\bf{G}})\times{\mathrm{\bf{b}}_{i}}|^{2}}{R_{\mathrm{eq}}^{2}|\mathrm{\bf{b}}_{i}|^{2}}-\frac{|(\mathrm{\bf{k}}-\mathrm{\bf{G}}-{\mathrm{\bf{b}}_{i}}/2)\cdot{\mathrm{\bf{b}}_{i}}|^{2}}{R_{\mathrm{pl}}^{2}|\mathrm{\bf{b}}_{i}|^{2}}\right)\,.$
(4)

Here, $\Theta(\cdot)$ is the Heaviside step function and summation is over
reciprocal lattice vectors $\mathrm{\bf{G}}$ and reciprocal primitive vectors
${\mathrm{\bf{b}}_{i}\in\\{\mathrm{\bf{b}}_{[100]},\mathrm{\bf{b}}_{[010]},\mathrm{\bf{b}}_{[001]}\\}}$.
Adding a flat background two-photon momentum density due to fully occupied
bands described by a third parameter $F_{\mathrm{backg}}$, projecting along
the given direction and convolving with the anisotropic Gaussian kernel that
models the experimental resolution, we fit the resulting simulated projections
to the correspondingly normalized LCW back-folded data for all measured
projection directions by varying the three free parameters. We obtain
${2R_{\mathrm{eq}}=0.662(2)}$ and ${2R_{\mathrm{pl}}=0.952(2)}$ in units of
the reciprocal lattice constant, where the reported uncertainties are the
statistical precision due to counting noise. We see that the equatorial
diameter $2R_{\mathrm{eq}}$ agrees well with the calculated ${X-M}$ length,
while the length along the ${X-\Gamma}$ direction is overestimated. We think
that this is due to positron wave-function and/or electron-positron
correlation effects, which are neglected in the LCW theorem [24]: While the
equatorial diameter is directly accessible in particular in the (001)
projection, where the contrast between the filled ellipsoids and the
background density is large, there is no projection direction that would
afford an unobstructed view on the polar termination of the ellipsoids. Thus,
a larger band weight of the conduction and/or core bands towards $\Gamma$
would necessarily result in an overestimated $R_{\text{pl}}$. Also, with the
current experimental resolution we cannot detect possible deviations of the
Fermi surface from the shape of a rotational ellipsoid in the vicinity of the
poles, or whether the ${X-M-R}$ cross-section of the FS shows four-fold
deviations from circular symmetry, as shown by the calculations (see Figure
3).

Figure 3: Fermi surface in the first BZ (a),(d) and ${\Gamma-X-M}$ cross-
section (${[010]-[001]}$ cut) (b),(c) from first-principle calculations for
$\mathrm{LaB}_{6}$ (a),(b) and for $\mathrm{CeB}_{6}$ (c),(d), respectively.
FS sheets, for two different spin projections, are depicted in the left and
the right half of the plots (c) and (d), respectively. The borders of the
first BZ (grey lines) and high-symmetry points ($\Gamma$, $X$, $M$, and $R$)
are shown as well.

###  3.2 $\mathrm{CeB}_{6}$: Theory

Among the most frequently used techniques to include short-range Coulomb
interactions between the electrons in the framework of DFT are the self-
interaction-corrected local-spin-density approximation (SIC-LSDA) [38] and
LDA+$U$ [3]. The LDA+$U$ method is particularly useful as it is
computationally less demanding and hence can be easily used for systems with
unit cells consisting of a considerable number of atoms. Furthermore, its
flexibility allows a systematic study of effects of the on-site repulsion $U$,
on the ground state properties such as equilibrium lattice parameter, magnetic
moment, and — in the current work — 2D-ACAR spectra. The orbital dependent
LDA+$U$ functional used in the present work is given by [3]:

$E_{\mathrm{LDA}+U}=E_{\mathrm{LDA}}+\frac{U-J}{2}\sum_{\sigma}[\tr\rho^{\sigma}-\tr(\rho^{\sigma}\rho^{\sigma})]$
(5)

Here, $\rho^{\sigma}$ (not to be confused with momentum densities) is the
density matrix for the $f$-states and $U$ and $J$ are the local Coulomb and
exchange Hund’s parameters. The lack of a diagrammatic expansion of the DFT
total energy makes it difficult to model the effect of the local Coulomb
interactions beyond the effects already captured by the exchange-correlation
functionals. To avoid the double counting of such effects several schemes have
been proposed for different kind of materials. One of them is the mean-field
approximation to the Hubbard correction, the so-called “fully-localized” limit
(FLL). In this scheme each localized (e.g., atomic) orbital is either full or
completely empty. This formulation of the double-counting term mimics an
expansion of the electronic energy around the strongly localized limit and
thus tends to work quite well for strongly correlated materials with localized
orbitals. For other systems, such as metals or “weakly correlated” materials,
the excessive stabilization of occupied states due to the “+$U$” potential may
lead to unphysical results, for example the enhancement of the Stoner factor
[29]. A different double-counting scheme, the “around mean-field” (AMF), was
introduced in ref [39] and further developed in ref [29]. In the present study
the AMF double counting scheme failed to produce the correct position of the
$4f$ bands, therefore FLL double-counting was employed in all presented
results.

Figure 4: Unconvolved radial anisotropy of the calculated ACAR spectra for
$\mathrm{LaB}_{6}$ (upper left) and $\mathrm{CeB}_{6}$ for ${U=0}$ (upper
right), ${U=4}$ eV (bottom left), and ${U=6}$ eV (bottom right). Results
correspond to LDA ($\mathrm{LaB}_{6}$ and $\mathrm{CeB}_{6}$ ${U=0}$) and
LDA+$U$ (the rest) calculations with the Drummond parametrization of the
electron-positron enhancement factor. The borders of the projected first BZ
(grey lines) in the repeated zone scheme are shown in all plots. From left to
right $(110)$, $(001)$, and $(111)$ high symmetry projections.

The DFT band structure of $\mathrm{CeB}_{6}$ produces the manifold of
$\mathrm{Ce}$-$4f$ orbitals in the close vicinity of the Fermi level.
According to the experimental study of Neupane et al. [17] they have to be
located at around $2$ eV in the conduction band. Therefore, for the FLL limit
a reasonable parameter for the effective Coulomb interaction would be $U=4$
eV. The band structure and the density of states (not shown) in the
ferromagnetic states are characterized by dispersive $5d$\- and flat
$4f$-bands along the ${M-X-M}$ and ${X-\Gamma-X}$ directions. The flat bands
are purely build from $\mathrm{Ce}$-$4f$ states. The dispersive $5d$ band
around the $X$-point are situated $2$ eV below $E_{\mathrm{F}}$ and touches
the $\mathrm{B}-2p$ states at the $\Gamma$ point. The position of these bands
agrees fairly well with experiments [17, 40].

In Figure 4 we present results for the calculated radial anisotropy of the
2D-ACAR spectra compared with $\mathrm{CeB}_{6}$. The electronic calculations
for $\mathrm{LaB}_{6}$ were done using DFT only. For $\mathrm{CeB}_{6}$ we
used LDA+$U$ by varying the value of $U$ between $0$ and $6$ eV. We also show
the corresponding LCW back-folded results in Figure 5. One can clearly see the
similarity between the $\mathrm{LaB}_{6}$ and $\mathrm{CeB}_{6}$, $U=4$ and
$6$ eV spectra. This can be attributed to the on-site Coulomb repulsion
(Hubbard $U$), which pushes unoccupied $4f$ -manifolds above the Fermi level
and localizes remained filled $4f$-orbitals. The localization of a single
electron in the $4f$-orbital causes the fully ferromagnetic ground state with
the magnetic moment of about $1\mu_{B}$.

Figure 5: Unconvolved LCW folded back $k$-space density of the calculated ACAR
spectra for $\mathrm{LaB}_{6}$ (upper left) and $\mathrm{CeB}_{6}$ for ${U=0}$
(upper right), ${U=4}$ eV (bottom left), and ${U=6}$ eV (bottom right).
Results correspond to LDA ($\mathrm{LaB}_{6}$ and $\mathrm{CeB}_{6}$ ${U=0}$)
and LDA+$U$ (the rest) calculations with the Drummond parametrization of the
electron-positron enhancement factor. The borders of the projected first BZ
(white lines) in the repeated zone scheme are shown in all plots. From left to
right $(110)$, $(001)$, and $(111)$ high symmetry projections.

The shown spectra agree with the 2D-ACAR-measurements presented in ref [19]
for both structures ($\mathrm{LaB}_{6}$ and $\mathrm{CeB}_{6}$). We observe
that the $\mathrm{CeB}_{6}$ spectra becomes similar to $\mathrm{LaB}_{6}$ for
$U=4$ and $6$ eV in all shown high-symmetry projections. This effect is
connected to the localization of the $4f$-orbital as its position in energy is
shifted further below the Fermi surface with increasing values of the Hubbard
$U$ parameter. On the central panel of the Figure 5 ($(001)$ projection),
where we show the LCW back-folded unconvolved data, one can clearly identify
the Fermi surface “ellipsoids” in ${[100]}$ and ${[010]}$ directions as well
as rounded-square cross-sections of the same ellipsoids along the ${[001]}$
direction. This demonstrates that the high-resolution 2D-ACAR spectra can be
used to deduce bulk FS features of the $R\mathrm{B}_{6}$ family of compounds,
which also including $4f$-valence electron bands.

##  4 Summary and Outlook

The 2D-ACAR experiments offer answers to some fundamental questions whether
the electron momentum density in a specific material possesses Fermi surface
breaks. Through this technique one may provide guidance in searching for a
suitable theoretical description of the electronic structure, and thus help in
developing an understanding of the mechanism responsible for the occurrence of
specific effects, among them also topological features. Alternatively, the
accuracy of any specific electronic structure model can in principle be
assessed by comparing the computed 2D-ACAR spectra with the corresponding
measurements. On the theoretical side, such comparisons have been already made
using the most advanced LDA+DMFT method for simple transition metals [26],
however its extension to correlated multi-orbital $4f$-systems with strong
spin-orbit coupling is not currently available.

In the present study we performed 2D-ACAR measurements on single crystalline
$\mathrm{LaB}_{6}$ along three high-symmetry directions. The experimental
spectra in both $p$\- and $k$-space were compared to LDA calculations and
showed good agreement. In particular the quantitative results extracted from
the experiment agree excellently with the corresponding theoretical values of
the FS along ${X-M}$, and ${X-R}$ directions; the higher experimental value
found for ${X-\Gamma}$, however, is attributed to the fact that this direction
is experimentally not directly accessible.

For another member of the rare-earth hexaboride class, the $\mathrm{CeB}_{6}$
compound with $4f$-states in the valence band, in absence of the recent
experimental measurements, the LDA+$U$ method has been applied with different
values of the Hubbard parameter $U$. The present study was motivated following
the proposal that $\mathrm{SmB}_{6}$, another member of the hexaboride family,
may be an interaction driven topological insulator [7, 41]. The evidence for a
metallic surface state in this material is largely accepted, and in the same
time low temperature transport [42, 43] and (spin-polarized) angle-resolved
photo-emission spectroscopy measurements [41] point towards a bulk small-gap
insulating state. Quantum oscillation measurements [36] found dHvA
frequencies, cyclotron masses, and amplitudes that are quite similar to other
metallic hexaborides such as $\mathrm{LaB}_{6}$ and $\mathrm{CeB}_{6}$ [44,
45, 46]. These two compounds exhibit a metallic ground state, involving
predominantly conduction electrons with low residual resistivity, and are
characterized according to the DFT (LDA) calculations by a multiply connected
Fermi surface of distorted ellipsoids. Band structure calculations for these
metallic compounds showed that this type of “ellipsoids” are universal Fermi
surface features for these materials. In addition, the same type of
calculations revealed similar features for $\mathrm{SmB}_{6}$, when the Fermi
level is shifted by hand from the computed position, which is in the gap,
either up in the conduction bands or down in the valence bands [36].

Based on the presented results and the combined theoretical and experimental
(for $\mathrm{LaB}_{6}$) analysis, we conclude that the “ellipsoid” cross-
sections ($\alpha$-orbit) and neck sizes increase, hence $\gamma$\- and
$\varepsilon$-hole orbit of $\mathrm{LaB}_{6}$ and $\mathrm{CeB}_{6}$ are
significantly reduced in size in comparison with $\mathrm{SmB}_{6}$. This is
in agreement with the quantum oscillation dHvA measurements, and support the
idea that these two compounds are topologically trivial but correlated metals.
Therefore, as a complementary approach to the present LDA(+$U$) treatment of
$\mathrm{LaB}_{6}$ and $\mathrm{CeB}_{6}$, we might consider the dynamical
mean field theory which is able to take into account the full local
correlation effects and its extensions also including inter-site correlations.
Here the role of the crystal electric-field becomes apparent. The sixfold
degenerate ${j=5/2}$ multiplet state of $\mathrm{Ce}$-$4f^{1}$ is split into a
$\Gamma_{8}$ quartet and the $\Gamma_{7}$ doublet. An effective Hamiltonian
formulation restricted to the $\Gamma_{8}$ quartet only, represents a
significant reduction of the local Hilbert space, and therefore makes the
many-body DMFT computation feasible. This will finally allow the description
of a truly paramagnetic metallic ground state and may reveal more subtle
effects connected to the shape of the Fermi surface.

Acknowledgements

Financial support by the Deutsche Forschungsgemeinschaft through TRR80
(project E2) Project number 107745057 is gratefully acknowledged. We
acknowledge the fruitful long-term collaboration with Stephen B. Dugdale. We
are very grateful to N. Shitsevalova, from Frantsevich Institute for Problems
of Materials Science, National Academy of Sciences of Ukraine, for providing
the $\mathrm{LaB}_{6}$ single crystal.

## References

  * [1] L. Degiorgi, E. Felder, H. R. Ott, J. L. Sarrao, Z. Fisk,  _Phys. Rev. Lett._ 1997, _79_ 5134\.
  * [2] P. S. Riseborough,  _Advances in Physics_ 2000, _49_ , 3 257.
  * [3] V. I. Anisimov, F. Aryasetiawan, A. I. Lichtenstein,  _Journal of Physics: Condensed Matter_ 1997, _9_ , 4 767.
  * [4] F. Aryasetiawan, O. Gunnarsson,  _Reports on Progress in Physics_ 1998, _61_ , 3 237\.
  * [5] G. Kotliar, D. Vollhardt,  _Physics Today_ 2004, _57_ , 3 53.
  * [6] G. Kotliar, S. Y. Savrasov, K. Haule, V. S. Oudovenko, O. Parcollet, C. A. Marianetti,  _Rev. Mod. Phys._ 2006, _78_ , 3 865.
  * [7] M. Dzero, K. Sun, V. Galitski, P. Coleman,  _Phys. Rev. Lett._ 2010, _104_ 106408\.
  * [8] T.-R. Chang, T. Das, P.-J. Chen, M. Neupane, S.-Y. Xu, M. Z. Hasan, H. Lin, H.-T. Jeng, A. Bansil,  _Phys. Rev. B_ 2015, _91_ 155151\.
  * [9] S.-H. Hung, H.-T. Jeng,  _Materials_ 2020, _13_ , 19.
  * [10] P. Hohenberg, W. Kohn,  _Phys. Rev._ 1964, _136_ B864.
  * [11] W. Kohn,  _Rev. Mod. Phys._ 1999, _71_ 1253\.
  * [12] R. O. Jones, O. Gunnarsson,  _Rev. Mod. Phys._ 1989, _61_ 689\.
  * [13] R. O. Jones,  _Rev. Mod. Phys._ 2015, _87_ 897\.
  * [14] K. Tanaka, Y. Ōnuki,  _Acta Crystallographica Section B_ 2002, _58_ , 3 Part 2 423.
  * [15] C.-H. Chen, T. Aizawa, N. Iyi, A. Sato, S. Otani,  _Journal of Alloys and Compounds_ 2004, _366_ , 1 L6.
  * [16] W. Joss, J. M. van Ruitenbeek, G. W. Crabtree, J. L. Tholence, A. P. J. van Deursen, Z. Fisk,  _Phys. Rev. Lett._ 1987, _59_ 1609\.
  * [17] M. Neupane, N. Alidoust, I. Belopolski, G. Bian, S.-Y. Xu, D.-J. Kim, P. P. Shibayev, D. S. Sanchez, H. Zheng, T.-R. Chang, H.-T. Jeng, P. S. Riseborough, H. Lin, A. Bansil, T. Durakiewicz, Z. Fisk, M. Z. Hasan,  _Phys. Rev. B_ 2015, _92_ 104420\.
  * [18] M. Biasini, M. Alam, S. Dugdale, H. Fretwell, Y. Kubo, H. Harima, N. Sato,  _Applied Surface Science_ 1997, _116_ 335, proceedings of the Seventh International Workshop on Slow-Positron Beam Techniques for Solids and Surfaces.
  * [19] M. Biasini, H. M. Fretwell, S. B. Dugdale, M. A. Alam, Y. Kubo, H. Harima, N. Sato,  _Phys. Rev. B_ 1997, _56_ 10192\.
  * [20] M. Biasini, M. Monge, G. Kontrym-Sznajd, M. Gemmi, N. Sato,  In _Positron Annihilation - ICPA-12_ , volume 363 of _Materials Science Forum_. Trans Tech Publications Ltd, 2001 582–584.
  * [21] J. A. Weber, A. Bauer, P. Böni, H. Ceeh, S. B. Dugdale, D. Ernsting, W. Kreuzpaintner, M. Leitner, C. Pfleiderer, C. Hugenschmidt,  _Phys. Rev. Lett._ 2015, _115_ 206404\.
  * [22] H. A. Ceeh, J. A. Weber, P. Böni, M. Leitner, D. Benea, L. Chioncel, H. Ebert, J. Minár, D. Vollhardt, C. Hugenschmidt,  _Scientific Reports_ 2016, _6_ 20898\.
  * [23] T. Jarlborg, A. K. Singh,  _Physical Review B_ 1987, _36_ 4660\.
  * [24] D. G. Lock, V. H. C. Crisp, R. N. West,  _Journal of Physics F: Metal Physics_ 1973, _3_ , 3 561\.
  * [25] M. Leitner, H. Ceeh, J.-A. Weber,  _New Journal of Physics_ 2012, _14_ , 12 123014.
  * [26] H. Ceeh, J. A. Weber, M. Leitner, P. Böni, C. Hugenschmidt,  _Rev. Sci. Instrum._ 2013, _84_ , 4 043905.
  * [27] J. K. Dewhurst, S. Sharma, L. Nordström, F. Cricchio, O. Granas, E. K. U. Gross,  ELK - electronic structure code - version 5.2.14,  URL http://elk.sourceforge.net.
  * [28] J. P. Perdew, Y. Wang,  _Phys. Rev. B_ 1992, _45_ 13244\.
  * [29] A. G. Petukhov, I. I. Mazin, L. Chioncel, A. I. Lichtenstein,  _Phys. Rev. B_ 2003, _67_ , 15 153106.
  * [30] E. Boroński, R. M. Nieminen,  _Phys. Rev. B_ 1986, _34_ , 6 3820.
  * [31] M. J. Puska, R. M. Nieminen,  _Rev. Mod. Phys._ 1994, _66_ 841\.
  * [32] D. Ernsting, D. Billington, T. Haynes, T. Millichamp, J. Taylor, J. Duffy, S. Giblin, J. Dewhurst, S. Dugdale,  _Journal of Physics: Condensed Matter_ 2014, _26_ , 49 495501.
  * [33] J. Laverock, T. D. Haynes, M. A. Alam, S. B. Dugdale,  _Phys. Rev. B_ 2010, _82_ 125127\.
  * [34] N. D. Drummond, P. López Ríos, C. J. Pickard, R. J. Needs,  _Phys. Rev. B_ 2010, _82_ 035107\.
  * [35] N. D. Drummond, P. López Ríos, R. J. Needs, C. J. Pickard,  _Phys. Rev. Lett._ 2011, _107_ 207402\.
  * [36] B. S. Tan, Y.-T. Hsu, B. Zeng, M. C. Hatnean, N. Harrison, Z. Zhu, M. Hartstein, M. Kiourlappou, A. Srivastava, M. D. Johannes, T. P. Murphy, J.-H. Park, L. Balicas, G. G. Lonzarich, G. Balakrishnan, S. E. Sebastian,  _Science_ 2015, _349_ , 6245 287.
  * [37] M. Leitner, J. A. Weber, H. Ceeh,  _New Journal of Physics_ 2016, _18_ , 6 063033.
  * [38] J. P. Perdew, A. Zunger,  _Phys. Rev. B_ 1981, _23_ 5048\.
  * [39] M. T. Czyzyk, G. A. Sawatzky,  _Phys. Rev. B_ 1994, _49_ , 20 14211.
  * [40] A. Koitzsch, N. Heming, M. Knupfer, B. Büchner, P. Y. Portnichenko, A. V. Dukhnenko, N. Y. Shitsevalova, V. B. Filipov, L. L. Lev, V. N. Strocov, J. Ollivier, D. S. Inosov,  _Nature Communications_ 2016, _7_ , 1 10876.
  * [41] M. Dzero, J. Xia, V. Galitski, P. Coleman,  _Annual Review of Condensed Matter Physics_ 2016, _7_ , 1 249.
  * [42] J. W. Allen, B. Batlogg, P. Wachter,  _Phys. Rev. B_ 1979, _20_ 4807\.
  * [43] X. Zhang, N. P. Butch, P. Syers, S. Ziemak, R. L. Greene, J. Paglione,  _Phys. Rev. X_ 2013, _3_ 011011\.
  * [44] Y. Ishizawa, T. Tanaka, E. Bannai, S. Kawai,  _Journal of the Physical Society of Japan_ 1977, _42_ , 1 112.
  * [45] H. Harima, O. Sakai, T. Kasuya, A. Yanase,  _Solid State Communications_ 1988, _66_ , 6 603.
  * [46] Y. Ōnuki, T. Komatsubara, P. H. P. Reinders, M. Springford,  _Journal of the Physical Society of Japan_ 1989, _58_ , 10 3698.