Zero Modes of Massive Fermions
Delocalize from Axion Strings

Hengameh Bagherian, Katherine Fraser, Samuel Homiller, and John Stout

Department of Physics, Harvard University, Cambridge, MA 02138, USA

Abstract
Massless chiral excitations can arise from the interactions between a fermion
and an axion string, propagating along the string and allowing it to
superconduct. The properties of these excitations, or zero modes, dictate how
the string interacts with light and can thus have important phenomenological
consequences. In this paper, we add a nowhere-vanishing Dirac mass for the
fermion in the usual model of axion electrodynamics. We find that the zero
modes exhibit an interesting phase structure in which they delocalize from the
string’s core as the mass increases, up until a critical value past which they
disappear. We study this structure from an analytic perspective, with explicit
numerical solutions, and via anomaly inflow arguments. Finally, we derive the
two-dimensional effective theory of the zero mode and its interactions with
the four-dimensional gauge field and show how this effective theory breaks
down as the zero modes delocalize.

###### Contents

  1. 1 Introduction
  2. 2 Axion Strings
    1. 2.1 The Axion String
    2. 2.2 Fermionic Zero Modes
  3. 3 Adding a Mass
    1. 3.1 The Small Mass Limit
    2. 3.2 Zero Modes on Domain Walls
    3. 3.3 The Critical Mass Case
  4. 4 Zero Modes from Anomaly Inflow
  5. 5 Low-Energy Effective Theory
  6. 6 Conclusions
  7. A Numerical Techniques
  8. B Multipole Moments

##  1 Introduction

Axions—scalar fields with compact field spaces—were originally proposed to
solve the strong CP problem [1, 2, 3, 4], but have since become ubiquitous in
high-energy physics. They have been postulated to serve many roles in beyond
the Standard Model physics: to be dark matter [5, 6, 7]; to be the inflaton
that generates exponential expansion in the early universe [8]; to give rise
to matter/antimatter asymmetry [9]; and many others. Axions are also a generic
prediction in string theory, where they arise copiously as the zero modes of
higher-form fields as a result of the non-trivial topology of the compactified
extra dimensions [10, 11, 12, 13, 14].

One particularly interesting aspect of axion phenomenology is their connection
to topological defects such as monopoles and strings. These connections are
interesting theoretically, as they can provide mechanisms for generating the
axion mass in Abelian gauge theories [15, 16]. They can also be important
phenomenologically: in post-inflationary axion models, for instance, the axion
abundance is determined by radiation from a network of cosmic axion strings
[17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33]. Other
recent work on potential signatures of axion string networks can be found in
[34, 35].

Recently, there has been renewed interest in the fact that axion strings can
be superconducting: they support charged zero modes localized to the string
core, which lead to a current on the string proportional to an applied
electric field. That axion strings are superconducting has been understood for
many years, as it is intimately related to the phenomenon of anomaly inflow
elucidated by Callan and Harvey [36], which we will review in what follows. A
great deal of the physics of anomaly inflow and axion strings was worked out
in subsequent years: see [37, 38, 39, 40, 41, 42, 43, 44, 45] for a collection
of early and important works on the subject. The contemporary interest in
axion string superconductivity is in part due to the realization that these
strings could interact with a primordial magnetic field, leading to striking
signatures due to the formation of vortons [46, 47] or bound states [48].

The existence of the charged zero modes is often explained in the context of
simple models with a classical PQ symmetry. In these models, the fermions
acquire mass from the vacuum expectation value of the radial mode of the
scalar field whose phase is the axion. In classical string configurations, the
vacuum expectation value goes to zero at the core of the string solution and
intuition strongly suggests that massless modes will exist in a region
localized at the string core. An analogous argument successfully explains the
existence of zero modes on domain walls, as we will review.

While this simple picture is intuitive, for axion strings, it is clearly
incomplete. For one, it pays no heed to the crucial fact that the zero modes
on axion strings are chiral—a property which distinguishes them from zero
modes on e.g., Witten strings [49]. Moreover, in more complicated models (for
instance, the DFSZ axion [50, 51]), there can exist string configurations in
which the vacuum expectation value of the scalar does not vanish at the string
core [52]. Despite this, there are arguments that these strings can be
superconducting as well.

This picture is also related to another puzzle. What happens to the zero mode
when this classical PQ symmetry is badly and explicitly broken? For instance,
if the fermions have a very large mass $m$, we should have no trouble
completely integrating them out and the axion string should be blind to their
existence. In this case, there is no anomaly to inflow, and so these zero
modes should not exist in the limit $m\to\infty$. In contrast, since they are
chiral we still expect them to exist for small but non-zero $m$. So, there
must be some critical value of the mass $m$ at which they cease to exist. What
happens to the zero modes near this critical mass?

The goal of this paper is to shed some light on these puzzles. We will do so
by studying a simple model of axion electrodynamics in which the classical PQ
symmetry is explicitly broken by a Dirac mass $m$ for the fermion. We
demonstrate explicitly that there exist zero mode solutions to the equations
of motion in the axion string background, and numerically solve for their
profile. We find that when the Dirac mass is roughly equal in size to the mass
$\mu$ induced by the scalar field which spontaneously breaks the PQ symmetry,
the profile functions change dramatically: the zero mode becomes completely
delocalized from the string and onto a semi-infinite wedge. We also revisit
the anomaly inflow story in the presence of the bulk mass term, clarifying the
topological origin of the zero modes even in the absence of a “classical” PQ
symmetry. Finally, we derive the low-energy two-dimensional effective theory
for the zero mode and calculate the leading higher-derivative interactions
with the bulk gauge field. As one would expect, we show that this effective
theory completely breaks down as $m\to\mu$. While this simple model is not
relevant phenomenologically, our hope is that this work can be applied to more
realistic models, with potential astrophysical or cosmological effects that
can be studied in future work.

It is worth emphasizing that these considerations are entirely distinct from
situations in which the zero modes localized on cosmic (non-axionic) strings
can acquire mass from the pairing of left- and right-moving modes [53], which
can be induced e.g., by finite temperature effects [54]. The mass term we
introduce is for the “full” theory of fermions propagating in four-dimensional
spacetime. The zero modes localized to the string remain massless on
topological grounds, as in the original example [36]. The zero modes we
discuss are also distinct from those found in the background of “$Z$-strings”
[55, 56, 57, 58]. These string configurations are not topologically stable and
there is no mechanism that protects the zero modes from pairing up and
acquiring a mass [59, 60].

Outline The rest of this paper is structured as follows. In §2, we review the
axion string solution in the usual case, with a global PQ symmetry, and solve
for the massless zero modes explicitly. In §3, we break the PQ symmetry with a
Dirac mass and demonstrate that the zero modes still exist. We discuss how
this can be understood by analogy to the existence of zero modes on domain
walls in $2+1$ dimensions, present numerical results for the profile of these
zero modes on the string, and discuss the behavior in the “critical” mass
case. In §4, we recap the original anomaly inflow story, with appropriate
modifications to account for the mass term. Finally, in §5, we derive the low-
energy theory of the zero modes on the string, and discuss how it is impacted
by the Dirac mass term. We conclude in §6. Two Appendices provide more details
on the numerical techniques used to solve for the zero modes using Chebyshev
interpolation, and on the derivation of the low-energy effective action.

##  2 Axion Strings

In this section, we will review the original ultraviolet completion of an
axion string studied by Callan and Harvey [36], with Lagrangian

$\mathcal{L}=-\frac{1}{4}F_{\mu\nu}F^{\mu\nu}+\bar{\psi}i\not{D}\psi+|\partial_{\mu}\Phi|^{2}+y\hskip
0.5pt\bar{\psi}\\!\hskip
0.5pt\left(\Phi_{1}+i\gamma^{5}\Phi_{2}\right)\\!\hskip 0.5pt\psi-V(\Phi)\,.$
(2.1)

This theory lives in four-dimensional Minkowski space $\mathcal{M}_{4}$ and
consists of an abelian gauge field $A_{\mu}$ with field strength
$F_{\mu\nu}=\partial_{\mu}A_{\nu}-\partial_{\nu}A_{\mu}$, an uncharged complex
scalar field $\Phi\equiv\Phi_{1}+i\Phi_{2}\equiv
f(x){{\mathrm{e}}}^{i\theta(x)}$ with potential

$V(\Phi)=\lambda\left(|\Phi|^{2}-v^{2}\right)^{2}\,,$ (2.2)

and a single charged Dirac fermion $\psi$, which chirally couples to $\Phi$
with strength $y$. We use $D_{\mu}=\partial_{\mu}-ieA_{\mu}$ to denote the
gauge covariant derivative, where $\not{D}\equiv\gamma^{\mu}D_{\mu}$, and
$\gamma^{\mu}$ are the standard Dirac gamma matrices with
$\\{\gamma^{\mu},\gamma^{\nu}\\}=2\eta^{\mu\nu}$ and $\gamma^{5}\equiv
i\gamma^{0}\gamma^{1}\gamma^{2}\gamma^{3}$. This theory enjoys a
${\mathrm{U}}(1)_{\textsc{pq}}$ global symmetry, commonly called the Peccei-
Quinn symmetry, under which the complex scalar and fermion transforms as
$\Phi\to{{\mathrm{e}}}^{i\alpha}\Phi$ and
$\psi\to{{\mathrm{e}}}^{{\scalebox{0.5}[0.85]{$-$}}i\gamma^{5}\alpha/2}\psi$,
respectively.

###  2.1 The Axion String

The potential (2.2) forces the scalar $\Phi$ to acquire a vacuum expectation
value $\langle\Phi\rangle=v$. It also allows for the existence of static
axionic string configurations,

$\Phi_{n}(x)=f(r)\hskip 0.5pt{{\mathrm{e}}}^{in\varphi}\,,$ (2.3)

characterized by an integer topological charge $n\in\mathbb{Z}$. We work in
standard cylindrical coordinates $(x,y,z)=(r\cos\varphi,r\sin\varphi,z)$, with
the string oriented along the $z$-axis. In these configurations, the phase of
the complex scalar field, i.e. the “axion”
$\theta(x)\equiv{\mathrm{arg}}\,\Phi(x)$, winds $n$ times around its field
space $\theta\sim\theta+2\pi$ as we move around the string. Said differently,

$n\equiv\frac{1}{2\pi}\oint_{\gamma}\mathrm{d}\theta\,,$ (2.4)

where $\gamma$ is a closed contour that encircles the string at $r=0$. The
“radial profile” of this string is determined by the real function $f(r)$. Far
from the string, the potential (2.2) forces the scalar to sit in the minimum
of its potential and so $f(r)\to v$ as $r\to\infty$. Likewise, regularity of
the solution forces $f(r)$ to vanish in the core of the string: $f(r)\to 0$ as
$r\to 0$. These solutions are also called “global vortices” in the literature,
since they are vortex solutions that are charged under the _global_
${\mathrm{U}}(1)_{\textsc{pq}}$ symmetry, as opposed to a gauge symmetry.

The tension of an axion string with charge $n$ is given by

$T_{n}=2\pi\int_{0}^{L}\\!\mathrm{d}r\,r\left[\left(\frac{\mathrm{d}f}{\mathrm{d}r}\right)^{2}\\!+\frac{n^{2}}{r^{2}}f^{2}+\lambda(f^{2}-v^{2})^{2}\right]\,.$
(2.5)

Since $f(r)\to v$ as $r\to\infty$, the tension diverges logarithmically with
the size of the system or IR cutoff $L$,

$T_{n}\sim 2\pi
v^{2}n^{2}\log\left(\frac{L}{r_{\mathrm{core}}}\right)\,,\mathrlap{\qquad
L\to\infty\,.}$ (2.6)

Here, we have introduced the “size” of the string $r_{\mathrm{core}}$, which
is of order
$r_{\mathrm{core}}=\big{(}\sqrt{\lambda}v\big{)}^{{\scalebox{0.5}[0.85]{$-$}}1}$.
This IR divergent total energy is common for global strings and is typically
regulated by either placing the theory in a box or by positing that there is
another string of opposite charge a distance $L$ away, as is often the case
in, e.g. cosmological simulations of axion string networks.

Minimizing the tension yields the equation of motion

$\frac{1}{r}\frac{\mathrm{d}}{\mathrm{d}r}\\!\left(r\hskip
0.5pt\frac{\mathrm{d}f}{\mathrm{d}r}\right)-\left[\frac{n^{2}}{r^{2}}+2\lambda(f^{2}-v^{2})\right]f=0\,.$
(2.7)

In general, this equation cannot be solved analytically but is amenable to
numerics. We will restrict our attention throughout to strings of charge
$n=+1$. An asymptotic analysis of (2.7) shows that the radial profile $f(r)$
behaves as

$f(r)\sim\begin{dcases}v\hskip 0.5pt\hskip 0.5ptC_{1}\hskip
0.5pt\big{(}\sqrt{\lambda}vr\big{)}+\cdots&r\to 0\\\
v\\!\left(1-\frac{1}{4\lambda
v^{2}r^{2}}\right)+\cdots&r\to\infty\end{dcases}\,,$ (2.8)

with $C_{1}$ an overall constant that can be analytically determined by
matching the two solutions in an intermediate region, or by numerics. It is
clear from the behavior as $r\to\infty$ that the string solution varies over
length scales
$r_{\mathrm{core}}=\big{(}\sqrt{\lambda}v\big{)}^{{\scalebox{0.5}[0.85]{$-$}}1}$.
We may then use the techniques outlined in Appendix A to numerically search
for the solution which varies over scales set by $r_{\mathrm{core}}$ with the
correct asymptotics (2.8). We show this profile in Figure 1.

Figure 1: The radial profile $f(r)$ [ blue], in units of $v$, for the axion
string of charge $\pm 1$ and corresponding (non-normalized) zero mode profile
$F(r)$ [ orange] with $y=\sqrt{\lambda}$.

###  2.2 Fermionic Zero Modes

We are interested in the low energy dynamics of (2.1) in the background of the
axion string. Absent the string, the fermion acquires a Dirac mass
$\mu=yv$—which we call the “Yukawa mass” since it arises from the spontaneous
breaking of the ${\mathrm{U}}(1)_{\textsc{pq}}$ symmetry—and so the low-energy
theory only contains the massless axion and the abelian gauge field. One might
think that the low-energy theory in the presence of the axionic string is
similar; far from the string, the fermions again get a mass $\mu$ and the low-
energy theory should be unchanged. However, as argued in [36] this reasoning
is incomplete: there are instead massless chiral zero modes which are
localized to the string. In this section, we will review the direct
construction of these zero modes. A common intuitive explanation is that these
massless modes exist because the radially-dependent “mass” of the bulk
fermions $yf(\rho)$ vanishes in the core of the string. However, we argue in
Section 3 that this is not a necessary condition, while in the Section 4 we
discuss why these modes are necessary from the perspective of anomaly inflow.

In the presence of the string of charge $n=+1$ and no background gauge field,
the fermions obey the equation of motion

$\left(i\not{\partial}+y\hskip 0.5ptf(r)\hskip
0.5pt{{\mathrm{e}}}^{i\gamma^{5}\varphi}\right)\\!\psi=0\,.$ (2.9)

For simplicity, we will restrict to solutions of the form

$\psi(x)=\psi(r,\varphi)\,{{\mathrm{e}}}^{{\scalebox{0.5}[0.85]{$-$}}ip(t+z)}\,,$
(2.10)

where $\psi(r,\varphi)=\psi_{\alpha}(r,\varphi)$ is a four-component spinor
function with $\alpha=0,\ldots,3$. This ansatz describes solutions that travel
in the $({\scalebox{0.75}[1.0]{$-$}}z)$-direction at the speed of light. In
the Weyl representation, (2.9) reduces to the set of coupled equations

$\begin{aligned} 0&=2p\hskip
0.5pt\psi_{2}+{{\mathrm{e}}}^{{\scalebox{0.5}[0.85]{$-$}}i\varphi}\left(yf(r)\hskip
0.5pt\psi_{0}+i\partial_{r}\psi_{3}+r^{{\scalebox{0.5}[0.85]{$-$}}1}\partial_{\varphi}\psi_{3}\right)\\\
0&=yf(r){{\mathrm{e}}}^{{\scalebox{0.5}[0.85]{$-$}}i\varphi}\hskip
0.5pt\psi_{1}+{{\mathrm{e}}}^{i\varphi}\left(i\partial_{r}\psi_{2}-r^{{\scalebox{0.5}[0.85]{$-$}}1}\partial_{\varphi}\psi_{2}\right)\\\
0&=yf(r){{\mathrm{e}}}^{i\varphi}\hskip
0.5pt\psi_{2}-{{\mathrm{e}}}^{{\scalebox{0.5}[0.85]{$-$}}i\varphi}\left(i\partial_{r}\psi_{1}+r^{{\scalebox{0.5}[0.85]{$-$}}1}\partial_{\varphi}\psi_{1}\right)\\\
0&=2p\hskip 0.5pt\psi_{1}+{{\mathrm{e}}}^{i\varphi}\left(yf(r)\hskip
0.5pt\psi_{3}-i\partial_{r}\psi_{0}+r^{{\scalebox{0.5}[0.85]{$-$}}1}\partial_{\varphi}\psi_{0}\right)\end{aligned}\,.$
(2.11)

Such solutions should have definite helicity, and so we may set
$\psi_{1}=\psi_{2}=0$. We then see that $\psi_{0}=\psi_{0}(r)$ and
$\psi_{3}=\psi_{3}(r)$ are purely radial functions that satisfy the coupled
system of equations

$\begin{aligned} 0&=yf(r)\hskip 0.5pt\psi_{0}+i\partial_{r}\psi_{3}\\\
0&=yf(r)\hskip 0.5pt\psi_{3}-i\partial_{r}\psi_{0}\end{aligned}\,.$ (2.12)

These equations can be easily solved to find two solutions of the form

$\psi_{p}(x)=\sqrt{\frac{p}{2}}\begin{pmatrix}1\\\ 0\\\ 0\\\ \mp
i\end{pmatrix}{{\mathrm{e}}}^{{\scalebox{0.5}[0.85]{$-$}}ip(t+z)}F(r)\,,$
(2.13)

where $F^{\prime}(r)=\pm yf(r)F(r)$ or

$F(r)=\mathcal{A}\exp\left(\pm
y\int_{0}^{r}\\!\mathrm{d}r^{\prime}\,f(r^{\prime})\right),$ (2.14)

with $\mathcal{A}$ an overall constant of integration. For this solution to be
appropriately normalizable,

$\int_{\mathbb{R}^{3}}\\!\mathrm{d}^{3}x\,\psi^{\dagger}_{p}(x)\psi_{p^{\prime}}(x)=2\pi|p|\hskip
0.5pt\hskip 0.5pt\delta(p-p^{\prime})\,,$ (2.15)

or specifically

$2\pi\int_{0}^{\infty}\\!\mathrm{d}r\,r\hskip 0.5pt\hskip 0.5ptF^{2}(r)=1\,,$
(2.16)

so we must select the negative sign in (2.14). We thus find that the axion
string with charge $+1$ supports a chiral fermionic zero mode that travels at
the speed of light in the $({\scalebox{0.75}[1.0]{$-$}}z)$-direction, with
form

$\psi(x)=\mathcal{A}\sqrt{\frac{p}{2}}\begin{pmatrix}1\\\ 0\\\ 0\\\
{\scalebox{0.5}[0.85]{$-$}}i\end{pmatrix}{{\mathrm{e}}}^{{\scalebox{0.5}[0.85]{$-$}}ip(t+z)}\exp\left({\scalebox{0.75}[1.0]{$-$}}y\int_{0}^{r}\\!\mathrm{d}r^{\prime}\,f(r^{\prime})\right)\,.$
(2.17)

If we were to repeat this analysis with the charge
${\scalebox{0.75}[1.0]{$-$}}1$ axion string, then it would support a
normalizable chiral zero mode that instead travels at the speed of light in
the $(+z)$-direction,

$\psi(x)=\mathcal{A}\sqrt{\frac{p}{2}}\begin{pmatrix}0\\\ 1\\\
{\scalebox{0.5}[0.85]{$-$}}i\\\
0\end{pmatrix}{{\mathrm{e}}}^{{\scalebox{0.5}[0.85]{$-$}}ip(t-z)}\exp\left({\scalebox{0.75}[1.0]{$-$}}y\int_{0}^{r}\\!\mathrm{d}r^{\prime}\,f(r^{\prime})\right)\,.$
(2.18)

In Figure 1, we show the radial profile $F(r)$ for $y=\sqrt{\lambda}$. These
zero modes are localized in a region about the string of size
$r_{\textsc{zm}}=r_{\mathrm{core}}/y$.

Before we move on, it will be useful to understand why the presence of the
axion string permits these normalizable modes to exist, solely from the
perspective of the differential equations. We will see that, contrary to naive
expectation, it is not because the fermion is “massless” in the core of the
string. With (2.10), we can decouple the equations in (2.12) to find two
copies of the second order equation

$\left[\frac{\mathrm{d}^{2}}{\mathrm{d}r^{2}}-\frac{f^{\prime}(r)}{f(r)}\frac{\mathrm{d}}{\mathrm{d}r}-y^{2}f^{2}(r)\right]\\!F(r)=0\,,$
(2.19)

with $\psi_{0}(r)=F(r)$ and $\psi_{3}(r)={\scalebox{0.75}[1.0]{$-$}}iF(r)$. As
is clear from Figure 1, the function $f(r)$ is smooth and non-singular for all
$r\in(0,\infty)$. This implies that the solutions to (2.19) are regular
everywhere, except possibly as $r\to 0$ and $r\to\infty$, which are regular
and irregular singular points, respectively. A normalizable solution must be
regular at both of these points and generally fails to exist because the
solution that is regular about one singularity is not regular about the other.

Using the asymptotics (2.8) of the background $f(r)$, far from the string
(2.19) reduces to

$\left[\frac{\mathrm{d}^{2}}{\mathrm{d}r^{2}}-\mu^{2}+\cdots\right]\\!F(r)=0\,,$
(2.20)

which has solutions that behave as

$F(r)=C_{\scriptscriptstyle-}F_{\scriptscriptstyle-}(r)+C_{\scriptscriptstyle+}F_{\scriptscriptstyle+}(r)\sim
C_{\scriptscriptstyle-}\hskip
0.5pt{{\mathrm{e}}}^{{\scalebox{0.5}[0.85]{$-$}}\mu
r}+C_{\scriptscriptstyle+}\hskip 0.5pt{{\mathrm{e}}}^{\mu
r}\,,\mathrlap{\qquad r\to\infty\,,}$ (2.21)

where we have defined two linearly-independent solutions $F_{\pm}(r)$ with
definite scaling as $r\to\infty$. In this limit, the fermion zero modes do not
“see” the string, but merely feel their acquired mass $\mu=yv$. However, near
the core of the string, (2.19) instead reduces to

$\left[\frac{\mathrm{d}^{2}}{\mathrm{d}r^{2}}-\frac{1}{r}\frac{\mathrm{d}}{\mathrm{d}r}+\cdots\right]\\!F(r)=0\,.$
(2.22)

As we move towards the core of string $r\to 0$, there is a balance between the
first and second terms in (2.19), the latter of which only appears due to the
existence of the string. The third term, due to the Yukawa interaction, can be
neglected. Solutions then behave as

$F(r)=C_{1}F_{1}(r)+C_{2}F_{2}(r)\sim C_{1}+C_{2}r^{2}\,,\mathrlap{\qquad r\to
0\,,}$ (2.23)

where we have again identified two solutions $F_{1,2}(r)$ which have definite
scaling as $r\to 0$. Importantly, both of these solutions are regular near the
origin.

We can understand why the axion string allows a normalizable zero mode to
exist as follows. Since (2.19) is a second-order differential equation, it has
only two linearly-independent solutions. Above, we found that a general
solution can either be expressed as a linear combination of $F_{1,2}(r)$ or
$F_{\pm}(r)$. Since these two sets of solutions are not linearly independent,
there always exists a linear map between them. For a solution to be
normalizable, it must decay as $r\to\infty$ and so we must have that
$C_{\scriptscriptstyle+}=0$; therefore we must be able to write this solution
as
$F(r)=C_{\scriptscriptstyle-}F_{\scriptscriptstyle-}=C_{{\scriptscriptstyle-},1}F_{1}(r)+C_{{\scriptscriptstyle-},2}F_{2}(r)$.
It is often the case that only one of the solutions $F_{1}(r)$ or $F_{2}(r)$
is well-behaved at the origin $r=0$ and, unless there is some special
structure that ensures it, the solution that is well-behaved as $r\to\infty$
will not be well-behaved as $r\to 0$. Fortunately, the presence of the string
ensures that _both_ linearly-independent solutions are regular as $r\to 0$.
Since $f(r)$ is a smooth function, we are also guaranteed that the solutions
are regular for positive $r$. Consequently, a normalizable solution to (2.19)
exists, regardless of its exact form.

An analogous situation occurs for the unbound states of the quantum mechanical
hydrogen atom, whose radial modes obey

$\left[-\frac{1}{2r^{2}}\frac{\mathrm{d}}{\mathrm{d}r}\left(r^{2}\frac{\mathrm{d}}{\mathrm{d}r}\right)+\frac{\ell(\ell+1)}{2r^{2}}-\frac{\alpha}{r}-E\right]\\!\psi(r)=0\,.$
(2.24)

Due to the centrifugal barrier, solutions behave as $\psi(r)\sim
C_{1}r^{\ell}+C_{2}r^{-(\ell+1)}$ as $r\to 0$ and $\psi(r)\sim
C_{\scriptscriptstyle-}{{\mathrm{e}}}^{{\scalebox{0.5}[0.85]{$-$}}\sqrt{{\scalebox{0.5}[0.85]{$-$}}2E}r}+C_{\scriptscriptstyle+}{{\mathrm{e}}}^{\sqrt{{\scalebox{0.5}[0.85]{$-$}}2E}r}$
as $r\to\infty$. Bound states, $E<0$, have a discrete spectrum because it is
impossible to simultaneously impose that the wavefunction be both regular at
the origin _and_ exponentially decay at spatial infinity, except at a discrete
set of energy eigenvalues. For all other energies, solutions with
$C_{\scriptscriptstyle+}=0$ necessarily have $C_{2}\neq 0$. For unbound
states, $E>0$, there is no restriction on the behavior as $r\to\infty$, and so
there exists a continuum of regular solutions with $C_{2}=0$, for all positive
energies.

Having reviewed how axion strings, and their associated fermionic zero modes,
arise in the Callan–Harvey model (2.1), we are now in a position to understand
how these zero modes behave as we deform (2.1). In the next section, we will
study what happens to these zero modes in the presence of a Dirac mass $m$ for
the four-dimensional fermion. Even though the four-dimensional fermion is
everywhere massive, these fermionic zero modes still exist as long as
$|m|<\mu$. Furthermore, we argue that there is an interesting “phase
structure” in which the zero modes become increasingly unbound from the string
as $|m|\to\mu$, eventually disappearing for $|m|>\mu$.

##  3 Adding a Mass

In the previous section, we reviewed how fermionic zero modes could arise
along an axion string in the context of a simple toy model. In that specific
UV completion, the axion string is a solitonic object in which the
${\mathrm{U}}(1)_{\textsc{pq}}$ symmetry that is spontaneously broken far from
the string is restored at the core. This causes the four-dimensional fermion
to see effectively zero mass near the core, and a popular refrain is that this
why we should expect to see massless fermionic zero modes localized on the
string. The goal of this section is to explain why this is not a necessary
condition and see that fermionic zero modes can arise even in theories where
bulk four-dimensional fermions have a nowhere vanishing mass.

We will deform (2.1) by adding a mass $m$ to the Dirac fermion $\psi$,

$\mathcal{L}=-\frac{1}{4}F_{\mu\nu}F^{\mu\nu}+\bar{\psi}(i\not{D}-m)\psi+|\partial_{\mu}\Phi|^{2}+y\hskip
0.5pt\bar{\psi}\left(\Phi_{1}+i\gamma^{5}\Phi_{2}\right)\psi-V(\Phi)\,,$ (3.1)

again with $V(\Phi)=\lambda(|\Phi|^{2}-v^{2})$. We restrict to the case where
the amplitude of $\psi$ is small. This theory permits the same axion string
solution (2.3) as before. However, the fermion now obeys a modified equation,

$\left(i\not{\partial}-M(r,\varphi)\hskip
0.5pt{{\mathrm{e}}}^{i\gamma^{5}\alpha(r,\varphi)}\right)\\!\psi=0\,,$ (3.2)

where $M(r,\varphi)\hskip
0.5pt{{\mathrm{e}}}^{i\gamma^{5}\alpha(r,\varphi)}=m-y\hskip 0.5ptf(r)\hskip
0.5pt\hskip 0.5pt{{\mathrm{e}}}^{i\gamma^{5}\varphi}$ or, explicitly,

$\begin{aligned}
M(r,\varphi)&=\sqrt{\big{(}m-yf(r)\cos\varphi\big{)}^{2}+y^{2}f^{2}(r)\sin^{2}\varphi}\\\
\alpha(r,\varphi)&={\mathrm{arg}}\hskip 0.5pt\hskip 0.5pt\big{(}m-yf(r)\hskip
0.5pt{{\mathrm{e}}}^{{\scalebox{0.5}[0.85]{$-$}}i\varphi}\big{)}=i\log\\!\left(\frac{m-yf(r){{\mathrm{e}}}^{{\scalebox{0.5}[0.85]{$-$}}i\varphi}}{M(r,\varphi)}\right)\end{aligned}\,\,.$
(3.3)

In this case, the fermion mass $M(r,\varphi)$ no longer goes to zero at the
core of the string even though $f(r)\to 0$, but instead approaches the “core
mass” $M(r=0,\varphi)=m$. Even so, we will still find that this string
supports fermionic zero modes as long as this core mass is less than the
Yukawa mass, $|m|<\mu$.

To find these zero modes, we again restrict to axion strings with charge $+1$
and search for zero modes of the form (2.10). The mass deformation
$m\bar{\psi}\psi$ explicitly breaks the chiral symmetry or, analogously, the
continuous axion shift symmetry $\varphi\to\varphi+\alpha$, and thus changes
the form of the axion string solution. However, classically this potential is
controlled by fermion’s number density $\bar{\psi}\gamma^{0}\psi$, and in the
limit of small amplitude—or, quantum mechanically, small occupation
numbers—this effect is subleading and can be ignored. Similarly, we will not
consider the electromagnetic field generated by the charged fermion, setting
the vector potential $A_{\mu}=0$, as this effect is also subleading in the
limit of small amplitudes.

It is worth noting that integrating out the fermion generates a potential for
the scalar that depends on $\varphi$. The leading Coleman–Weinberg estimate
$\propto M^{4}(r,\varphi)\log(M(r,\varphi)/\mu)$ [61] is minimized along
$\varphi=0$ and leads to a domain wall for the axion that emanates off the
string and extends in this direction. In what follows, we are explicitly
assuming that this contribution is fine-tuned away, though we will see that
there is a remnant of the axion domain wall in the orientation of the zero
mode profiles.

Figure 2: The profile functions of the zero mode wavefunction
$|\psi_{0}(r,\varphi)|$ in the presence of a Dirac core mass $m$ and
normalized according to (2.15), for various values of $m/\mu$ and
$\sqrt{\lambda}=y=\sqrt{0.1}$, with contours in units of $\sqrt{\omega}\mu$.
The other component satisfies
$\psi_{3}(r,\varphi)={\scalebox{0.75}[1.0]{$-$}}i\psi_{0}(r,{\scalebox{0.75}[1.0]{$-$}}\varphi)$
and is thus identical in structure. The white dot denotes the origin of the
string, $r=0$. As $m\to\mu$, the zero mode profiles spread out and become
unbound from the string—note the difference in scales of the contours between
the two rows. For $m>\mu$, these zero modes cease to exist. Qualitatively
similar results apply when $\sqrt{\lambda}\neq y$.

By Lorentz symmetry, we can set $\psi_{1}(r,\varphi)=\psi_{2}(r,\varphi)=0$
since the modes we are searching for propagate in the
$({\scalebox{0.75}[1.0]{$-$}}z)$-direction at the speed of light, and so in
the Weyl representation (3.2) reduces to

$\begin{aligned}
0&=-m\psi_{0}+{{\mathrm{e}}}^{{\scalebox{0.5}[0.85]{$-$}}i\varphi}\\!\left(yf(r)\hskip
0.5pt\psi_{0}+i\partial_{r}\psi_{3}+r^{{\scalebox{0.5}[0.85]{$-$}}1}\partial_{\varphi}\psi_{3}\right)\\\
0&=-m\psi_{3}+{{\mathrm{e}}}^{i\varphi}\\!\left(yf(r)\hskip
0.5pt\psi_{3}-i\partial_{r}\psi_{0}+r^{{\scalebox{0.5}[0.85]{$-$}}1}\partial_{\varphi}\psi_{0}\right)\end{aligned}\,.$
(3.4)

We can solve (3.4) numerically using Chebyshev interpolation, whose details we
describe in Appendix A. The resulting profile functions are shown in Figure 2
for various values of $m/\mu$. In the plot, we see the normalized zero mode as
a function of dimensionless $x$ and $y$ variables, with a white dot
corresponding to the location of the origin and hence the axion string.
Surprisingly, we find that this zero mode becomes less localized to the string
the larger the core mass $m$ is. For large $m\lesssim\mu$, it is no longer
centered at the core of the axion string and instead stretches out along the
positive $x$-axis. This direction is determined by the orientation of the
string and specifically the axis along which $\varphi=0$, corresponding to the
minimum of $M(r,\varphi)$. As $m\to\mu$, we find that there is a sort of phase
transition in which case the zero mode completely delocalizes onto a two-
dimensional wedge. Finally, the mode completely disappears for $m>\mu$.

To better understand the counterintuitive results shown in Figure 2, in §3.1
we first solve (3.4) by treating $m$ as a small perturbation. As in §2.2, we
will see that it is the presence of the axion string, and not a vanishing mass
at the core, which allows for this zero mode to exist. In §3.2, we describe a
simple analog of this mechanism for Dirac fermions in the three-dimensional
half-space $\mathbb{R}^{2,1}_{+}$. In section §3.3, we analyze the zero mode
for the critical case $m=\mu$.

###  3.1 The Small Mass Limit

In this section, we will argue that a normalizable zero mode solution to (3.4)
exists in the small Dirac mass limit by analyzing the structure of (3.4) when
$m\ll\mu$. Expanding each spinor into Fourier modes,

$\psi_{\alpha}(r,\varphi)=\sum_{\ell\in\mathbb{Z}}\psi_{\alpha,\ell}(r){{\mathrm{e}}}^{i\ell\varphi}\,,$
(3.5)

the zero mode equation (3.4) reduces to an infinite set of coupled equations

$\begin{aligned} 0&=-m\psi_{0,\ell-1}+yf(r)\hskip
0.5pt\psi_{0,\ell}+i\partial_{r}\psi_{3,\ell}+i\ell
r^{{\scalebox{0.5}[0.85]{$-$}}1}\psi_{3,\ell}\\\
0&=-m\psi_{3,\ell+1}+yf(r)\hskip
0.5pt\psi_{3,\ell}-i\partial_{r}\psi_{0,\ell}+i\ell
r^{{\scalebox{0.5}[0.85]{$-$}}1}\psi_{0,\ell}\,.\end{aligned}\,$ (3.6)

When $m=0$, the length scale over which the zero mode varies is set by the
Yukawa mass $\mu=yv$, which dominates the system of equations (3.6) as
$r\to\infty$. The mass $m$ will also be relevant as $r\to\infty$, and so our
perturbative expansion will be in powers of $m/\mu$.

We will thus search for solutions of the form

$\psi(x)=\sqrt{\frac{p}{2}}\left[\begin{pmatrix}1\\\ 0\\\ 0\\\
{\scalebox{0.5}[0.85]{$-$}}i\end{pmatrix}F(r)+\delta\psi(r,\varphi)\right]{{\mathrm{e}}}^{{\scalebox{0.5}[0.85]{$-$}}ip(t+z)}\,,$
(3.7)

and work to first order in $\delta\psi=\mathcal{O}(m/\mu)$. In this limit,
(3.6) reduces to two identical sets of differential equations,

$\begin{aligned}
\left(\frac{\mathrm{d}}{\mathrm{d}r}-\frac{1}{r}\right)\\!G(r)&=yf(r)H(r)\\\
\left(\frac{\mathrm{d}}{\mathrm{d}r}+\frac{1}{r}\right)\\!H(r)&=yf(r)G(r)-mF(r)\end{aligned}\,,$
(3.8)

where $(G,H)=(\delta\psi_{0,1},{\scalebox{0.75}[1.0]{$-$}}i\delta\psi_{3,1})$
or
$(G,H)=(i\delta\psi_{3,{\scalebox{0.5}[0.85]{$-$}}1},{\scalebox{0.75}[1.0]{$-$}}\delta\psi_{0,{\scalebox{0.5}[0.85]{$-$}}1})$.
All other equations in (3.6) are trivially satisfied by setting the other
$\delta\psi_{\alpha,\ell}(r)=0$. The function $G(r)$ satisfies the
inhomogeneous equation

$\displaystyle\left[\frac{\mathrm{d}^{2}}{\mathrm{d}r^{2}}-\frac{f^{\prime}(r)}{f(r)}\frac{\mathrm{d}}{\mathrm{d}r}+\frac{f^{\prime}(r)}{rf(r)}-y^{2}f^{2}(r)\right]G(r)$
$\displaystyle=-myf(r)F(r)\,,$ (3.9)

which can be solved using solutions to the homogeneous equation (with $m=0$)
via variation of parameters [62]. These homogeneous solutions will determine
how the inhomogeneous solution behaves, and so it will be helpful to
understand their asymptotic behavior.

Far from the string, at $r\to\infty$, the homogeneous equation reduces to

$\left[\frac{\mathrm{d}^{2}}{\mathrm{d}r^{2}}-\mu^{2}+\cdots\right]\\!G(r)=0\,,$
(3.10)

in which case the general solution behaves as

$G(r)=C_{\scriptscriptstyle-}G_{\scriptscriptstyle-}(r)+C_{\scriptscriptstyle+}G_{\scriptscriptstyle+}(r)\sim
C_{\scriptscriptstyle-}{{\mathrm{e}}}^{{\scalebox{0.5}[0.85]{$-$}}\mu
r}+C_{\scriptscriptstyle+}{{\mathrm{e}}}^{\hskip 0.5pt\mu
r}\,,\mathrlap{\qquad r\to\infty\,,}$ (3.11)

where we have again introduced the two linearly-independent solutions
$G_{\pm}(r)$ with definite scaling as $r\to\infty$. We will choose their
normalization such that the Wronskian is

$W(r)=2yf(r)\,.$ (3.12)

Likewise, as $r\to 0$, the homogeneous equation simplifies to

$\left[\frac{\mathrm{d}^{2}}{\mathrm{d}r^{2}}-\frac{1}{r}\frac{\mathrm{d}}{\mathrm{d}r}+\frac{1}{r^{2}}+\cdots\right]\\!G(r)=0\,.$
(3.13)

As before, the presence of the axion string introduces a regular singularity
at $r=0$, and solutions take the form

$G_{\scalebox{0.75}[0.75]{$\pm$}}(r)=C_{{\scalebox{0.75}[0.75]{$\pm$}},1}\sum_{k=0}^{\infty}a^{{\scalebox{0.75}[0.75]{$\pm$}}}_{k}r^{k+1}+C_{{\scalebox{0.75}[0.75]{$\pm$}},2}\sum_{k=0}^{\infty}b^{{\scalebox{0.75}[0.75]{$\pm$}}}_{k}r^{k+1}\log
r\,.$ (3.14)

Crucially, both of the linearly-independent homogeneous solutions are regular
as $r\to 0$, regardless of the overall coefficients
$C_{{\scalebox{0.75}[0.75]{$\pm$}},1}$ or
$C_{{\scalebox{0.75}[0.75]{$\pm$}},2}$. The most general solution to (3.9)
that decays as $r\to\infty$ is

$G(r)=\frac{1}{2}m\hskip
0.5ptC_{\scriptscriptstyle-}G_{\scriptscriptstyle-}(r)+\frac{1}{2}m\,G_{\scriptscriptstyle+}(r)\\!\int^{\infty}_{r}\\!\mathrm{d}r^{\prime}\,F(r^{\prime})\hskip
0.5ptG_{\scriptscriptstyle-}(r^{\prime})+\frac{1}{2}m\,G_{\scriptscriptstyle-}(r)\\!\int^{r}_{0}\\!\mathrm{d}r^{\prime}\,F(r^{\prime})\hskip
0.5ptG_{\scriptscriptstyle+}(r^{\prime})\,,$ (3.15)

for arbitrary values of $C_{\scriptscriptstyle-}$. Since $F(r)$ and
$G_{\scriptscriptstyle-}(r)$ both decay as
${{\mathrm{e}}}^{{\scalebox{0.5}[0.85]{$-$}}\mu r}$, while
$G_{\scriptscriptstyle+}(r^{\prime})$ grows as ${{\mathrm{e}}}^{\hskip
0.5pt\mu r}$, all three terms in (3.15) decay as
${{\mathrm{e}}}^{{\scalebox{0.5}[0.85]{$-$}}\mu r}$ as $r\to\infty$. Likewise,
both $F(r)$ and the homogeneous solutions $G_{\pm}(r)$ are regular as $r\to
0$, and thus so is (3.15). This solution is thus normalizable.

Given a solution for $G(r)$, we can solve for the other profile function
$H(r)$ via

$H(r)=\frac{1}{yf(r)}\left(\frac{\mathrm{d}}{\mathrm{d}r}-\frac{1}{r}\right)\\!G(r)\,.$
(3.16)

This equation should cause some alarm: since $f(r)\propto\sqrt{\lambda}vr$
diverges as $r\to 0$, we might worry that $H(r)$ is _not_ regular as $r\to 0$
even though $G(r)$ is. However, the free coefficient $C_{\scriptscriptstyle-}$
in (3.15) can generally be chosen to yield a normalizable $H(r)$. Writing
(3.16) as

$\begin{aligned}
H(r)&=\frac{m}{2yf(r)}\Bigg{[}C_{\scriptscriptstyle-}\,r\hskip
0.5pt\frac{\mathrm{d}}{\mathrm{d}r}\\!\left(\frac{G_{\scriptscriptstyle-}(r)}{r}\right)+r\hskip
0.5pt\frac{\mathrm{d}}{\mathrm{d}r}\\!\left(\frac{G_{\scriptscriptstyle+}(r)}{r}\right)\int_{0}^{\infty}\\!\mathrm{d}r^{\prime}\,F(r^{\prime})G_{\scriptscriptstyle-}(r^{\prime})\\\
&\qquad-r\hskip
0.5pt\frac{\mathrm{d}}{\mathrm{d}r}\\!\left(\frac{G_{\scriptscriptstyle+}(r)}{r}\right)\int_{0}^{r}\\!\mathrm{d}r^{\prime}\,F(r^{\prime})G_{\scriptscriptstyle-}(r^{\prime})+r\hskip
0.5pt\frac{\mathrm{d}}{\mathrm{d}r}\\!\left(\frac{G_{\scriptscriptstyle-}(r)}{r}\right)\\!\int_{0}^{r}\\!\mathrm{d}r^{\prime}\,F(r^{\prime})G_{\scriptscriptstyle+}(r^{\prime})\Bigg{]}\end{aligned}\,\,,$
(3.17)

the combinations

$r\frac{\mathrm{d}}{\mathrm{d}r}\\!\left(\frac{G_{\pm}(r)}{r}\right)=C_{{\scalebox{0.75}[0.75]{$\pm$}},1}\sum_{k=0}^{\infty}a^{\pm}_{k}kr^{k}+C_{{\scalebox{0.75}[0.75]{$\pm$}},2}\sum_{k=0}^{\infty}b^{\pm}_{k}(1+k\log
r)r^{k}\,,$ (3.18)

are crucially regular as $r\to 0$, approaching a constant $b^{\pm}_{0}$.
Furthermore, since each of the integrals in the second line of (3.17) decay as
$\propto r^{2}\log r$ as $r\to 0$, each of those terms vanish as $r\to 0$ and
thus do not generate a non-normalizable contribution to $H(r)$.

In contrast, the terms in the first line are dangerous: they can contain a
constant piece which, when divided by $f(r)$, would cause $H(r)\sim H_{0}/r$
as $r\to 0$ for some constant $H_{0}$. Note, however, that we can always111The
only way this strategy could fail is if both $G_{\pm}(r)$ also have definite
scaling behavior as $r\to 0$, such that either
$a_{k}^{\scriptscriptstyle+}=b_{k}^{\scriptscriptstyle-}=0$. In this case, we
cannot tune $C_{\scriptscriptstyle-}$ to cancel the constant and the solution
is non-normalizable. While this behavior is extremely non-generic, we would
need the connection formulae for (3.9) to prove it does not happen. We take
our numerical solutions to be proof that it does not. choose
$C_{\scriptscriptstyle-}$ so that the overall constant term in the square
braces of (3.17) cancels. In this case, $H(r)\sim H_{0}\log r$ as $r\to 0$,
for some constant $H_{0}$. But, this is still square-integrable, and so we
conclude that we can construct, at least perturbatively, a normalizable zero
mode solution in the presence of a small, non-zero Dirac mass $m$. As in the
unperturbed case reviewed in Section 2, the axion string modifies the wave
equation so that both linearly-independent solutions are regular as $r\to 0$,
allowing for a normalizable zero mode.

The numerical solutions presented in Figure 2 demonstrate that these solutions
exist non-perturbatively as well, as long as $m<\mu$. As $m\to\mu$, these
solutions become unbound from the string and occupy a two-dimensional “wedge”
in the plane orthogonal to the string. To better understand the critical case
$m=\mu$, and why these solutions exist at all, it will be helpful to first
study a simpler, albeit analogous, system in which normalizable zero modes
exist even in the presence of a nowhere vanishing gap.

###  3.2 Zero Modes on Domain Walls

That a massless field can emerge from one that is everywhere massive is
counterintuitive, so it will be helpful to consider a simpler system where it
also occurs: the Dirac fermion on the half-plane
$\mathbb{R}_{\scriptscriptstyle+}^{2,1}$ [36, 63, 64, 65]. Studying this
system222This system naturally appears, for instance, in the study of the
quantum anomalous hall effect [66] along interfaces between topological and
regular insulators [65], and in the study of charge fractionalization in
polyacetylene chains [67, 68]. will also help us understand why these zero
modes completely delocalize onto a two-dimensional wedge when $m=\mu$, which
we will discuss in detail in §3.3.

Let us first consider a single Dirac fermion in the full space
$\mathbb{R}^{2,1}$ with a spatially-varying mass, and Lagrangian

$\mathcal{L}=\bar{\psi}\\!\left(i\not{\partial}-m(x)\right)\\!\psi\,.$ (3.19)

Here,
$(\gamma^{0},\gamma^{1},\gamma^{2})=(\sigma_{1},i\sigma_{3},{\scalebox{0.75}[1.0]{$-$}}i\sigma_{2})$
are the three-dimensional $\gamma$-matrices and we use the Cartesian
coordinates $(t,x,z)$. We take the mass term $m(x)$ to only depend on the
coordinate $x$. We can search for zero modes that move at the speed of light
in the $(+z)$-direction by assuming that there exists a solution of the form

$\psi(x)=\sqrt{p}\,{{\mathrm{e}}}^{{\scalebox{0.5}[0.85]{$-$}}ip(t-z)}\\!\begin{pmatrix}\eta(x)\\\
\chi(x)\end{pmatrix}\,.$ (3.20)

with $\eta(x)$ and $\chi(x)$ functions of $x$. With this ansatz, (3.19)
reduces to

$\begin{aligned} 0&=\eta^{\prime}+m(x)\eta-2p\chi\\\
0&=\chi^{\prime}-m(x)\chi\end{aligned}\,,$ (3.21)

where the ′ denotes differentiation with respect to $x$. This solution must
have definite helicity, with $\chi(x)=0$, so that solutions take the form

$\eta(x)=\mathcal{A}\hskip
0.5pt\exp\\!\left(-\\!\int_{0}^{x}\\!\mathrm{d}x^{\prime}\,m(x^{\prime})\right)\,,$
(3.22)

where the coefficient $\mathcal{A}$ is determined by a normalization condition
analogous to (2.15).

For (3.22) to yield a normalizable zero mode propagating in the
$(+z)$-direction, the mass $m(x)$ must be positive as $x\to\infty$ and
negative as $x\to{\scalebox{0.75}[1.0]{$-$}}\infty$, and thus _vanish_ for
some value of $x$. Likewise, if the mass is instead negative as $x\to\infty$
and positive as $x\to{\scalebox{0.75}[1.0]{$-$}}\infty$, there is a
normalizable zero mode propagating in the
$({\scalebox{0.75}[1.0]{$-$}}z)$-direction. Thus, if there exists a _domain
wall_ in which the fermion mass crosses through zero, there will be chiral
fermionic zero modes—often called chiral domain wall fermions—localized and
propagating along it.333This is a well-known fact which has been exploited
[63] to simulate chiral fermions on the lattice. Many of the phenomena we find
in this paper have analogs there. For instance, these chiral modes only exist
for momenta ${\mathbf{p}}$ in a particular region of the Brillouin zone [69],
which we can denote as ${\mathbf{p}}\in\mathcal{C}$. For
${\mathbf{p}}\notin\mathcal{C}$, the chiral modes cease to exist. In analogy
with Figure 2, the chiral modes have wavefunctions that are well-localized to
the domain wall for momenta comfortably inside $\mathcal{C}$, while they
completely delocalize from the domain wall as ${\mathbf{p}}$ approaches the
boundary $\partial\mathcal{C}$, and are no longer normalizable for
${\mathbf{p}}\notin\mathcal{C}$. We see that in $\mathbb{R}^{2,1}$, the
fermion cannot be everywhere gapped and still yield a normalizable zero mode.
This is one case in which vanishing mass and the existence of fermionic zero
modes are inextricably linked.

However, the situation changes if we restrict to the half-space
$\mathbb{R}^{2,1}_{\scriptscriptstyle+}$, with $x\geq 0$. In this case,
constant mass $m(x)=\mu$ _does_ yield a normalizable zero mode which is
exponentially localized to the boundary at $x=0$. For example, by decoupling
(3.21) we find that $\eta$ obeys

$\eta^{\prime\prime}+\big{(}m^{\prime}(x)-m^{2}(x)\big{)}\eta=0\,.$ (3.23)

To match onto the analysis of §2.2 and §3.1, we note that, since $m(x)$
approaches a constant $\mu$ at infinity, then $\eta(x)\sim
C_{\scriptscriptstyle-}{{\mathrm{e}}}^{{\scalebox{0.5}[0.85]{$-$}}\mu
x}+C_{\scriptscriptstyle+}{{\mathrm{e}}}^{\hskip 0.5pt\mu x}$ as $x\to\infty$.
Likewise, as long as $m(x)$ is regular as $x\to 0$, both linearly-independent
solutions will be regular, and so we will always be able to construct a
normalizable solution to (3.23). When $m(x)=\mu$ is constant everywhere on
$\mathbb{R}^{2,1}_{{\scriptscriptstyle+}}$, this solution is just
$\eta(x)=\mathcal{C}_{\scriptscriptstyle-}{{\mathrm{e}}}^{{\scalebox{0.5}[0.85]{$-$}}\mu
x}$. In contrast, this fails on the full-space $\mathbb{R}^{2,1}$ because the
solution that exponentially decays as $x\to\infty$ never matches onto the one
which decays as $x\to{\scalebox{0.75}[1.0]{$-$}}\infty$, unless the mass
$m(x)$ switches sign for some $x$.

To connect this simple system to the one we are interested in (3.2), we can
picture the polar radial coordinate in $\mathbb{R}^{3,1}$ as the analog of the
$x$ coordinate on the half-space $\mathbb{R}^{2,1}_{\scriptscriptstyle+}$,
while the axion string itself provides an effective “mass” $m(r)$ which is
regular as $r\to 0$, as are the solutions to (3.2). The same mechanism is at
work for both cases: both linearly-independent solutions are regular at the
core of the string or wall, and thus the solution that decays at spatial
infinity is necessarily normalizable. As we discuss in the next section, this
simple system is also useful for understanding the critical case in which core
and Yukawa masses are equal, $m=\mu$, and the zero modes become completely
delocalized from the string.

###  3.3 The Critical Mass Case

From the numerical results presented in Figure 2, we found that there is a
phase transition that occurs as we tune the core mass to the Yukawa mass,
$m\to\mu$, wherein the fermionic zero mode seemingly delocalizes from the
axion string. Beyond this critical mass, $m>\mu$, the zero mode ceases to
exist. Since this critical case turns out to be very difficult to study
numerically, it will instead be helpful to attack this case analytically to
understand exactly how these modes behave when $m=\mu$.

First, however, it will be useful to qualitatively understand how solutions to
(3.2) should behave in the limit $m\to\mu$. In the previous section, we
described how a fermionic zero mode can arise whenever a fermion’s spatially-
dependent mass crosses through zero. At criticality and far from the string,
the mass (3.3) approaches

$\lim_{r\to\infty}M(r,\varphi)=\sqrt{(m-\mu\cos\varphi)^{2}+\mu^{2}\sin^{2}\varphi}\,,$
(3.24)

which crosses through zero at $\varphi=0$ when (and only when) $m=\mu$. From
the logic of the previous section, we should then expect that the zero modes
are no longer localized along the axion string, but are instead allowed to
propagate freely along—and are localized to—the plane defined by $\varphi=0$.
Thus, we expect that these zero modes may have momentum along both the $x$\-
and $z$-directions, which we denote $p_{x}$ and $p_{z}$, respectively, with
frequency $\omega=\sqrt{p_{x}^{2}+p_{z}^{2}}$.

We can exhibit these solutions in the critical case by searching for solutions
to the equation of motion (3.2) far from the string along the positive
$x$-axis. In this limit, the equations of motion become approximately
translationally invariant in both the $x$\- and $z$-directions, and so we may
search for solutions of the form

$\psi(x)=\psi(x,y)\hskip 0.5pt\hskip
0.5pt{{\mathrm{e}}}^{{\scalebox{0.5}[0.85]{$-$}}i\omega t+ip_{x}x+ip_{z}z}\,,$
(3.25)

where $\psi(x,y)$ varies slowly along the $x$-direction. Since these modes
must have definite helicity, we will assume an ansatz of the form

$\psi(x,y)\approx\sqrt{\frac{\omega}{2}}\begin{pmatrix}\cos\alpha\\\
\sin\alpha\\\ i\sin\alpha\\\
{\scalebox{0.75}[1.0]{$-$}}i\cos\alpha\end{pmatrix}F(x,y)$ (3.26)

with $\cos 2\alpha={\scalebox{0.75}[1.0]{$-$}}p_{z}/\omega$, and $\sin
2\alpha={\scalebox{0.75}[1.0]{$-$}}p_{x}/\omega$. With this ansatz, when far
from the string (3.2) reduces to

$\partial_{y}F(x,y)\sim-\mu(y/x)F(x,y)\,,\mathrlap{\qquad x\gg|y|\,,}$ (3.27)

where we have dropped derivatives with respect to $x$ and approximated
${{\mathrm{e}}}^{i\varphi}\approx 1+i(y/x)$. This can be solved to find

$F(x,y)\approx\mathcal{A}(x)\,\exp\\!\hskip 0.5pt\left(-\\!\hskip
0.5pt\int_{0}^{|y|}\\!\mathrm{d}y^{\prime}\,M(x,y)\right)\sim\left(\frac{\mu}{\pi
x}\right)^{\frac{1}{4}}\exp\\!\hskip 0.5pt\left(-\frac{\mu y^{2}}{2x}\right),$
(3.28)

for $x\gg|y|$, where
$M(x,y)=\mu\sqrt{(1-\cos\varphi)^{2}+\sin^{2}\varphi}\approx\mu(y/x)$ is the
mass (3.24) far from the string, where $\mu x\gg 1$. The overall amplitude
$\mathcal{A}(x)$ is determined by imposing the normalization condition,

$\int\\!\mathrm{d}^{3}x\,\psi_{{\mathbf{p}}}^{\dagger}(x)\psi_{{\mathbf{p}}^{\prime}}(x)=(2\pi)^{2}\omega\,\delta^{(2)}({\mathbf{p}}-{\mathbf{p}}^{\prime})$
(3.29)

at any constant $t$, the direct analog of (2.15). For the ansatz (3.26), this
translates into the requirement that $\int\\!\mathrm{d}y\,|F(x,y)|^{2}=1$.

Figure 3: Plot of the approximate critical wavefunction $|F(x,y)|$, trustable
for large positive $x$, with contours displayed in units of $\sqrt{\mu}$. When
$p_{x}=0$ and $p_{z}={\scalebox{0.75}[1.0]{$-$}}\omega$,
$|F(x,y)|=|\psi_{0}(x,y)|$, and so this figure can be directly compared to the
bottom left panel of Figure 2.

We plot this approximate critical wavefunction (3.28) in Figure 3 and find
that it matches with the $m\to\mu$ behavior we observe in the numerical
results shown in Figure 2. As our qualitative arguments suggested, we find
that the zero modes delocalize from the string once $m=\mu$ and instead are
free to move along the half-plane defined by $\varphi=0$. They are relatively
well-localized along $y=\varphi=0$, but become more and more spread out the
further we get from the string, with an approximate width of $\Delta
y\sim\sqrt{x}/\mu$.

To summarize, in this section we studied how fermionic zero modes along an
axion string respond to the addition of a non-zero Dirac mass $m$. Absent this
core mass, the four-dimensional fermion sees effectively zero mass at the core
of the string, and so one might expect that this vanishing mass explains why
the axion string supports fermionic zero modes. However, we showed that these
zero modes exist even when the four-dimensional fermion has a core mass $m$
and is thus everywhere massive, and we explained what it is about the string
that allows these modes to exist. Furthermore, we found that these zero modes
completely delocalize from the string as $m\to\mu$, and cease to exist for
$m>\mu$. In the next section, we explain this phase structure from the
perspective of anomaly inflow, in which these zero modes are necessary to
render the low-energy effective theory consistent.

##  4 Zero Modes from Anomaly Inflow

The existence of these zero modes can also be inferred on topological grounds,
based on the logic of anomaly inflow [36]. Under an infinitesimal
${\mathrm{U}}(1)$ gauge transformation, the fields in (3.1) transform as

$\delta_{\Lambda}\psi=ie\Lambda(x)\psi\,,\qquad\delta_{\Lambda}A_{\mu}=\partial_{\mu}\Lambda(x)\,,$
(4.1)

where $\Lambda(x)$ is the ${\mathrm{U}}(1)$ gauge parameter. In a
topologically trivial background, the action (3.1) is manifestly invariant
under such a transformation. In the presence of the axion string, however,
gauge invariance is more subtle. This symmetry is _anomalous_ in the presence
of the string, such that electric charge is not conserved in a region
localized to the string. The theory is then inconsistent unless there are some
anomalous degrees of freedom, i.e. charged chiral excitations, that are
localized to the string and can cancel this anomaly and carry away electric
charge. These zero modes then allow electric charge and the anomaly to “flow”
out of the bulk spacetime and onto the string, ensuring that the full theory
remains consistent.

The general strategy is as follows. We will attempt to construct a low-energy
effective field theory of the axion and abelian gauge field by integrating out
both the fermion $\psi(x)$ and radial mode $f(x)$ of the complex scalar. We do
this by first _ignoring_ the contribution from possible fermionic zero modes
and by carefully considering the effect of a gauge transformation on the
effective theory. Our arguments are similar to the original story, which can
be found in [36, 38, 43], except that the ${\mathrm{U}}(1)_{{\textsc{pq}}}$
breaking mass leads to a subtlety in the identification of the phase which
should be rotated away. Depending on the relative size of $m$ and $\mu$, we
find a gauge anomaly localized to the string. This inconsistency then forces
us to include an additional contribution from fermionic zero modes localized
to the string which cancels this anomaly.

We are interested in the gauge invariance of the low-energy effective theory
for the axion field after integrating out the fermions. In the string
background (2.3), the fermion sector of the theory can be written

$\mathcal{L}\supset\bar{\psi}\big{(}i\not{D}-m+yf(r){{\mathrm{e}}}^{i\gamma^{5}\theta(x)}\big{)}\psi,$
(4.2)

with $\theta(x)$ the axion field, and so we wish to compute

$\mathcal{Z}_{\psi}(\theta)=\int\\!\mathcal{D}\bar{\psi}\mathcal{D}\psi\,\exp\\!\bigg{[}\hskip
0.5pt\hskip 0.5pti\\!\int\\!{\rm
d}^{4}x\,\bar{\psi}\big{(}i\not{D}-M(r,\theta)\hskip 0.5pt\hskip
0.5pt{{\mathrm{e}}}^{i\gamma^{5}\alpha(r,\theta)}\big{)}\psi\bigg{]}\,,$ (4.3)

with $M(r,\theta)$ and $\alpha(r,\theta)$ defined in (3.3). Performing the
path integral by introducing a Pauli-Villars regulator with mass $\tilde{M}$,
we find

$\mathcal{Z}_{\psi}(\theta)=\frac{\det(i\not{D}-M(r,\theta){{\mathrm{e}}}^{i\gamma^{5}\alpha(r,\theta)})}{\det(i\not{D}-\tilde{M})}\,.$
(4.4)

Now we perform a spatially-dependent chiral field redefinition,
$\psi\to{{\mathrm{e}}}^{{\scalebox{0.5}[0.85]{$-$}}i\gamma^{5}\alpha(r,\theta)/2}\psi$
to try and remove $\mathcal{Z}_{\psi}(\theta)$’s dependence on
$\alpha(r,\theta)$. However, due to the chiral anomaly, this transformation
introduces a Jacobian factor and we thus have

$\mathcal{Z}_{\psi}(\theta)=\frac{\det(i\not{D}-M(r,\theta))}{\det(i\not{D}-\tilde{M})}\hskip
0.5pt\exp\bigg{[}\frac{i}{8\pi^{2}}\\!\int\\!\alpha(r,\theta)F\wedge
F\bigg{]},$ (4.5)

which can be written in the form of an effective action for the axion $\theta$
interacting with the gauge field $A_{\mu}$. Note that in the $m\to 0$ limit,
in which there is a classical $U(1)_{\textsc{pq}}$ symmetry,
$\alpha(r,\theta)\to\theta$ and this Jacobian factor reduces to the usual
quantized coupling of the axion to the gauge field, $\propto\theta F\wedge F$.
Furthermore, integrating out the fermions leads to a nontrivial effective
potential for the axion via the $\theta$-dependence in $M(r,\theta)$.

Figure 4: The phase of fermion “mass” (4.6) far from the string. The phase is
double-valued at $\theta=2\pi\sim 0$ for $m\leq\mu$ [ blue,  orange,  green].
The zero mode solutions vanish as $m>\mu$ [ red] once the phase becomes
single-valued.

It is this additional $\alpha(r,\theta)F\wedge F$ term in the effective action
that leads to a potential issue with gauge invariance. Whether or not this
term is gauge invariant is determined by whether the phase far from the
string,

$\alpha(\theta)=\lim_{r\to\infty}\alpha(r,\theta)=\arg(m-\mu\hskip
0.5pt{{\mathrm{e}}}^{{\scalebox{0.5}[0.85]{$-$}}i\theta})\,,$ (4.6)

is single- or multi-valued as a function of $\theta$, and this depends on the
relative size of $m$ and $\mu$ as illustrated in Figure 4. For instance, in
the limit that $m\ll\mu$, we find

$\alpha(\theta)\approx\pi-\theta,$ (4.7)

which is clearly not single-valued at $\theta=2\pi\sim 0$, and this applies
more generally for all $m<\mu$. In this case, we must integrate the effective
action by parts to make it well-defined,

$\frac{1}{8\pi^{2}}\\!\int\\!\alpha(\theta)F\wedge
F\to-\frac{1}{8\pi^{2}}\\!\int\\!{\rm d}\alpha(\theta)\wedge A\wedge F\,.$
(4.8)

However, under the gauge transformation (4.1), this transforms as

$\delta_{\Lambda}\bigg{[}-\frac{1}{8\pi^{2}}\\!\int\\!{\rm
d}\alpha(\theta)\wedge A\wedge F\bigg{]}=-\frac{1}{8\pi^{2}}\\!\int\\!{\rm
d}\alpha(\theta)\wedge{\rm d}(\hskip 0.5pt\Lambda F\hskip
0.5pt)=+\frac{1}{8\pi^{2}}\\!\int\\!{\rm d}^{2}\alpha(\theta)\wedge(\hskip
0.5pt\Lambda F\hskip 0.5pt)\,,$ (4.9)

where we have integrated by parts again in the last equality, neglecting any
potential surface terms at the string core. Since $\alpha(\theta)$ is not
single-valued, ${\rm d}^{2}\alpha(\theta)\neq 0$ but is instead localized
along the string, and we find an apparent gauge anomaly.

The reason for this apparent anomaly is that we have heretofore incorrectly
assumed that there are no zero modes at the core of the string. Their
existence modifies the chiral transformation in (4.5) and, with a careful
treatment of the variation at the core,444 Somewhat famously, arriving at the
exact value of the gauge anomaly, which matches the contribution from two-
dimensional chiral fermions localized on the string worldsheet, requires some
care. A proper treatment involves introducing a “bump-form” which extends the
validity of the effective theory in (4.5) to all of spacetime and properly
includes the zero modes on the string core so that all surface contributions
vanish [43]. This bump form can be computed explicitly in terms of the zero
modes profiles, and the calculation of [43] is unchanged in the case $m\neq
0$. we find a violation of gauge invariance precisely equal and opposite to
the gauge anomaly due to a massless chiral fermion localized to the string
worldsheet,

$\delta_{\Lambda}\bigg{[}-\frac{1}{8\pi^{2}}\\!\int\\!{\rm
d}\alpha(\theta)\wedge A\wedge F\bigg{]}=\frac{1}{4\pi}\Lambda\,F\,.$ (4.10)

On the other hand, it is precisely when $m>\mu$ that the phase
$\alpha(\theta)$ becomes single-valued since, near the end points at
$\theta=0$ and $2\pi$, we have $\alpha(0)=\alpha(2\pi)=\arg(m-\mu)=0$. In this
case the interaction $\propto\alpha(\theta)\,F\wedge F$ is both well-defined
and single-valued, and gauge invariance is maintained without the need for
additional degrees of freedom along the string. This analysis breaks down as
$m\to\mu$ because it no longer makes sense to define an effective field theory
for the axion and gauge field “far” from the string since, as evident from our
numerical results in Figure 2 and discussed in §3.3, the fermionic zero modes
are no longer localized at small $r$.

Finally, it is illuminating to compare our effective action to the effective
action computed via the method of Goldstone and Wilczek [70], as in [36, 38].
In that case, one computes the current at a point $x_{0}$ far from the string
core, in response to a background electromagnetic field, in the long
wavelength limit. We choose $\theta(x_{0})=0$, so that the fermion-scalar
interaction can be expanded as

$y\bar{\psi}f(r)e^{i\gamma_{5}\theta(x)}\psi\simeq
i\mu\,\partial_{\lambda}\theta(x_{0})(w-x_{0})^{\lambda}\bar{\psi}\gamma^{5}\psi.$
(4.11)

The leading contribution comes from the diagram in Figure 5. For the theory
(3.1), this yields

$\displaystyle\langle
J^{\mu}(z)\rangle=e^{2}\mu\,\partial_{\lambda}\theta(x_{0})\\!\int\frac{{\rm
d}^{4}p}{(2\pi)^{4}}$ $\displaystyle\frac{{\rm d}^{4}q}{(2\pi)^{4}}\frac{{\rm
d}^{4}k}{(2\pi)^{4}}\,{\rm d}^{4}y\,{\rm
d}^{4}w\,(w-x_{0})^{\lambda}A_{\nu}(y)$ (4.12)
$\displaystyle\times{{\mathrm{e}}}^{{\scalebox{0.5}[0.85]{$-$}}i[p(z-y)+(p+q)(y-w)+(p+q+k)(w-z)]}$
$\displaystyle\times{\rm
Tr}\bigg{[}\frac{1}{\not{p}-M}\gamma^{\nu}\frac{1}{\not{p}+\not{q}-M}\gamma^{5}\frac{1}{\not{p}+\not{q}+\not{k}-M}\gamma^{\mu}$
$\displaystyle\qquad\quad+\frac{1}{\not{p}-M}\gamma^{5}\frac{1}{\not{p}+\not{k}-M}\gamma^{\nu}\frac{1}{\not{p}+\not{k}+\not{q}-M}\gamma^{\mu}\bigg{]}$

and evaluates to

$\langle
J^{\mu}(z)\rangle=\frac{e^{2}}{8\pi}\frac{\mu}{\mu-m}\epsilon^{\mu\nu\kappa\lambda}\partial_{\nu}\theta(x_{0})F_{\kappa\lambda}\,.$
(4.13)

Following [38], we can posit an effective action (valid far from the string),

$S_{\mathrm{eff}}=\frac{1}{8\pi^{2}}\frac{\mu}{\mu-m}\int\\!\hskip
0.5pt\mathrm{d}\theta\wedge A\wedge F\,,$ (4.14)

from which the current (4.13) can be computed. When $m=0$, integrating (4.14)
by parts yields the usual, properly quantized, coupling of the axion to the
gauge field. When $m\neq 0$, the $\mu/(\mu-m)$ prefactor would appear to
violate the quantization condition for the axion–gauge field coupling.
However, expanding for small values of the axion field,

${\rm
d}\alpha(\theta)=\left[\frac{\mu}{\mu-m}+\mathcal{O}(\theta^{2})\right]{\rm
d}\theta\,,$ (4.15)

and we see that (4.14) is precisely the leading term of (4.8), which is
properly quantized.

$\partial_{\nu}\theta(x_{0})$$A_{\lambda}(y)$$J_{\mu}(z)$ Figure 5: Leading
order diagram for the current $\langle J^{\mu}(z)\rangle$ produced by a
spatial variation of the axion $\partial_{\nu}\theta(x_{0})$ in the presence
of a background electromagnetic field $A_{\lambda}$, used by Goldstone and
Wilczek to compute the effective action.

Let us note that there is an ambiguity in which field one calls the “axion,”
$\theta$ or $\alpha(\theta)$. Both are compact fields, $\theta\sim\theta+2\pi$
and $\alpha\sim\alpha+2\pi$. If one works in the “charge-quantized” basis in
which the coefficient of the topological coupling (4.8) is an integer, the
axion $\alpha$ has a (possibly highly) non-trivial metric on its field space
given by

$g_{\alpha\alpha}=v^{2}(\theta^{\prime}(\alpha))^{2}=v^{2}\left[1+\frac{m}{\mu}\frac{\cos\alpha}{\sqrt{1-(m/\mu)^{2}\sin^{2}\alpha}}\right]^{2}\,,$
(4.16)

and so around $\theta\approx 0$ or $\alpha\approx\pi$ the decay constant is
effectively reduced by a factor of $(\mu-m)/\mu$, which is identical to the
factor that appears in (4.14). If one instead works with $\theta$, the
(classical) metric on field space is flat but the coefficient of the
topological coupling (4.14) is not an integer. This story is similar in spirit
to the reason for the non-quantized coupling of the axion to photons, due to
mixing with the neutral pion.

To summarize, we have shown that—with a careful treatment of the axion’s
periodicity—the arguments for anomaly inflow, and the resulting necessity of
chiral zero modes on the string, persist in the presence of a field-
independent mass, as long as this “bulk” mass is smaller than the mass
generated by the chirally coupled scalar field. This makes it clear that it is
the topology of the axion string, and not the dynamics of the scalar field
which spontaneously breaks the classical ${\mathrm{U}}(1)_{{\textsc{pq}}}$
symmetry, which is responsible for the existence of the fermion zero modes.

##  5 Low-Energy Effective Theory

In the previous section, we found that the full four-dimensional effective
theory was rendered consistent when $m<\mu$ by including chiral fermionic zero
mode localized to the axion string, while they disappeared for $m>\mu$. Our
numerical results in Figure 2 show that they disappear because the zero mode
completely delocalizes from the string as $m\to\mu$. It will be useful to
understand how this process appears from the zero mode’s two-dimensional
effective theory on the string worldsheet. At leading order in
$\mu^{{\scalebox{0.5}[0.85]{$-$}}1}$, this action is completely constrained by
symmetry and does not depend on the ratio $m/\mu$. However, we will show that
the Wilson coefficients of the higher derivative interactions between the zero
mode and the gauge field depend sensitively on $m/\mu$ and diverge as
$m\to\mu$, causing the two-dimensional effective theory to break down.

We again will specialize to an axion string with charge $+1$, and start by
introducing the (quantized) mode expansion for the fermion

$\psi(x)=\int_{0}^{\infty}\\!\frac{{\rm
d}p}{2\pi}\frac{1}{\sqrt{2|p|}}\,\hat{a}_{p}\hskip
0.5pt\psi_{p}(r,\varphi)\hskip
0.5pt{{\mathrm{e}}}^{{\scalebox{0.5}[0.85]{$-$}}ip(t+z)}+\int_{-\infty}^{0}\frac{{\rm
d}p}{2\pi}\frac{1}{\sqrt{2|p|}}\hat{b}_{{\scalebox{0.5}[0.85]{$-$}}p}^{\dagger}\hskip
0.5pt\psi_{p}(r,\varphi){{\mathrm{e}}}^{{\scalebox{0.5}[0.85]{$-$}}ip(t+z)}+\cdots$
(5.1)

where the $\cdots$ denote the non-zero modes of the massive four-dimensional
fermion about the axion string. Here, $\psi_{p}(r,\varphi)$ is the zero mode
wavefunction (2.10) found in Section 3, subject to the normalization condition
(2.15). The creation and annihilation operators, $\hat{b}^{\dagger}_{p}$,
$\hat{a}_{p}$ satisfy the canonical commutation relations appropriate for a
two-dimensional fermion,

$\\{\hat{a}_{p}^{\vphantom{\dagger}},\hat{a}_{p\smash{{}^{\prime}}}^{\dagger}\\}=2\pi\,\delta(p-p^{\prime}),\qquad\\{\hat{b}^{\vphantom{\dagger}}_{p},\hat{b}_{p\smash{{}^{\prime}}}^{\dagger}\\}=2\pi\,\delta(p-p^{\prime}).$
(5.2)

Since the only place the momentum $p$ enters into $\psi_{p}(r,\varphi)$ is in
its overall normalization, we can define the spinor
$\mathcal{F}(r,\varphi)=\psi_{p}(r,\varphi)/\sqrt{p}$ and write (5.1) as
$\psi(x)=\chi_{\scalebox{0.5}[0.85]{$-$}}(t,z)\mathcal{F}(r,\varphi)+\cdots$,
where

$\chi_{\scalebox{0.5}[0.85]{$-$}}(t,z)=\int_{0}^{\infty}\\!\frac{\mathrm{d}p}{2\pi}\,\Big{[}\hat{a}_{p}\hskip
0.5pt\hskip
0.5pt{{\mathrm{e}}}^{{\scalebox{0.5}[0.85]{$-$}}ip(t+z)}+\hat{b}^{\dagger}_{p}\hskip
0.5pt\hskip 0.5pt{{\mathrm{e}}}^{ip(t+z)}\Big{]}$ (5.3)

is a canonically normalized two-dimensional fermion field operator with
negative chirality.

It will be convenient to repackage $\chi_{\scalebox{0.5}[0.85]{$-$}}(t,z)$
into a two-dimensional Dirac fermion $\chi(t,z)$ living on the axion string
worldsheet

$\chi(t,z)=\begin{pmatrix}\chi_{\scriptscriptstyle-}(t,z)\\\
\chi_{\scriptscriptstyle+}(t,z)\end{pmatrix}\,,$ (5.4)

by grouping it with a positive chirality fermion
$\chi_{\scriptscriptstyle+}(t,z)$, which we later set to zero. We will use
$a,b,\ldots=0,3$ to denote worldsheet indices, and define our worldsheet
$\gamma$-matrices as
$\gamma^{a}=(\tilde{\gamma}^{0},\tilde{\gamma}^{3})=(\sigma_{1},i\sigma_{2})$555We
use $\tilde{\gamma}^{0}$ and $\tilde{\gamma}^{3}$ to distinguish these two-
dimensional $\gamma$-matrices from the four-dimensional $\gamma^{0}$ and
$\gamma^{3}$, used for the four-dimensional Dirac fermion $\psi(x)$.  with
chiral projectors
$\gamma_{\scalebox{0.75}[0.75]{$\pm$}}=\frac{1}{2}(\mathbbm{1}\pm\tilde{\gamma}^{0}\tilde{\gamma}^{3})$
such that
$\gamma_{\scalebox{0.75}[0.75]{$\pm$}}\chi=\chi_{{\scalebox{0.75}[0.75]{$\pm$}}}$.
We will thus impose the constraint $\gamma_{\scriptscriptstyle+}\chi=0$ to
remove $\chi_{\scriptscriptstyle+}$ from the theory. The kinetic term for the
two-dimensional fermion then takes the standard Dirac form,

$S_{\chi}=\int\\!\mathrm{d}^{2}\sigma\,i\bar{\chi}\gamma^{a}\partial_{a}\chi+\cdots=\int\\!\mathrm{d}^{2}z\,i\chi^{\dagger}_{\scriptscriptstyle-}\partial_{\scriptscriptstyle+}\chi_{\scriptscriptstyle-}+\cdots$
(5.5)

where we denote $\partial_{\pm}=\partial_{t}\mp\partial_{z}$ and use
$\mathrm{d}^{2}\sigma=\mathrm{d}t\,\mathrm{d}z$ to denote the volume element
on the string worldsheet, while we will use
$\mathrm{d}^{2}r=r\,\mathrm{d}r\,\mathrm{d}\varphi=\mathrm{d}x\,\mathrm{d}y$
to denote the volume element orthogonal to the string. The $\cdots$ denote
terms in the effective action that encode how the chiral zero mode
$\chi_{\scriptscriptstyle-}$ interacts with the gauge field $A_{\mu}$, which
we now derive.

The zero mode generates a four-dimensional current that points along the axion
string

$j^{\mu}=e\bar{\psi}\gamma^{\mu}\psi=\begin{dcases}e|\mathcal{F}(r,\varphi)|^{2}(\bar{\chi}\gamma^{a}\chi)(t,z)&\mu=a=0,3\\\
0&\mu=1,2\end{dcases}\,,$ (5.6)

where
$|\mathcal{F}(r,\varphi)|^{2}=(\mathcal{F}^{\dagger}\mathcal{F})(r,\varphi)$,
and so its interaction with the gauge field $A_{\mu}$ is determined by

$S_{\textsc{em}}=-\\!\int\\!\mathrm{d}^{4}x\,j^{\mu}A_{\mu}=-\\!\int\\!\mathrm{d}^{4}x\,e|\mathcal{F}(r,\varphi)|^{2}(\bar{\chi}\gamma^{a}\chi)(t,z)A_{\mu}(x)\,.$
(5.7)

Our goal is to encode the information contained in (5.7) into a series of
effective interactions in the two-dimensional effective theory between $\chi$
and (derivatives of) the gauge field pulled back onto the string worldsheet.
This will encode the zero mode wavefunction’s multipolar structure into a set
of effective two-dimensional interactions and will define a set of multipole
moments for the current (5.6).

To derive these effective interactions, we perform a long-wavelength expansion
of the gauge field about the axion string $x=(t,0,0,z)$ following a procedure
similar to the one outlined in [71] and explained in more detail in Appendix
B. After using current conservation and integration by parts, (5.7) can be
written as

$\displaystyle S_{\textsc{em}}$
$\displaystyle=-\int\\!\mathrm{d}^{2}\sigma\big{[}e(\bar{\chi}\gamma^{a}\chi)A_{a}\big{]}\\!\hskip
0.5pt(t,z)$ (5.8)
$\displaystyle-\sum_{n=1}^{\infty}\frac{1}{n!}\\!\left[\int\\!{\rm
d}^{2}r\,x^{k_{1}}\cdots
x^{k_{n}}\,|\mathcal{F}(r,\varphi)|^{2}\right]\\!\int\\!\mathrm{d}^{2}\sigma\big{[}e(\bar{\chi}\gamma^{a}\chi)\partial_{k_{2}}\cdots\partial_{k_{n}}F_{k_{1}a}\big{]}\\!\hskip
0.5pt(t,z),$

where we use $i,j,k_{1},\cdots=1,2$ indices to denote directions orthgonal to
the string, i.e. ${x^{i}=(x,y)}$. It is more illuminating to rewrite the
higher order terms in this expansion in terms of irreducible
${\mathrm{SO}}(2)$ tensors. This introduces terms $\propto\partial^{k}F_{ka}$,
which can be rewritten using the equations of motion that follow from (4.5),

$\partial^{k}F_{ka}={\scalebox{0.75}[1.0]{$-$}}\partial^{b}F_{ba}+\frac{e^{2}}{16\pi^{2}}\partial^{\mu}(\alpha(\theta)\tilde{F}_{a\mu})+j_{a}\,,$
(5.9)

where
$\tilde{F}_{\mu\nu}=\tfrac{1}{2}\epsilon_{\mu\nu\rho\sigma}F^{\rho\sigma}$ is
the dual field strength. To second order in the derivative expansion, we find
that the two-dimensional current $\tilde{\jmath}^{\hskip 0.5pt\hskip
0.5pta}\equiv\bar{\chi}\gamma^{a}\chi$ couples to

$\displaystyle
S_{\textsc{em}}\supset\int\\!\mathrm{d}^{2}\sigma\\!\left[-\tilde{\jmath}^{\hskip
0.5pt\hskip 0.5pta}A_{a}+I_{1}^{i}\,\tilde{\jmath}^{\hskip 0.5pt\hskip
0.5pta}F_{ia}+\frac{1}{4}I_{2}^{ij}\tilde{\jmath}^{\hskip 0.5pt\hskip
0.5pta}\partial_{(i}F_{j)a}-\frac{1}{4}I_{2}^{\mathrm{tr}}\tilde{\jmath}^{\hskip
0.5pt\hskip
0.5pta}\partial^{b}F_{ba}+\frac{e^{2}}{64\pi^{2}}I_{2}^{\mathrm{tr}}\tilde{\jmath}^{\hskip
0.5pt\hskip 0.5pta}\partial^{\rho}(\theta\tilde{F}_{a\rho})\right]$ (5.10)

where

$I_{1}^{i}=\int\\!{\rm d}^{2}r\,x^{i}|\mathcal{F}(r,\varphi)|^{2}$ (5.11)

is the zero mode current’s “dipole,” while

$\begin{gathered}\quad I_{2}^{ij}=\int\\!{\rm
d}^{2}r\,(x^{i}x^{j}-\tfrac{1}{2}\delta^{ij}x_{k}x^{k})|\mathcal{F}(r,\varphi)|^{2},\qquad
I_{2}^{\mathrm{tr}}=\int\\!\mathrm{d}^{2}r\,r^{2}|\mathcal{F}(r,\varphi)|^{2}\end{gathered}\,,$
(5.12)

are the current’s “quadrupole” and variance. Each of the operators in (5.10)
are operators on the string worldsheet, while the dimensionless coefficients
$c_{1}^{i}\equiv\mu I^{i}_{1}$, $c_{2}^{ij}\equiv\mu^{2}I^{ij}_{2}$,
$c_{2}^{\mathrm{tr}}\equiv\mu^{2}I_{2}^{\mathrm{tr}}$ can be understood as
Wilson coefficients in the effective theory on the string worldsheet.666There
are only three independent Wilson coefficients that appear at each order in
this multipole expansion. This can be most easily seen by remembering that the
irreducible tensors of ${\mathrm{SO}}(2)$ with weight $\ell$ are in one-to-one
correspondence with the two Fourier modes ${{\mathrm{e}}}^{\pm i\ell\varphi}$
and so, in general, correspond to two complex (or four real) degrees of
freedom, which we may denote as $f_{\ell}$ and
$f_{{\scalebox{0.5}[0.85]{$-$}}\ell}$. Since the distribution is purely real,
these coefficients obey $f_{\ell}=f_{{\scalebox{0.5}[0.85]{$-$}}\ell}^{*}$,
and so the distribution is described by only two real degrees of freedom at
each weight $\ell$. To this, we must add the pure trace term at each $\ell$,
e.g. $I_{2}^{\mathrm{tr}}$, which increases the number of moments at each
order to three.

Figure 6: Plots of the Wilson coefficients $c_{1}^{x}$, $c_{2}^{xx}$, and
$c_{2}^{\mathrm{tr}}$ vs. $m/\mu$. The coefficients $c_{1}^{y}$ and
$c_{2}^{xy}=c_{2}^{yx}$ vanish due to the zero mode’s reflection symmetry
about the $x$-axis, as evident in Figure 2, while the other coefficient
$c_{2}^{yy}$ is completely determined by
$c_{2}^{yy}={\scalebox{0.75}[1.0]{$-$}}c_{2}^{xx}$.

We plot these Wilson coefficients in Figure 6 as functions of $m/\mu$. Two
features are immediately obvious: at $m=0$, the dipole and quadrupole moments
vanish due to the cylindrical symmetry of the zero mode solution, as expected.
The trace term $c_{2}^{{\mathrm{tr}}}$ does not vanish but instead approaches
a value $c_{2}^{\mathrm{tr}}\to 2.19$ as $m\to 0$. Such trace terms also
appear in the generic multipole expansion of a point-like source [71] where,
as here, they multiply _time_ derivatives of the field strength and thus do
not appear in the static multipole expansion.

As $m$ approaches the critical value, $m=\mu$, all three Wilson coefficients
quickly blow up, reflecting the breakdown of the two-dimensional effective
theory. This is also expected since, as seen in Figure 2, the zero modes
delocalize from the string as $m\to\mu$ and must completely disappear for
$m>\mu$, and so the worldsheet effective theory must fail. As is usual, the
effective theory signals its own demise through Wilson coefficients that
become uncontrollably large.

So, far from a charge $+1$ axion string the theory (3.1) is well-described by
the effective action

$\begin{aligned}
S&=\int\\!\mathrm{d}^{4}x\left[-\frac{1}{4}F_{\mu\nu}F^{\mu\nu}+v^{2}(\partial\theta)^{2}+\frac{\alpha(\theta)}{16\pi^{2}}F_{\mu\nu}\tilde{F}^{\mu\nu}+\cdots\right]\\\
&+\int\\!\mathrm{d}^{2}\sigma\left[\bar{\chi}i\not{D}\chi+\frac{c_{1}^{i}}{\mu}\tilde{\jmath}^{\hskip
0.5pt\hskip 0.5pta}F_{ia}+\frac{c_{2}^{ij}}{4\mu^{2}}\tilde{\jmath}^{\hskip
0.5pt\hskip
0.5pta}\partial_{(i}F_{j)a}-\frac{c_{2}^{\mathrm{tr}}}{4\mu^{2}}\tilde{\jmath}^{\hskip
0.5pt\hskip
0.5pta}\partial^{b}F_{ba}+\frac{e^{2}c_{2}^{\mathrm{tr}}}{64\pi^{2}\mu^{2}}\tilde{\jmath}^{\hskip
0.5pt\hskip
0.5pta}\partial^{\rho}(\theta\tilde{F}_{a\rho})+\cdots\right]\end{aligned}\,.$
(5.13)

The first line describes the dynamics of the gauge field and axion, where the
$\cdots$ denote interactions that are induced by, for instance, integrating
out the scalar’s radial mode $f$ and the non-zero mode fluctuations of the
fermions and include a potential for the axion $\theta$. The second line
describes the chiral zero mode living on the string and its interactions with
the gauge field and axion. Here, the $\cdots$ denote terms that are either
subleading in the derivative expansion or are nonlinear in the current
$\tilde{\jmath}^{\hskip 0.5pt\hskip 0.5pta}$, the latter of which are also
generated by the long-wavelength expansion described above and in Appendix B
and upon integrating out the heavy scalar and fermionic modes.

##  6 Conclusions

In this work, we revisited the superconductivity of string solutions in a
simple model of axion electrodynamics in which the
${\mathrm{U}}(1)_{\textsc{pq}}$ symmetry is explicitly broken by a mass term
for the fermions. When the PQ-breaking mass $m$ is smaller than the asymptotic
mass $\mu$ acquired from the radial mode of the scalar whose phase is the
axion, the fermionic zero modes responsible for superconductivity persist. We
demonstrated the existence of these zero modes both by studying the asymptotic
behavior of the equations of motion and by solving for their profiles
numerically. We also demonstrated how their existence can be understood from
anomaly inflow, with some modifications to the original arguments by Callan
and Harvey.

For $m>\mu$, the zero modes cease to exist. As $m$ approaches $\mu$, however,
the zero modes exhibit an interesting critical behavior, in which they
delocalize from the core of the string and propagate along a two-dimensional
wedge.777Amusingly, this behavior is very similar to how light fermions behave
in the presence of a magnetic monopole [72], but in reverse. There, the light
fermionic modes become more and more delocalized from the monopole core as
their mass _decreases_ , ultimately explaining why the $\theta$-angle becomes
a redundant parameter for dyon physics once $m=0$. Here, counterintuitively,
the light fermionic modes delocalize as their mass _increases_. We further
studied the effective theory of the zero modes on the worldsheet and their
interactions with external gauge fields. As $m$ approaches its critical value,
we demonstrated explicitly how this effective theory breaks down as a result
of the zero modes delocalizing.

While we have focused on a simple “minimal” model of axion superconductivity
here, there are several aspects which may be phenomenologically relevant in
more realistic theories that warrant further exploration. First, as mentioned
in the introduction, DFSZ axion models have multiple distinct, topologically
stable string configurations and at least one of these (the “Type-C strings”
in [52]) has a scalar profile which does not restore electroweak symmetry at
the string core. It is possible that the nonzero vacuum expectation value at
the core in these configurations may play a similar role to the explicit PQ-
breaking mass term we have considered, which helps explain how these solutions
may still be superconducting (as anomaly inflow arguments suggest they must
be). This should be studied in more detail.

It is also worth recalling that in a high-temperature background, fermions
acquire a Debye mass $\sim gT$ which has the same effect as the explicit mass
$m$. In an expanding universe, if the initial temperature is large enough that
the Debye mass is larger than the Yukawa mass, a phase transition to the
superconducting case may actually occur. It may thus be of interest to study
if the delocalization of the zero modes found in the simple model studied here
can have consequences for the evolution of string networks in cosmology.

Finally, we should emphasize that our numerical results were derived under the
assumption that the fermion amplitude is small. It would be worthwhile to
understand how the profiles of the zero modes change beyond the limit of small
fluctuations. In the same vein, it would be enlightening to understand the
higher-order source terms in the worldsheet effective action that we have
dropped. We hope to return to some of these topics in the future.

## Acknowledgments

We are grateful to Clay Cordova, Nicolas Fernandez, Junwu Huang, Austin Joyce,
Seth Koren, Ho Tat Lam, Cody Long, Rashmish Mishra, Matthew Reece, Shu-Heng
Shao and Tanmay Vachaspati for useful conversations. We are also especially
grateful to Matthew Reece for comments on the manuscript. HB, KF, SH and JS
are supported in part by the DOE grant DE-SC0013607. KF and SH are also
supported in part by the Alfred P. Sloan Foundation Grant No. G-2019-12504.
HB, KF and JS are also supported in part by NASA Grant 80NSSC20K0506.

##  Appendix A Numerical Techniques

We rely on Chebyshev interpolation to numerically solve the zero mode
equations of motion (3.4). This method is one of several numerical approaches
used for collocation: the discretization of partial differential equations
(PDEs) to transform them into matrix equations. Here we give a brief
introduction; more details can be found in [73, 74].

Chebyshev interpolation provides a method of approximating a smooth function
$f(\zeta)$ on the interval $\zeta\in[{\scalebox{0.75}[1.0]{$-$}}1,1]$ using
its values at a finite and discrete set of points. That is, given the function
$f(\zeta)$ evaluated at the set of _Cheybshev points_ ,

$\zeta_{k}=\cos\left(\frac{\pi(2k+1)}{2(N+1)}\right)\,,\mathrlap{\qquad
k=0,\ldots,N}\,,$ (A.1)

we may approximate it at _any_ point in
$\zeta\in[{\scalebox{0.75}[1.0]{$-$}}1,1]$ via the interpolant

$f_{N}(\zeta)=\sum_{k=0}^{N}f(\zeta_{k})p_{k}(\zeta)$ (A.2)

where the functions $p_{k}(\zeta)$, defined such that
$p_{k}(\zeta_{n})=\delta_{kn}$, are the degree-$N$ Lagrange polynomials
associated to the points (A.1). Importantly, the Lagrange polynomials can be
reliably evaluated away from the interpolation points using the _second
barycentric form_ [75, 74],

$p_{k}(\zeta)=\frac{w_{k}}{\zeta-\zeta_{k}}\,\Bigg{/}\sum_{n=0}^{N}\frac{w_{n}}{\zeta-\zeta_{n}}\,,$
(A.3)

where we introduce the _Chebyshev weights_

$w_{k}=\sin\left(\frac{2\pi k(N+2)+\pi}{2(N+1)}\right)\,,\mathrlap{\qquad
k=0,\ldots,N\,.}$ (A.4)

The interpolant (A.2) then provides a degree-$N$ polynomial approximation that
exactly agrees with $f(\zeta)$ at the points $\zeta_{k}$ and has an error
which decays exponentially in $N$ for all other points on the interval.

The interpolant (A.2) allows us to represent the function $f(\zeta)$ as a
$(N+1)$-component vector $\vec{f}=\\{f(\zeta_{k})\\}$. We can similarly
represent derivatives of this function working with the derivative matrices
$p^{\prime}_{k}(\zeta_{n})$ and $p^{\prime\prime}_{k}(\zeta_{n})$, defined by
[74, 76]

$\displaystyle p^{\prime}_{k}(\zeta_{n})$
$\displaystyle=\begin{dcases}\frac{w_{k}/w_{n}}{\zeta_{n}-\zeta_{k}}&n\neq
k\\\ -\sum_{k\neq n}p^{\prime}_{k}(\zeta_{n})&n=k\end{dcases}\,\,\,,$ (A.5)
$\displaystyle p^{\prime\prime}_{k}(\zeta_{n})$
$\displaystyle=\begin{dcases}2p^{\prime}_{k}(\zeta_{n})p^{\prime}_{n}(\zeta_{n})-\frac{2p^{\prime}_{k}(\zeta_{n})}{\zeta_{n}-\zeta_{k}}&n\neq
k\\\ -\sum_{k\neq n}p^{\prime\prime}_{k}(\zeta_{n})&n=k\end{dcases}\,\,\,.$

Any differential equation for $f(\zeta)$ on the interval
$\zeta\in[{\scalebox{0.75}[1.0]{$-$}}1,1]$ can then be approximated by a
matrix equation for the function values $\\{f(\zeta_{k})\\}$ by substituting
(A.2) for $f(\zeta)$, evaluating the equation at each Chebshev point
$\zeta_{k}$, and using (A.5) to convert differentiation into matrix
multiplication. This matrix equation can then be solved by any number of
standard numerical techniques.

For example, the equation

$\left[\frac{\mathrm{d}^{2}}{\mathrm{d}\zeta^{2}}+g(\zeta)\right]\\!f(\zeta)=0\,,$
(A.6)

for any regular function $g(\zeta)$ on
$\zeta\in[{\scalebox{0.75}[1.0]{$-$}}1,1]$ can be approximated by the matrix
equation

$\sum_{k=0}^{N}\left(\mathcal{D}_{2}+\mathcal{G}\right)_{nk}f(\zeta_{k})=0\,,$
(A.7)

where $(\mathcal{G})_{nk}=g(\zeta_{k})\delta_{nk}$,
$(\mathcal{D}_{2})_{nk}=p_{k}^{\prime\prime}(\zeta_{n})$, which is defined in
(A.5), and the number of nodes is chosen to be $N+1$. The solutions to (A.6)
are then well-approximated by the interpolant (A.2) constructed from the null
space of the matrix $(\mathcal{D}_{2}+\mathcal{G})$. As long as $N$ is large
enough, this interpolant then provides an approximate solution to (A.6) whose
pointwise error decays exponentially in $N$.

When the equation we aim to solve is defined on the semi-infinite line
$r\in[0,\infty)$, we must first map it to the interval
$\zeta\in[{\scalebox{0.75}[1.0]{$-$}}1,1]$ to apply Chebyshev interpolation.
We find that the simplest map, the linear fractional transformation

$\zeta(r)=\frac{r-L}{r+L},$ (A.8)

works well, where $L$ is a tuneable parameter that broadly controls which $r$
values are sampled by the corresponding Chebyshev nodes (A.1). In practice, we
take $L$ to be an $\mathcal{O}(1)$ number times the scale over which the
solution varies naturally (set by $\mu^{{\scalebox{0.5}[0.85]{$-$}}1}$ for the
equations we are interested in) though we ensure that our solutions do not
depend on the precise choice of $L$.

For two-dimensional equations of motion like (3.4), we have a choice in how to
apply Chebyshev interpolation. The simplest approach is to expand in Fourier
modes, as in §3.1, to reduce the problem to a system of coupled ordinary
differential equations. That is, we expand each spinor component into their
lowest $2J+1$ Fourier modes

$\psi_{\alpha}(r,\varphi)=\sum_{|\ell|\leq
J}\psi_{\alpha,\ell}(r)\,{{\mathrm{e}}}^{i\ell\varphi}\,,$ (A.9)

such that (3.4) reduces to a finite set of coupled equations (3.6) involving
$2\times(2J+1)$ mode functions. We can then interpolate each radial function
$\psi_{\alpha,\ell}(r)$ with $N+1$ nodes to convert (3.6) into a
$2\times(2J+1)\times(N+1)$-dimensional linear equation. The profiles displayed
in Figure 2 were generated with $J=41$ and $N=61$, though our results are
robust to changing both $J$ and $N$.

Another approach is to interpolate the spinor components $\psi_{\alpha}(x,y)$
using a two-dimensional Cartesian grid. This requires first mapping both $x$
and $y$ from $({\scalebox{0.75}[1.0]{$-$}}\infty,\infty)$ to
$[{\scalebox{0.75}[1.0]{$-$}}1,1]$. We do so via two different maps

$\zeta(x_{i})=\begin{dcases}\arctan\hskip 0.5pt({x_{i}}/{L_{i}})&0\leq
m/\mu\leq 0.8\\\ \tanh\hskip 0.5pt({x_{i}}/{L_{i}})&0.4\leq m/\mu\leq
0.95\end{dcases}\,,$ (A.10)

depending on the value of $m/\mu$, where $x_{i}=\\{x,y\\}$. Here, the
parameters $L_{i}=\\{L_{x},L_{y}\\}$ determine the relevant scales in the
original $x$ or $y$ coordinate. These maps efficiently sample small $(x,y)$
values while avoiding numerical artifacts at large $x$ and $y$. Since the
solutions elongate in the $x$-direction as $m\to\mu$, we used
$L_{x}\simeq\mathcal{O}(1)\times 1/(\mu-m)$ and
$L_{y}\simeq\mathcal{O}(1)\times\mu^{{\scalebox{0.5}[0.85]{$-$}}1}$. We find
that with these ranges of $L_{x}$ and $L_{y}$, our Cartesian results agree
with each other for overlapping $m/\mu$ and that our solutions are independent
of our specific choice. For large numbers of nodes, $N_{x},N_{y}\geq 90$, our
results obtained with a Cartesian grid agree with those obtained using a
Fourier mode expansion, but formulating the problem in the Fourier basis is
significantly faster.

Finally, to compute the moments of our solutions and properly normalize them,
we must numerically evaluate integrals involving these profiles. We do so via
quadrature, in which the integral is approximated by a weighted sum of the
function’s values on the Chebyshev grid. Explicitly, an arbitrary integral of
the density
$|\psi(r,\varphi)|^{2}=|\psi_{0}(r,\varphi)|^{2}+|\psi_{3}(r,\varphi)|^{2}$
can be expressed as

$\displaystyle\int\\!\mathrm{d}^{2}r\,a(r,\varphi)\hskip 0.5pt\hskip
0.5pt|\psi(r,\varphi)|^{2}$
$\displaystyle=\sum_{\alpha,n,\ell,\ell^{\prime}}{}\int_{{\scalebox{0.5}[0.85]{$-$}}1}^{1}\\!\mathrm{d}\zeta\,\int_{0}^{2\pi}\\!\mathrm{d}\varphi\,r^{\prime}(\zeta)a_{n}(\zeta)\hskip
0.5pt{{\mathrm{e}}}^{in\varphi}{{\mathrm{e}}}^{i(\ell-\ell^{\prime})\varphi}\hskip
0.5pt\psi^{\phantom{*}}_{\alpha,\ell}(\zeta)\hskip
0.5pt\psi^{*}_{\alpha,\ell^{\prime}}(\zeta)$ (A.11)
$\displaystyle=2\pi\\!\sum_{\alpha,k,\ell,n}q_{k}\hskip 0.5pt\hskip
0.5ptr^{\prime}(\zeta_{k})\hskip 0.5pta_{n}(\zeta_{k})\hskip
0.5pt\psi^{\phantom{*}}_{\alpha,\ell}(\zeta_{k})\hskip
0.5pt\psi^{*}_{\alpha,\ell+n}(\zeta_{k})\,$

where

$q_{k}=\int_{{\scalebox{0.5}[0.85]{$-$}}1}^{1}\\!\mathrm{d}\zeta\,p_{k}(\zeta)$
(A.12)

is the quadrature weight, which can be efficiently numerically evaluated, and

$a(r,\varphi)=\sum_{|n|\leq
J,k}a_{n}(\zeta_{k})p_{k}(\zeta){{\mathrm{e}}}^{in\varphi}$ (A.13)

is an arbitrary weight function, whose $n$’th Fourier mode has value
$a_{n}(\zeta_{k})$ at $r(\zeta_{k})$. For instance, to normalize the profile
functions in Figure 2, we evaluate (A.11) with $a(r,\varphi)=1$, while to
compute the dipole (5.11) we choose $a(r,\varphi)\propto r\cos\varphi$ and
$a(r,\varphi)\propto r\sin\varphi$.

##  Appendix B Multipole Moments

In this Appendix we provide more details on the derivation of the worldsheet
effective action given in (5.8). We follow the method in [71], adapted to the
case of a string. Our starting point is the term in the action that couples
the gauge field to the four-dimensional current from (5.6), which is pointing
along the axion string,

$S_{\textsc{em}}=-\\!\int\\!\mathrm{d}^{4}x\;j^{\mu}A_{\mu}.$ (B.1)

As in the main text, we take $\mu=0,1,2,3$, $a=0,3$, and $i,j,k_{i}=1,2$.

Our goal is an effective action describing these interactions which is valid
when the wavelength of the gauge field is much larger than the characteristic
size of the source. In this regime, we can Taylor expand the gauge field about
the string source at the origin,

$A_{\mu}(t,\vec{x})=-\\!\sum_{n=0}^{\infty}\frac{1}{n!}\;x^{k_{1}}\\!\cdots
x^{k_{n}}\big{[}\partial_{k_{1}}\\!\cdots\partial_{k_{n}}A_{\mu}\big{]}_{x=(t,0,0,z)}.$
(B.2)

Since our source of interest points along the string (see 5.6), we can drop
the terms with $\mu=1,2$. We are left with the action

$S_{\textsc{em}}=-\\!\sum_{n=0}^{\infty}\int\\!\mathrm{d}^{2}\sigma\\!\int\\!{\rm
d}^{2}r\,\frac{1}{n!}x^{k_{1}}\\!\cdots x^{k_{n}}j^{a}\\!\hskip
0.5pt(t,\vec{x})\big{[}\partial_{k_{1}}\\!\cdots\partial_{k_{n}}A_{a}\big{]}\\!\hskip
0.5pt(t,z).$ (B.3)

It is useful to rewrite the action in terms of gauge invariant quantities. We
first consider the $A^{0}$ term, and note that
$E^{k_{i}}=F^{0k_{i}}=\dot{A}^{k_{i}}+\partial_{k_{i}}A^{0}$. This allows us
to exchange radial derivatives of $A^{0}$ for time-derivatives of $A^{k_{i}}$
and terms involving the electric field, $E$. We find,

$\begin{split}S_{A_{0}}=&-\int\\!\mathrm{d}^{2}\sigma\int\\!\mathrm{d}^{2}r\;j^{0}\big{[}A^{0}\big{]}\\!\hskip
0.5pt(t,z)-\sum_{n=1}^{\infty}\frac{1}{n!}\int\\!\mathrm{d}^{2}\sigma\int\\!{\rm
d}^{2}r\;j^{0}x^{k_{1}}\\!\cdots
x^{k_{n}}\big{[}\partial_{k_{1}}\\!\cdots\partial_{k_{n-1}}E^{k_{n}}\big{]}\\!\hskip
0.5pt(t,z)\\\
&+\sum_{n=1}^{\infty}\frac{1}{n!}\int\\!\mathrm{d}^{2}\sigma\int\\!\mathrm{d}^{2}r\;j^{0}x^{k_{1}}\\!\cdots
x^{k_{n}}\big{[}\partial_{k_{1}}\\!\cdots\partial_{k_{n-1}}\dot{A}^{k_{n}}\big{]}\\!\hskip
0.5pt(t,z)\end{split}$ (B.4)

The last term in this expression can be combined with the $A^{3}$ term of
(B.3) to get a gauge invariant term involving the magnetic field. To see this,
we integrate by parts within the ${\rm d}^{2}\sigma$ integral, and use current
conservation: $\partial_{0}j^{0}=-\partial_{3}j^{3}$. This allows us to
rewrite,

$\begin{split}\left[\int\\!{\rm d}^{2}r\;j^{0}x^{k_{1}}\\!\cdots
x^{k_{n}}\right]\partial_{k_{1}}\\!\cdots\partial_{k_{n-1}}\dot{A}^{k_{n}}&=\left[\int\\!{\rm
d}^{2}r\;(\partial_{3}j^{3})\;x^{k_{1}}\\!\cdots
x^{k_{n}}\right]\partial_{k_{1}}\\!\cdots\partial_{k_{n-1}}A^{k_{n}}\\\
&=-\left[\int\\!{\rm d}^{2}r\;j^{3}x^{k_{1}}\\!\cdots
x^{k_{n}}\right]\partial_{k_{1}}\\!\cdots\partial_{k_{n-1}}\partial_{3}A^{k_{n}},\end{split}$
(B.5)

where in the last equality we have integrated by parts again. With this
replacement, the sum of the last term of (B.4) and the $A^{3}$ term from (B.3)
is

$\begin{split}S_{B}=&-\int\\!\mathrm{d}^{2}\sigma\int\\!{\rm
d}^{2}r\;j^{3}A_{3}-\sum_{n=1}^{\infty}\frac{1}{n!}\int\\!\mathrm{d}^{2}\sigma\int\\!{\rm
d}^{2}r\;j^{3}x^{k_{1}}\\!\cdots
x^{k_{n}}\partial_{k_{1}}\\!\cdots\partial_{k_{n-1}}\big{(}\partial_{3}A^{k_{n}}\big{)}\\\
&-\sum_{n=1}^{\infty}\frac{1}{n!}\int\\!\mathrm{d}^{2}\sigma\int\\!{\rm
d}^{2}r\;j^{3}x^{k_{1}}\\!\cdots
x^{k_{n}}\partial_{k_{1}}\\!\cdots\partial{k_{n-1}}\big{(}\partial_{k_{n}}A_{3}\big{)}\\\
=&-\int\\!\mathrm{d}^{2}\sigma\int\\!{\rm
d}^{2}r\;j^{3}A_{3}-\sum_{n=1}^{\infty}\frac{1}{n!}\int\\!\mathrm{d}^{2}\sigma\int\\!{\rm
d}^{2}r\;j^{3}x^{k_{1}}\\!\cdots
x^{k_{n}}\partial_{k_{1}}\\!\cdots\partial_{k_{n-1}}\epsilon^{3k_{n}j}B_{j},\end{split}$
(B.6)

where we have used
$\partial_{i}A_{3}-\partial_{3}A_{i}=F^{i3}=\epsilon^{3ij}B_{j}$. Adding back
the remaining terms in $S_{A^{0}}$ and relabeling indices, we obtain the full
action

$\begin{split}S_{\textsc{em}}=&-\int\\!\mathrm{d}^{2}\sigma\int\\!{\rm
d}^{2}r\;(j^{0}A_{0}+j^{3}A_{3})-\sum_{n=1}^{\infty}\frac{1}{n!}\int\\!\mathrm{d}^{2}\sigma\left[\int\\!{\rm
d}^{2}rj^{0}x^{k_{1}}\\!\cdots
x^{k_{n}}\right]\partial_{k_{2}}\\!\cdots\partial_{k_{n}}E^{k_{1}}\\\
&-\sum_{n=1}^{\infty}\frac{1}{n!}\int\\!\mathrm{d}^{2}\sigma\left[\int\\!{\rm
d}^{2}r\;j^{3}x^{k_{1}}\\!\cdots
x^{k_{n}}\right]\partial_{k_{2}}\\!\cdots\partial_{k_{n}}\epsilon^{3k_{1}j}B_{j},\end{split}$
(B.7)

which can be written in terms of the field strength as

$\begin{split}S_{\textsc{em}}=&-\int\\!\mathrm{d}^{2}\sigma\int\\!{\rm
d}^{2}r\;(j^{0}A_{0}+j^{3}A_{3})-\sum_{n=1}^{\infty}\frac{1}{n!}\int\\!\mathrm{d}^{2}\sigma\left[\int\\!{\rm
d}^{2}rj^{a}x^{k_{1}}\\!\cdots
x^{k_{n}}\right]\partial_{k_{2}}\\!\cdots\partial_{k_{n}}F_{k_{1}a}.\end{split}$
(B.8)

This action reduces to (5.8) using the current (5.6) and the normalization
condition for $|\mathcal{F}(r,\varphi)|^{2}$.

## References

  * [1] R. D. Peccei and H. R. Quinn, “Constraints Imposed by CP Conservation in the Presence of Instantons,” Phys. Rev. D 16 (1977) 1791–1797.
  * [2] R. D. Peccei and H. R. Quinn, “CP Conservation in the Presence of Instantons,” Phys. Rev. Lett. 38 (1977) 1440–1443.
  * [3] S. Weinberg, “A New Light Boson?,” Phys. Rev. Lett. 40 (1978) 223–226.
  * [4] F. Wilczek, “Problem of Strong $P$ and $T$ Invariance in the Presence of Instantons,” Phys. Rev. Lett. 40 (1978) 279–282.
  * [5] J. Preskill, M. B. Wise, and F. Wilczek, “Cosmology of the Invisible Axion,” Phys. Lett. B 120 (1983) 127–132.
  * [6] M. Dine and W. Fischler, “The Not So Harmless Axion,” Phys. Lett. B 120 (1983) 137–141.
  * [7] L. F. Abbott and P. Sikivie, “A Cosmological Bound on the Invisible Axion,” Phys. Lett. B 120 (1983) 133–136.
  * [8] K. Freese, J. A. Frieman, and A. V. Olinto, “Natural inflation with pseudo - Nambu-Goldstone bosons,” Phys. Rev. Lett. 65 (1990) 3233–3236.
  * [9] S. H.-S. Alexander, M. E. Peskin, and M. M. Sheikh-Jabbari, “Leptogenesis from gravity waves in models of inflation,” Phys. Rev. Lett. 96 (2006) 081301, arXiv:hep-th/0403069.
  * [10] E. Witten, “Some Properties of O(32) Superstrings,” Phys. Lett. B 149 (1984) 351–356.
  * [11] P. Svrcek and E. Witten, “Axions In String Theory,” JHEP 06 (2006) 051, arXiv:hep-th/0605206.
  * [12] A. Arvanitaki, S. Dimopoulos, S. Dubovsky, N. Kaloper, and J. March-Russell, “String Axiverse,” Phys. Rev. D 81 (2010) 123530, arXiv:0905.4720 [hep-th].
  * [13] M. Demirtas, C. Long, L. McAllister, and M. Stillman, “The Kreuzer-Skarke Axiverse,” JHEP 04 (2020) 138, arXiv:1808.01282 [hep-th].
  * [14] M. Demirtas, N. Gendler, C. Long, L. McAllister, and J. Moritz, “PQ axiverse,” JHEP 06 (2023) 092, arXiv:2112.04503 [hep-th].
  * [15] J. Stout, “Instanton expansions and phase transitions,” JHEP 05 (2022) 168, arXiv:2012.11605 [hep-th].
  * [16] J. Fan, K. Fraser, M. Reece, and J. Stout, “Axion Mass from Magnetic Monopole Loops,” Phys. Rev. Lett. 127 no. 13, (2021) 131602, arXiv:2105.09950 [hep-ph].
  * [17] A. Vilenkin and A. E. Everett, “Cosmic Strings and Domain Walls in Models with Goldstone and PseudoGoldstone Bosons,” Phys. Rev. Lett. 48 (1982) 1867–1870.
  * [18] P. Sikivie, “Of Axions, Domain Walls and the Early Universe,” Phys. Rev. Lett. 48 (1982) 1156–1159.
  * [19] R. L. Davis, “Cosmic Axions from Cosmic Strings,” Phys. Lett. B 180 (1986) 225–230.
  * [20] D. Harari and P. Sikivie, “On the Evolution of Global Strings in the Early Universe,” Phys. Lett. B 195 (1987) 361–365.
  * [21] E. P. S. Shellard, “Cosmic String Interactions,” Nucl. Phys. B 283 (1987) 624–656.
  * [22] R. L. Davis and E. P. S. Shellard, “Do Axions Need Inflation?,” Nucl. Phys. B 324 (1989) 167–186.
  * [23] C. Hagmann and P. Sikivie, “Computer simulations of the motion and decay of global strings,” Nucl. Phys. B 363 (1991) 247–280.
  * [24] R. A. Battye and E. P. S. Shellard, “Global string radiation,” Nucl. Phys. B 423 (1994) 260–304, arXiv:astro-ph/9311017.
  * [25] R. A. Battye and E. P. S. Shellard, “Axion string constraints,” Phys. Rev. Lett. 73 (1994) 2954–2957, arXiv:astro-ph/9403018. [Erratum: Phys.Rev.Lett. 76, 2203–2204 (1996)].
  * [26] M. Yamaguchi, M. Kawasaki, and J. Yokoyama, “Evolution of axionic strings and spectrum of axions radiated from them,” Phys. Rev. Lett. 82 (1999) 4578–4581, arXiv:hep-ph/9811311.
  * [27] V. B. . Klaer and G. D. Moore, “The dark-matter axion mass,” JCAP 11 (2017) 049, arXiv:1708.07521 [hep-ph].
  * [28] M. Gorghetto, E. Hardy, and G. Villadoro, “Axions from Strings: the Attractive Solution,” JHEP 07 (2018) 151, arXiv:1806.04677 [hep-ph].
  * [29] A. Vaquero, J. Redondo, and J. Stadler, “Early seeds of axion miniclusters,” JCAP 04 (2019) 012, arXiv:1809.09241 [astro-ph.CO].
  * [30] M. Buschmann, J. W. Foster, and B. R. Safdi, “Early-Universe Simulations of the Cosmological Axion,” Phys. Rev. Lett. 124 no. 16, (2020) 161103, arXiv:1906.00967 [astro-ph.CO].
  * [31] M. Gorghetto, E. Hardy, and G. Villadoro, “More axions from strings,” SciPost Phys. 10 no. 2, (2021) 050, arXiv:2007.04990 [hep-ph].
  * [32] M. Dine, N. Fernandez, A. Ghalsasi, and H. H. Patel, “Comments on axions, domain walls, and cosmic strings,” JCAP 11 (2021) 041, arXiv:2012.13065 [hep-ph].
  * [33] M. Buschmann, J. W. Foster, A. Hook, A. Peterson, D. E. Willcox, W. Zhang, and B. R. Safdi, “Dark matter from axion strings with adaptive mesh refinement,” Nature Commun. 13 no. 1, (2022) 1049, arXiv:2108.05368 [hep-ph].
  * [34] P. Agrawal, A. Hook, and J. Huang, “A CMB Millikan experiment with cosmic axiverse strings,” JHEP 07 (2020) 138, arXiv:1912.02823 [astro-ph.CO].
  * [35] J. N. Benabou, M. Buschmann, S. Kumar, Y. Park, and B. R. Safdi, “Signatures of Primordial Energy Injection from Axion Strings,” arXiv:2308.01334 [hep-ph].
  * [36] C. G. Callan, Jr. and J. A. Harvey, “Anomalies and Fermion Zero Modes on Strings and Domain Walls,” Nucl. Phys. B 250 (1985) 427–436.
  * [37] D. B. Kaplan and A. Manohar, “Anomalous Vortices and Electromagnetism,” Nucl. Phys. B 302 (1988) 280–290.
  * [38] S. G. Naculich, “Axionic Strings: Covariant Anomalies and Bosonization of Chiral Zero Modes,” Nucl. Phys. B 296 (1988) 837–867.
  * [39] A. Manohar, “Anomalous Vortices and Electromagnetism. 2.,” Phys. Lett. B 206 (1988) 276. [Erratum: Phys.Lett. 209, 543 (1988)].
  * [40] J. A. Harvey and S. G. Naculich, “Cosmic Strings From Pseudoanomalous U(1)$s$,” Phys. Lett. B 217 (1989) 231–237.
  * [41] D. Harari and P. Sikivie, “Effects of a Nambu-Goldstone boson on the polarization of radio galaxies and the cosmic microwave background,” Phys. Lett. B 289 (1992) 67–72.
  * [42] J. D. Blum and J. A. Harvey, “Anomaly inflow for gauge defects,” Nucl. Phys. B 416 (1994) 119–136, arXiv:hep-th/9310035.
  * [43] J. A. Harvey and O. Ruchayskiy, “The Local structure of anomaly inflow,” JHEP 06 (2001) 044, arXiv:hep-th/0007037.
  * [44] B. Heidenreich, M. Reece, and T. Rudelius, “The Weak Gravity Conjecture and axion strings,” JHEP 11 (2021) 004, arXiv:2108.11383 [hep-th].
  * [45] H. Fukuda and K. Yonekura, “Witten effect, anomaly inflow, and charge teleportation,” JHEP 01 (2021) 119, arXiv:2010.02221 [hep-th].
  * [46] H. Fukuda, A. V. Manohar, H. Murayama, and O. Telem, “Axion strings are superconducting,” JHEP 06 (2021) 052, arXiv:2010.02763 [hep-ph].
  * [47] M. Ibe, S. Kobayashi, Y. Nakayama, and S. Shirai, “On Stability of Fermionic Superconducting Current in Cosmic String,” JHEP 05 (2021) 217, arXiv:2102.05412 [hep-ph].
  * [48] P. Agrawal, A. Hook, J. Huang, and G. Marques-Tavares, “Axion string signatures: a cosmological plasma collider,” JHEP 01 (2022) 103, arXiv:2010.15848 [hep-ph].
  * [49] E. Witten, “Superconducting Strings,” Nucl. Phys. B 249 (1985) 557–592.
  * [50] M. Dine, W. Fischler, and M. Srednicki, “A Simple Solution to the Strong CP Problem with a Harmless Axion,” Phys. Lett. B 104 (1981) 199–202.
  * [51] A. R. Zhitnitsky, “On Possible Suppression of the Axion Hadron Interactions. (In Russian),” Sov. J. Nucl. Phys. 31 (1980) 260.
  * [52] Y. Abe, Y. Hamada, and K. Yoshioka, “Electroweak axion string and superconductivity,” JHEP 06 (2021) 172, arXiv:2010.02834 [hep-ph].
  * [53] C. T. Hill and L. M. Widrow, “Superconducting Cosmic Strings with Massive Fermions,” Phys. Lett. B 189 (1987) 17–22.
  * [54] H. A. Weldon, “Effective Fermion Masses of Order gT in High Temperature Gauge Theories with Exact Chiral Invariance,” Phys. Rev. D 26 (1982) 2789.
  * [55] A. Achucarro and T. Vachaspati, “Semilocal and electroweak strings,” Phys. Rept. 327 (2000) 347–426, arXiv:hep-ph/9904229.
  * [56] R. Jackiw and P. Rossi, “Zero Modes of the Vortex - Fermion System,” Nucl. Phys. B 190 (1981) 681–691.
  * [57] G. D. Starkman, D. Stojkovic, and T. Vachaspati, “Neutrino zero modes on electroweak strings,” Phys. Rev. D 63 (2001) 085011, arXiv:hep-ph/0007071.
  * [58] G. Starkman, D. Stojkovic, and T. Vachaspati, “Zero modes of fermions with a general mass matrix,” Phys. Rev. D 65 (2002) 065003, arXiv:hep-th/0103039.
  * [59] S. G. Naculich, “Fermions destabilize electroweak strings,” Phys. Rev. Lett. 75 (1995) 998–1001, arXiv:hep-ph/9501388.
  * [60] H. Liu and T. Vachaspati, “Perturbed electroweak strings and fermion zero modes,” Nucl. Phys. B 470 (1996) 176–194, arXiv:hep-ph/9511216.
  * [61] S. R. Coleman and E. J. Weinberg, “Radiative Corrections as the Origin of Spontaneous Symmetry Breaking,” Phys. Rev. D 7 (1973) 1888–1910.
  * [62] C. Bender and S. Orszag, Advanced Mathematical Methods for Scientists and Engineers I: Asymptotic Methods and Perturbation Theory.  Springer, 1999.
  * [63] D. B. Kaplan, “Chiral fermions on the lattice,” Nucl. Phys. B Proc. Suppl. 30 (1993) 597–600.
  * [64] K. Aitken, A. Baumgartner, A. Karch, and B. Robinson, “3d Abelian Dualities with Boundaries,” JHEP 03 (2018) 053, arXiv:1712.02801 [hep-th].
  * [65] E. H. Fradkin, Field Theories of Condensed Matter Physics, vol. 82.  Cambridge Univ. Press, Cambridge, UK, 2, 2013.
  * [66] F. D. M. Haldane, “Model for a Quantum Hall Effect without Landau Levels: Condensed-Matter Realization of the ’Parity Anomaly’,” Phys. Rev. Lett. 61 (1988) 2015–2018.
  * [67] W. P. Su, J. R. Schrieffer, and A. J. Heeger, “Solitons in polyacetylene,” Phys. Rev. Lett. 42 (1979) 1698–1701.
  * [68] R. Jackiw and J. R. Schrieffer, “Solitons with Fermion Number 1/2 in Condensed Matter and Relativistic Field Theories,” Nucl. Phys. B 190 (1981) 253–265.
  * [69] K. Jansen and M. Schmaltz, “Critical momenta of lattice chiral fermions,” Phys. Lett. B 296 (1992) 374–378, arXiv:hep-lat/9209002.
  * [70] J. Goldstone and F. Wilczek, “Fractional Quantum Numbers on Solitons,” Phys. Rev. Lett. 47 (1981) 986–989.
  * [71] A. Ross, “Multipole expansion at the level of the action,” Phys. Rev. D 85 (2012) 125033, arXiv:1202.4750 [gr-qc].
  * [72] C. G. Callan, Jr., “Disappearing Dyons,” Phys. Rev. D 25 (1982) 2141.
  * [73] L. Trefethen, Approximation Theory and Approximation Practice, vol. 128.  Siam, 2013.
  * [74] J.-P. Berrut and L. Trefethen, “Barycentric Lagrange Interpolation,” SIAM review 46 (2004) 501–517.
  * [75] N. Higham, “The Numerical Stability of Barycentric Lagrange Interpolation,” IMA Journal of Numerical Analysis 24 (2004) 547–556.
  * [76] R. Baltensperger, “Improving the Accuracy of the Matrix Differentiation Method for Arbitrary Collocation Points,” Applied Numerical Mathematics 33 (2000) 143–149.