††thanks: Orcid ID https://orcid.org/0000-0002-6229-5620 ††thanks: Orcid ID
https://orcid.org/0000-0002-0209-7358

# Equilibrium parameters in coupled storage ring lattices and practical
applications

V. Ziemann<EMAIL_ADDRESS>Department of Physics and Astronomy,
Uppsala University, Uppsala, Sweden  A. Streun<EMAIL_ADDRESS>Paul
Scherrer Institut, CH-5232 Villigen, Switzerland

###### Abstract

We calculate equilibrium emittances and damping times due to the emission of
synchrotron radiation for coupled storage ring lattices by evaluating the
projections of the commonly used synchrotron radiation integrals onto the
normal modes of the coupled motion. Orbit distortion is included by
calculating off-axis contributions to the radiation integrals. We provide
explicit formulae for fast forward calculation, which have been implemented
into the interactive lattice design code OPA.

###### pacs:

29.20.Dh, 29.27.Bd, 41.60.Ap

##  I Introduction

The emittance is one of the crucial quantities that determines the transverse
beam dimensions and thus the performance of storage rings. In electron rings
it is determined by the balance of damping and excitation due to the emission
of synchrotron radiation Sands (1969). For rings with negligible coupling
between the transverse planes, the effect is quantified by the classical
synchrotron radiation integrals Helm et al. (1973); Edwards and Sypers (1993).
Coupling was initially described by a heuristic emittance coupling constant
$\kappa=\varepsilon_{y}/\varepsilon_{x}$ where $\varepsilon_{x}$ and
$\varepsilon_{y}$ are the horizontal and vertical emittance, respectively.
With the increased performance of accelerators elaborate methods to calculate
the emittance coupling were developed Chao (1979a); Ohmi et al. (1994); Wolski
(2006). We complement these often elegant, but sometimes less intuitive
methods, by a generalization of the radiation integrals to a coupled lattice.

The beam optics code OPA OPA , written and maintained by one of the authors,
is tailored to interactively design beam optical systems. Especially for the
initial design phase, estimating the equilibrium beam properties from the
optics of a non-periodic system, helps to guide the optimization; and this
requires the evaluation of the radiation integrals, especially for
intrinsically coupled lattices, such as Spiral-COSAMI Wrulich et al. (2021).

We consider the normal mode parametrization of the coupled one-turn transfer
matrix as introduced by Edward and Teng Edwards and Teng (1973) and extended
by Sagan and Rubin Sagan and Rubin (1999), and calculate the projection of the
radiation integrals onto the corresponding normal modes. Note that we also
implemented the algorithm in MATLAB using the software from Ziemann (2019)
which is useful to validate and cross-check the software.

In the discussion and also with regard to practical implementation, we assume
that all linear beam optical elements such as quadrupoles or dipoles are
defined by their transfer matrices for the un-coupled elements and sandwiched
between coordinate rotation matrices that describe their misalignment
pertaining to coupling of the transverse planes (i.e. a skew quadrupole is
represented by a quadrupole sandwiched between $\pm 45\mbox{${}^{\circ}$}$
rotations).

We only consider transverse coupling and assume the longitudinal dynamics to
be decoupled by adiabatic approximation, i.e. the synchrotron oscillation are
considered slow compared to the betatron oscillations, and a particle’s
longitudinal momentum is assumed as a constant. Furthermore our analysis is
based on highly relativistic, paraxial and large bending radius approximations
as commonly used for high energy storage rings, where beam properties are
determined by synchrotron radiation effects.

We start with a brief reminder of the required beam dynamics formalism, then
we calculate the radiation integral $I_{5}$ (in the notation of ref Helm et
al. (1973)) responsible for excitation of betatron oscillations followed by
the evaluation of $I_{4}$ that describes damping. From these follow the
equilibrium beam parameters of a coupled lattice. We proceed further including
contributions from orbit excursions, which may be due to energy offsets or
lattice imperfections. Finally we discuss implementation issues and present
some applications.

##  II Normal modes and dispersions

Transverse beam dynamics using coordinates $(x,x^{\prime},y,y^{\prime},\Delta
E/E)$ is covered by $5\times 5$ transfer matrices for elements and [circular]
concatenations of elements:

$\left(\begin{array}[]{c|c}R&\vec{V}\\\ 0&1\end{array}\right)$ (1)

with $R$ the $4\times 4$ transverse transfer matrix and $\vec{V}$ the
4-dimensional vector of dispersion production.

If $R$ is the one-turn matrix of a storage ring, starting at an arbitrary
longitudinal position, a normal mode decomposition may be found Edwards and
Teng (1973); Sagan and Rubin (1999):

$R=T^{-1}{\cal A}^{-1}{\cal O}{\cal A}T$ (2)

where $\cal O$ is the matrix containing the eigen-tunes

${\cal O}=\left(\begin{array}[]{cc}O_{a}&0\\\
0&O_{b}\end{array}\right)\quad\mathrm{with}\quad
O_{a,b}=\left(\begin{array}[]{rc}\cos\mu_{a,b}&\sin\mu_{a,b}\\\
-\sin\mu_{a,b}&\cos\mu_{a,b}\end{array}\right)$ (3)

and ${\cal A}$ is the matrix that contains the local normal-mode beta
functions. It is given by

${\cal A}=\left(\begin{array}[]{cc}{\cal A}_{a}&0\\\ 0&{\cal
A}_{b}\end{array}\right)\quad\mathrm{with}\quad{\cal
A}_{a,b}=\left(\begin{array}[]{cc}\frac{1}{\sqrt{\beta_{a,b}}}&0\\\
\frac{\alpha_{a,b}}{\sqrt{\beta_{a,b}}}&\sqrt{\beta_{a,b}}\end{array}\right)$
(4)

where the indices $a$ and $b$ label the two eigen modes. The coupling matrix
$T$ and its inverse $T^{-1}$ are given by Sagan and Rubin (1999)

$T=\left(\begin{array}[]{cc}gI&-C\\\
C^{+}&gI\end{array}\right)\qquad\mathrm{and}\qquad
T^{-1}=\left(\begin{array}[]{cc}gI&C\\\ -C^{+}&gI\end{array}\right),$ (5)

with the $2\times 2$ identity matrix $I$ and the $2\times 2$ coupling matrix
$C$. Its symplectic conjugate $C^{+}$ and the scalar $g$ are given by

$C^{+}=C^{-1}\det C,\quad\mbox{\rm and}\quad g^{2}=1-\det C.$ (6)

The periodic solution for the four-dimensional dispersion
$\vec{D}=(\vec{D}_{x},\vec{D}_{y})^{T}=(D_{x},D^{\prime}_{x},D_{y},D^{\prime}_{y})^{T}$
is constrained by

$\vec{D}=R\vec{D}+\vec{V}$ (7)

which can be solved for $\vec{D}$ with the result

$\vec{D}=(1-R)^{-1}\vec{V}\ .$ (8)

Thus $\vec{D}$ describes the periodic dispersion in physical space. The normal
mode dispersion is given by

$\vec{\cal D}=(\vec{\cal D}_{a},\vec{\cal D}_{b})^{T}=({\cal D}_{a},{\cal
D}_{a}^{\prime},{\cal D}_{b},{\cal D}_{b}^{\prime})^{T}=T\vec{D}.$ (9)

Further application of matrix $\cal{A}$ gives the normal mode dispersion in
normalized phase space:

$\vec{\tilde{D}}=(\tilde{D}_{1},\tilde{D}_{2},\tilde{D}_{3},\tilde{D}_{4})^{T}={\cal
A}\vec{\cal D}=S^{-1}\vec{D}$ (10) $\mbox{\rm with~{}~{}~{}~{}}S=({\cal
A}T)^{-1}=T^{-1}{\cal A}^{-1}.$ (11)

We thus defined dispersion in three different coordinate systems, which we
will use as it is convenient, and which are related as

$\vec{\tilde{D}}_{\mbox{\rm~{}(normalized)}}\quad\stackrel{{\scriptstyle\cal
A}}{{\longleftarrow}}\quad\vec{\cal
D}_{\mbox{\rm~{}(decoupled)}}\quad\stackrel{{\scriptstyle
T}}{{\longleftarrow}}\quad\vec{D}_{\mbox{\rm~{}(real space)}}$

Having laid out the theoretical framework and definition of notations we
proceed to calculate the effect of quantum excitations on the normal modes.

##  III Quantum excitation

The stochastic nature of the emission of synchrotron radiation causes heating
of the beam if the photons are emitted at locations with dispersion. In a
planar un-coupled lattice the effect of the lattice parameters, such as beta
function $\beta$ and dispersion $D$ is described by

${\cal H}=\beta_{x}D_{x}^{\prime 2}+2\alpha_{x}D_{x}D_{x}^{\prime}+\gamma
D_{x}^{2}$ (12)

with $\gamma=(1+\alpha^{2})/\beta.$ The synchrotron radiation integral $I_{5}$
is closely related to ${\cal H}$ by

$I_{5}=\oint{\cal H}|h|^{3}\,ds\ ,$ (13)

where $h=1/\rho$ is the local curvature in the bending magnets, and the
integral extends over one turn in the ring as discussed in chapter 7 of
Edwards and Sypers (1993). We follow the same logic to generalize the
derivation to the coupled case and assume to know the dispersion in real space
coordinates, $\vec{D}$ at every point in the ring and especially inside the
dipoles, where we know that the emission of photons happens. This emission is
characterized by an average loss, which causes damping and a fluctuating part
$u$ that has average energy value zero $\langle u\rangle=0$ and squared
expectation value $u^{2}.$ Here $u=\Delta E/E$ is the relative energy
deviation of a particle. We consider only one location and in a particular
realization of the random process the energy loss at a given turn $i$ will be
denoted by $u_{i}.$ We will use the statistical properties later on. In order
to understand the excitation process we consider a single electron that
travels on the closed orbit appropriate for its energy. After the emission
process, the electron is still at the same place, but the orbit appropriate
for its new energy is different, because of the finite dispersion at the
location of emission. Consequently the electron will start oscillation around
the new equilibrium orbit. The change of coordinates $\Delta\vec{x}$ at which
it starts the oscillation are given by

$\Delta\vec{x}=-u_{i}\vec{D}$ (14)

and after traveling $n$ turns with the chance of an emission process on each
turn we expect the electron’s position to be

$\vec{x}_{n}=-\sum_{i=1}^{n}R^{n-i}u_{i}\vec{D}\ .$ (15)

Using the parametrization from eq. 2 it is obvious that powers of the transfer
matrix $R^{m}$ can be evaluated by

$R^{m}=T^{-1}{\cal A}^{-1}{\cal O}^{m}{\cal A}T$ (16)

which can be interpreted in the following way. Reading from right to left: $T$
removes the coupling, ${\cal A}$ turns the ellipses into circles, and then
${\cal O}$ takes care of the phase advance or tune for the specified number of
turns $m.$ Finally we move back from normalized phase space to real space by
the inverses of ${\cal A}$ and $T.$

Returning to eq. 15 we simplify it by mapping the entire equation into
normalized phase space by left multiplying with ${\cal A}T$. We use the
notation $\vec{\tilde{x}}={\cal A}T\vec{x}$ and write

$\vec{\tilde{x}}_{n}=-{\cal
A}T\left(\sum_{i=1}^{n}R^{n-i}u_{i}\right)\vec{D}=-\left(\sum_{i=1}^{n}{\cal
O}^{n-i}u_{i}\right){\cal A}T\vec{D}=-\left(\sum_{i=1}^{n}{\cal
O}^{n-i}u_{i}\right)\vec{\tilde{D}}$ (17)

with $\vec{\tilde{D}}$ the dispersion vector mapped into normalized phase
space from eq. 10. Equation 17 has a nice intuitive interpretation. Instead of
every energy loss $u_{i}$ producing a kick of magnitude $u_{i}\vec{D}$ in real
space, it produces a kick of magnitude $u_{i}\vec{\tilde{D}}$ in normalized
phase space that propagates by multiplying it with the rotation matrix ${\cal
O}$ for a given number of turns and eventually all the kicks from the separate
turns are summed up.

In order to calculate the emittance growth due to such a sequence of kicks we
first need to specify which emittance we really mean. In the following we will
use the emittances of the two normal modes. In the case of an uncoupled
lattice this will revert to the common definition of emittances. Moreover, the
emittances are the ensemble averages over the Courant-Snyder action variables
$J_{a}$ and $J_{b}$ for the respective normal modes, labeled $a$ and $b.$ We
start by writing the previous equation in $2\times 2$ block matrix form

$\left(\begin{array}[]{c}\vec{\tilde{x}}_{a,n}\\\
\vec{\tilde{x}}_{b,n}\end{array}\right)=-\sum_{i=1}^{n}u_{i}\left(\begin{array}[]{cc}O_{a}^{n-i}&0\\\
0&O_{b}^{n-i}\end{array}\right)\left(\begin{array}[]{c}\vec{\tilde{D}}_{a}\\\
\vec{\tilde{D}}_{b}\end{array}\right)=-\sum_{i=1}^{n}u_{i}\left(\begin{array}[]{c}O_{a}^{n-i}\vec{\tilde{D}}_{a}\\\
O_{b}^{n-i}\vec{\tilde{D}}_{b}\end{array}\right)$ (18)

where $\vec{\tilde{x}}=(\vec{\tilde{x}}_{a},\vec{\tilde{x}}_{b})^{T}.$ Here
$\vec{\tilde{x}}_{a}$ are just the first two components of the vector
$\vec{\tilde{x}}$ and $\vec{\tilde{x}}_{b}$ components three and four. The
action variable for the first normal-mode after $n$ turns is thus given by

$J_{a,n}=\frac{1}{2}\vec{\tilde{x}}_{a,n}^{T}\vec{\tilde{x}}_{a,n}\ .$ (19)

We now insert $\vec{\tilde{x}}_{a,n}$ from eq. 18 and arrive at

$\displaystyle J_{a,n}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\sum_{i=1}^{n}\sum_{j=1}^{n}u_{i}u_{j}\vec{\tilde{D}}_{a}^{T}\left(O_{a}^{n-j}\right)^{T}O_{a}^{n-i}\vec{\tilde{D}}_{a}$
$\displaystyle=$
$\displaystyle\frac{1}{2}\sum_{i=1}^{n}\sum_{j=1}^{n}u_{i}u_{j}\vec{\tilde{D}}_{a}^{T}\left(O_{a}^{T}\right)^{n-j}O_{a}^{n-i}\vec{\tilde{D}}_{a}$
$\displaystyle=$
$\displaystyle\frac{1}{2}\vec{\tilde{D}}_{a}^{T}\sum_{i=1}^{n}\sum_{j=1}^{n}O_{a}^{j-i}u_{i}u_{j}\vec{\tilde{D}}_{a}$

where we use that the power of a transpose matrix equals the transpose of the
matrix to the same power. Furthermore, the matrix $O_{a}$ is a rotation matrix
and therefore orthogonal. This implies that its transpose equals its inverse.

Now we consider the ensemble average over many particles. Since the emission
of photons from turn to turn is uncorrelated and has rms magnitude $u^{2}$ as
discussed before we have $\langle u_{i}u_{j}\rangle=u^{2}\delta_{ij}$ which
allows us to evaluate one sum and the other one sums over a constant $u^{2}.$
We obtain

$\langle
J_{a,n}\rangle=\frac{1}{2}\vec{\tilde{D}}_{a}^{T}\sum_{i=1}^{n}\sum_{j=1}^{n}u^{2}\delta_{ij}O^{j-i}_{a}\vec{\tilde{D}}_{a}=\frac{1}{2}\vec{\tilde{D}}_{a}^{T}\vec{\tilde{D}}_{a}nu^{2}$
(21)

and observe that the average action variable grows linearly with the number of
turns $n$.

The magnitude of the term $u^{2}$ is the rms energy kick received by the beam.
It is then given by the ratio of the rms spread of the emitted photons
$\dot{N}\langle\varepsilon^{2}\rangle$ and the beam energy $E$ and is given by
(Section 3.1.4 in Zimmermann et al. (2013))

$u^{2}=\dot{N}\langle\varepsilon^{2}\rangle=\frac{55}{24\sqrt{3}}\frac{\varepsilon_{c}P_{\gamma}}{E^{2}}=\frac{4cr_{e}}{3}C_{q}\frac{\gamma^{5}}{|\rho^{3}|}$
(22)

with $r_{e}$ being the classical electron radius, the constant $C_{q}$, and
the critical photon energy $\varepsilon_{c}$ given by

$C_{q}=\frac{55\,\hbar
c}{32\sqrt{3}\,mc^{2}}\qquad\mathrm{and}\qquad\varepsilon_{c}=\frac{3\,\hbar
c\gamma^{3}}{2\,|\rho|}$ (23)

The radiated power is given by

$P_{\gamma}=\frac{e^{2}c^{3}}{2\pi}C_{\gamma}E^{2}B^{2}=\frac{cC_{\gamma}}{2\pi}\frac{E^{4}}{\rho^{2}}\qquad\mbox{\rm
with}\quad C_{\gamma}=\frac{4\pi r_{e}}{3(mc^{2})^{3}}\ .$ (24)

Integrating the action variable $J_{a}$ over the circumference of the ring we
collect the contributions of all the dipole magnets in the ring. This allows
us to define a growth rate for the normal-mode action variable

$\frac{d\langle
J_{a}\rangle}{dt}=\frac{2}{3}C_{q}r_{e}\gamma^{5}\oint\frac{\vec{\tilde{D}}_{a}^{T}\vec{\tilde{D}}_{a}}{|\rho^{3}|}ds=\frac{2}{3}C_{q}r_{e}\gamma^{5}I_{5a}$
(25)

and similarly for the other normal mode $b.$ Here we implicitly define the
radiation integral $I_{5a}$ for the first normal mode. Comparing with the
expressions in refs. Edwards and Sypers (1993); Zimmermann et al. (2013); Lee
(2004) we observe that ${\cal
H}_{a}=\vec{\tilde{D}}_{a}^{T}\vec{\tilde{D}}_{a}$ takes the place of the
commonly used symbol ${\cal H}$ that occurs in the description of emittance
growth in uncoupled rings. It is important to note, that only the absolute
curvature of the bending field is used but not the direction of deflection,
because it provides the energy fluctuations as the source of emittance, and
the normal mode dispersions project this noise on to the transverse eigen-
modes.

For an uncoupled lattice the coupling matrix $T$ from eq. 5 is unity and the
matrix ${\cal A}$ contains the conventional un-coupled beta functions. Using
the first two components of eq. 10 we can write

$\displaystyle\left(\begin{array}[]{c}\tilde{D}_{x}\\\
\tilde{D}^{\prime}_{x}\end{array}\right)$ $\displaystyle=$ $\displaystyle{\cal
A}_{x}\left(\begin{array}[]{c}D_{x}\\\
D^{\prime}_{x}\end{array}\right)=\left(\begin{array}[]{cc}1/\sqrt{\beta_{x}}&0\\\
\alpha_{x}/\sqrt{\beta_{x}}&\sqrt{\beta_{x}}\end{array}\right)\left(\begin{array}[]{c}D_{x}\\\
D^{\prime}_{x}\end{array}\right)$ (34) $\displaystyle=$
$\displaystyle\left(\begin{array}[]{c}D_{x}/\sqrt{\beta_{x}}\\\
\alpha_{x}D_{x}/\sqrt{\beta_{x}}+\sqrt{\beta_{x}}D^{\prime}_{x}\end{array}\right)$
(37)

and multiplying with its transpose, ${\cal
H}_{x}=\vec{\tilde{D}}_{x}^{T}\vec{\tilde{D}}_{x}$, yields eq. 12.

##  IV Damping

The second consequence of the emission of synchrotron radiation is the damping
of transverse oscillations. In dipole magnets synchrotron radiation is emitted
and that emission is energy dependent with the recoil from the photons along
the direction of motion of the electrons. Subsequently the energy is restored
in a radio-frequency cavity, but only the longitudinal component of the
momentum vector is increased. The joint effect of emitting and restoring the
energy leads to overall damping. We thus proceed to determine the damping
times of the two normal modes $a$ and $b$ and loosely follow the description
from ref. Lee (2004).

Figure 1: Bending magnet geometry

In dipoles we denote the average emitted energy by $\Delta\bar{u}$ which is
related to the power of the emitted synchrotron radiation $P_{\gamma}$ by
$\Delta\bar{u}=P_{\gamma}dl/cE$ where $dl/c$ is a segment of the arc an
electron travels in a magnet, see fig. 1, and $P_{\gamma}$ is given in eq. 24.

If an electron suddenly changes its energy by an amount $\Delta\bar{u}$ due to
photon emission, at a location with dispersion $\vec{D}$, its reference orbit
will jump away from the electron and it will perform betatron oscillations
starting with $\Delta\vec{x}=-\vec{D}\Delta\bar{u}$ around the new reference
orbit. Multiplying both $\Delta\vec{x}$ and the dispersion $\vec{D}$ by ${\cal
A}T$ we obtain the effect on the normal modes

$\Delta\vec{\tilde{x}}=-\vec{\tilde{D}}\Delta\bar{u}$ (38)

and the action variables $J_{a}$ and $J_{b}$ will change by

$\displaystyle\Delta J_{a}$ $\displaystyle=$
$\displaystyle\tilde{x}_{1}\Delta\tilde{x}_{1}+\tilde{x}_{2}\Delta\tilde{x}_{2}=-\left(\tilde{x}_{1}\tilde{D}_{1}+\tilde{x}_{2}\tilde{D}_{2}\right)\Delta\bar{u},$
$\displaystyle\Delta J_{b}$ $\displaystyle=$
$\displaystyle\tilde{x}_{3}\Delta\tilde{x}_{3}+\tilde{x}_{4}\Delta\tilde{x}_{4}=-\left(\tilde{x}_{3}\tilde{D}_{3}+\tilde{x}_{4}\tilde{D}_{4}\right)\Delta\bar{u}.$
(39)

We now need to determine the dependence of the energy loss $\Delta\bar{u}$ on
the position in a combined function magnet with position-dependent magnetic
field. We only consider an upright combined function magnet with only a
horizontal dependence of the vertical field component. Other orientations or
rotated magnets can easily be accommodated by coordinate rotations which enter
the analysis by their influence on the coupling matrix $T.$ Thus, without loss
of generality, we describe the total average energy of the emitted photons by
Lee (2004)

$\displaystyle\Delta\bar{u}$ $\displaystyle=$
$\displaystyle-\frac{P_{\gamma}}{E}\frac{dl}{c}=-\frac{P_{0}}{cE}\left(1+\frac{2}{B}\frac{\partial
B}{\partial x}x\right)\left(1+\frac{x}{\rho}\right)ds$ (40)
$\displaystyle\approx$
$\displaystyle-\frac{P_{0}}{cE}\left(1+\frac{2}{B}\frac{\partial B}{\partial
x}x+\frac{x}{\rho}\right)ds$ (41)

where $x$ denotes the horizontal coordinate in the lab system, and $P_{0}$ is
the power radiated on the design orbit. The term with the derivative of the
magnetic field arises from the quadratic dependence of the emitted power on
the magnetic field, see eq. 24, and the second term with $x/\rho$ describes
the longer path of the electron in the magnet if it is further out with larger
$x$, see fig. 1. In dipoles containing also a skew gradient we would need to
add $(dBy/dy)y$ to $(dBy/dx)x$ and consistently carry through all following
steps.

Expressing the horizontal coordinate $x$ through the normalized normal mode
coordinates $\tilde{x}$ by using the matrix $S$ from eq. 11,

$x=\sum_{i=1}^{4}S_{1i}\tilde{x}_{i}$ (42)

allows us to write the change in the action variables as

$\displaystyle\Delta J_{a}$ $\displaystyle=$
$\displaystyle\frac{P_{0}}{cE}\left(1+\sum_{i=1}^{4}\left(\tilde{x}_{1}\tilde{D}_{1}+\tilde{x}_{2}\tilde{D}_{2}\right)S_{1i}\left(\frac{2}{B}\frac{\partial
B}{\partial x}\tilde{x}_{i}+\frac{\tilde{x}_{i}}{\rho}\right)\right)ds.$ (43)

Averaging over the normal-mode phases denoted by angle brackets, we find

$\displaystyle\langle\Delta J_{a}\rangle$ $\displaystyle=$
$\displaystyle\frac{P_{0}ds}{cE}\left(\frac{2}{B}\frac{\partial B}{\partial
x}+\frac{1}{\rho}\right)\left(S_{11}\langle\tilde{x}_{1}^{2}\rangle\tilde{D}_{1}+S_{12}\langle\tilde{x}_{2}^{2}\rangle\tilde{D}_{2}\right)\,,$
(44)

where we use the fact that $\langle\tilde{x}_{i}\rangle=0$ and
$\langle\tilde{x}_{i}\tilde{x}_{j}\rangle=0$ with $i\neq j$, since the phases
are evenly distributed for a large number of particles. Moreover, the
expressions
$\langle\tilde{x}_{1}^{2}\rangle=\langle\tilde{x}_{2}^{2}\rangle=\langle
J_{a}\rangle$ are valid for normalized coordinates such that we finally obtain

$\displaystyle\langle\Delta J_{a}\rangle$ $\displaystyle=$
$\displaystyle\frac{P_{0}ds}{cE}\left(\frac{2}{B}\frac{\partial B}{\partial
x}+\frac{1}{\rho}\right)\left[S_{11}\tilde{D}_{1}+S_{12}\tilde{D}_{2}\right]\langle
J_{a}\rangle$ $\displaystyle\langle\Delta J_{b}\rangle$ $\displaystyle=$
$\displaystyle\frac{P_{0}ds}{cE}\left(\frac{2}{B}\frac{\partial B}{\partial
x}+\frac{1}{\rho}\right)\left[S_{13}\tilde{D}_{3}+S_{14}\tilde{D}_{4}\right]\langle
J_{b}\rangle\,,$ (45)

where we recover the uncoupled case by setting $S_{11}=1$ and all other
$S_{1i}$ to zero. The expressions with the square brackets and the dispersions
$\tilde{D}$ are the projections of the coupled dispersions on the normal modes
and thus intuitively generalize the uncoupled formalism to the transversely
coupled case. The first two elements of the first row of the matrix $S$ from
eqs. 4–11 can be evaluated to be $S_{11}=g\sqrt{\beta_{a}}$ and $S_{12}=0$
such that we have

$S_{11}\tilde{D}_{1}+S_{12}\tilde{D}_{2}=g{\cal D}_{a}\ .$ (46)

with ${\cal D}_{a}$ from eq. 9. Since $S$ maps back from normalized phase
space to real space eq. 42 applies to the physical dispersion as well, and we
have

$D_{x}=\sum_{i=1}^{4}S_{1i}\tilde{D}_{i}=g{\cal
D}_{a}+S_{13}\tilde{D}_{3}+S_{14}\tilde{D}_{4}.$ (47)

${\cal D}_{a}$ is the (non-normalized) $a$-mode dispersion, the $b$-mode
dispersion is not needed. We use eq. 46 immediately to simplify the equation
for $\langle\Delta J_{a}\rangle$: we insert $P_{\gamma}$ from eq. 24 and the
definitions of curvature and focusing strength, $h=1/\rho=B/(B\rho)$ and
$k=(\partial B/\partial x)/(B\rho)$, where $k>0$ identifies a horizontally
focusing magnet.

$\frac{\langle\Delta J_{a}\rangle}{\langle
J_{a}\rangle}=\frac{C_{\gamma}E^{4}}{2\pi}\frac{1}{E}g{\cal
D}_{a}(2kh+h^{3})\,ds\,.$ (48)

We proceed to express the constant $C_{\gamma}E^{4}/2\pi$ by the total emitted
power $U_{0}$ and the second radiation integral, defined as $I_{2}=\int
h^{2}ds$

$U_{0}=\frac{C_{\gamma}E^{4}}{2\pi}\oint
h^{2}\,ds=\frac{C_{\gamma}E^{4}}{2\pi}I_{2}$ (49)

and arrive at

$\displaystyle\frac{\langle\Delta J_{a}\rangle}{\langle
J_{a}\rangle}=\frac{U_{0}}{E}\frac{1}{I_{2}}g{\cal D}_{a}(2kh+h^{3})\,ds$ (50)

for the first normal mode. Considering eq. 47, the relative damping of the
second normal mode amplitude is conveniently expressed by the damping of the
first mode amplitude and the damping in the uncoupled case.

$\frac{\langle\Delta J_{b}\rangle}{\langle J_{b}\rangle}=\frac{\langle\Delta
J_{x}\rangle}{\langle J_{x}\rangle}-\frac{\Delta J_{a}}{J_{a}}\quad\mbox{\rm
with}\quad\frac{\langle\Delta J_{x}\rangle}{\langle
J_{x}\rangle}=\frac{U_{0}}{E}\frac{1}{I_{2}}D_{x}(2kh+h^{3})\,ds$ (51)

where $D_{x}$ is the first component of the real-space dispersion from eq. 8.

So far we have considered sector dipoles such that damping of betatron
amplitudes is given by integrating eqs. 50–51 over the length of the magnet.
If, on the other hand, a magnet edge is rotated by an angle $\theta$ relative
to the edge of a pure sector magnet, with $\theta>0$ approaching a rectangular
bend (i.e. a bend where both edge angles are half the bend angle:
$\theta_{1}=\theta_{2}=\Phi/2$), the path of a particle at position $x$ is
shortened by $\Delta L=-x\tan\theta$, see fig. 1, thus its energy loss is
reduced compared to the energy loss of a particle on the design orbit by

$\Delta\bar{u}=-\frac{P_{0}}{cE}\left(1+2k\rho x\right)\Delta
L\approx\frac{P_{0}}{cE}x\tan\theta\ .$ (52)

The corresponding change in the action variable is (cf. eq. 43)

$\Delta
J_{a}=-\frac{P_{\gamma}}{cE}(\tilde{x}_{1}\tilde{D}_{1}+\tilde{x}_{2}\tilde{D}_{2})\sum_{i=1}^{4}S_{1i}\tilde{x}_{i}\tan\theta.$
(53)

Repeating the manipulations from eq. 43 to 50 gives

$\frac{\langle\Delta J_{a}\rangle}{\langle
J_{a}\rangle}=-\frac{U_{0}}{E}\frac{1}{I_{2}}g{\cal D}_{a}h^{2}\tan\theta\ .$
(54)

The relative change of the action variable in a combined function magnet is
obtained by integrating eq. 50 along the design orbit and adding the edge-
effects from eq. 54. After summing over all bending magnets in the lattice we
write

$\frac{\langle\Delta J_{a,b}\rangle}{\langle
J_{a,b}\rangle}=\frac{U_{0}}{E}\frac{I_{4{a,b}}}{I_{2}}\,,$ (55)

where we have introduced two variants of the fourth radiation integral
$I_{4a}$ and $I_{4b}$, one for each normal mode, which are given by an
integral over all bending sectors and a sum over all edges:

$\displaystyle I_{4a}$ $\displaystyle=$ $\displaystyle\oint g(s){\cal
D}_{a}(s)(2h(s)k(s)+h(s)^{3})\,ds-\sum_{k}^{\rm edges}g_{k}{\cal
D}_{ak}h_{k}^{2}\tan\theta_{k}$ (56) $\displaystyle I_{4b}$ $\displaystyle=$
$\displaystyle I_{4}-I_{4a}\ .$

$I_{4}$ is the well-known integral without coupling, and given by $I_{4a}$ for
$g=1$:

$I_{4}=\oint D_{x}(s)(2k(s)h(s)+h(s)^{3})\,ds-\sum_{k}^{\rm
edges}D_{xk}h_{k}^{2}\tan\theta_{k}$

which reproduces the well-known case (Section 3.1.4 in Zimmermann et al.
(2013)). For a small rectangular magnet ($\theta_{1}=\theta_{2}=hL/2$) with
constant field and length $L$, where the optical functions do not change much
over the magnet, eq. 56 simplifies to $\Delta I_{4a}\approx 2g{\cal
D}_{a}hkL\,.$

The second, and dominant, contribution to transverse damping comes from the
restoration of the longitudinal momentum in the radio-frequency cavities.
There, on average, the momentum corresponding to the total energy lost in one
turn $\Delta p=-U_{0}/c$ is added to the longitudinal momentum of the
electron, whereas the transverse momenta remain unchanged. Since the
transverse angles are the ratio of transverse to the much larger longitudinal
momentum, we find Lee (2004) that both transverse angles are reduced by
$-U_{0}/E.$ This effect affects both transverse coordinates equally and
therefore also the coordinates of normalized phase space $\tilde{x}$, and we
have

$\Delta\tilde{x}_{2}=-\tilde{x}_{2}\frac{U_{0}}{E}\qquad\mathrm{and}\qquad\Delta\tilde{x}_{4}=-\tilde{x}_{4}\frac{U_{0}}{E}\,,$
(57)

which causes a change in the action variables $\Delta J_{a}$ and $\Delta
J_{b}$ that, upon averaging over betatron phases, leads to

$\left.\frac{\langle\Delta J_{a}\rangle}{\langle
J_{a}\rangle}\right|_{RF}=\left.\frac{\langle\Delta J_{b}\rangle}{\langle
J_{b}\rangle}\right|_{RF}=-\frac{U_{0}}{E}$ (58)

for the relative variation of the action variable due to the acceleration in
the cavities.

Combining the effect of acceleration and emission of synchrotron radiation
from eqs. 50, 51 leads to the following expressions

$\displaystyle\left.\frac{\langle\Delta J_{a}\rangle}{\langle
J_{a}\rangle}\right|_{tot}$ $\displaystyle=$
$\displaystyle-\frac{U_{0}}{E}\left(1-\frac{I_{4a}}{I_{2}}\right)$
$\displaystyle\left.\frac{\langle\Delta J_{b}\rangle}{\langle
J_{b}\rangle}\right|_{tot}$ $\displaystyle=$
$\displaystyle-\frac{U_{0}}{E}\left(1-\frac{I_{4b}}{I_{2}}\right)\,,$ (59)

where we have to keep in mind that this is the variation of the action during
one turn in the storage ring which has the duration $T_{0}.$ Thus we find for
the damping times $\tau_{a}$ and $\tau_{b}$ for the two normal modes

$\displaystyle\frac{d\langle J_{a}\rangle}{dt}$ $\displaystyle\approx$
$\displaystyle\frac{\langle\Delta
J_{a}\rangle}{T_{0}}=-\frac{U_{0}}{T_{0}E}\left(1-\frac{I_{4a}}{I_{2}}\right)\langle
J_{a}\rangle=-\frac{2}{\tau_{a}}\langle J_{a}\rangle$
$\displaystyle\frac{d\langle J_{b}\rangle}{dt}$ $\displaystyle\approx$
$\displaystyle\frac{\langle\Delta
J_{b}\rangle}{T_{0}}=-\frac{U_{0}}{T_{0}E}\left(1-\frac{I_{4b}}{I_{2}}\right)\langle
J_{b}\rangle=-\frac{2}{\tau_{b}}\langle J_{b}\rangle$ (60)

and we recover the well-known expression Sands (1969) for the damping time,
but generalized to describe the damping of the normal modes. Note that the
conventional damping time refers to the damping of the amplitude of betatron
oscillations, whereas in eq. IV we calculate the reduction of the action
variable, which accounts for the factor two in the definition of the damping
time. The terms in brackets are commonly called the damping partition numbers
$\cal J$, which tell how the total damping is distributed among the two
transverse and the longitudinal mode.

${\cal J}_{a,b}=1-\frac{I_{4a,b}}{I_{2}}.$ (61)

As proven under very general conditions, which include coupling Robinson
(1958), the total damping is constant, and as a consequence the damping
partition numbers always sum up to 4. For the longitudinal damping partition
and the three damping times we thus get:

${\cal J}_{e}=4-{\cal J}_{a}-{\cal
J}_{b}\qquad\mathrm{and}\qquad\tau_{a,b,e}=\frac{2ET_{0}}{U_{0}{\cal
J}_{a,b,e}}$ (62)

Moreover our finding $I_{4a}+I_{4b}=I_{4}$ implies that ${\cal J}_{e}$ does
not change with transverse coupling.

##  V Emittance

In sections III and IV we derived the generalized quantum excitation and
damping integrals $I_{5}$ and $I_{4}$ including coupling. The other three
integrals $I_{1}$, $I_{2}$, $I_{3}$ do not depend on coupling. The integral
$I_{1}=\oint h(s)D_{x}(s)\,ds$ determines the energy dependent path length
(momentum compaction). Here only the dispersion in the plane of deflection is
required; the dispersion in the orthogonal plane does not affect the path
length in first order. A coordinate rotation may perform the transformation to
entry or exit of a bend oriented differently. The integrals $I_{2}=\oint
h(s)^{2}\,ds$ and $I_{3}=\oint|h(s)|^{3}\,ds$ do not depend on any optical
functions and are therefore not affected by coupling.

The equilibrium normal-mode emittances, which we identify as the ensemble
average of the action variables $\varepsilon_{a,b}=\langle J_{a,b}\rangle$,
are given by the balance of heating, as described by eq. 25 and damping in eq.
IV with the result

$\varepsilon_{a}=C_{q}\gamma^{2}\frac{I_{5a}}{I_{2}-I_{4a}}\quad\mathrm{and}\quad\varepsilon_{b}=C_{q}\gamma^{2}\frac{I_{5b}}{I_{2}-I_{4b}}\
.$ (63)

These are the same expressions found for uncoupled rings with the exception
that here we use the generalized radiation integrals that describe the
projections of excitation and damping on to the normal modes.

Turning back to the general coupled case and given the normal-mode emittances
and the decomposition of the transfer matrix it is straightforward to
determine the beam matrix in real space from that in normalized phase space.
With the normal-mode emittances $\varepsilon_{a}$ and $\varepsilon_{b}$ the
beam sigma matrix in normalized phase space $\tilde{\sigma}$ is given by

$\tilde{\sigma}=\left(\begin{array}[]{cccc}\varepsilon_{a}&0&0&0\\\
0&\varepsilon_{a}&0&0\\\ 0&0&\varepsilon_{b}&0\\\
0&0&0&\varepsilon_{b}\end{array}\right)$ (64)

and the transfer matrix $S=T^{-1}{\cal A}^{-1}$ transports the particles from
normalized phase space back to real space such that the sigma matrix $\sigma$
becomes

$\sigma=S\tilde{\sigma}S^{T}$ (65)

and the projected emittances can be extracted from $\sigma$ by taking the
determinant of the $2\times 2$ blocks on the diagonal which are the square of
the projected emittances $\hat{\varepsilon}_{x/y}.$

A parameter quantifying the coupling between the planes can be introduced
either by taking the ratio of the the normal-mode emittances
$\tilde{\kappa}=\varepsilon_{b}/\varepsilon_{a}$ or the projected emittances
$\kappa=\hat{\varepsilon}_{y}/\hat{\varepsilon}_{x}.$ Which to take in a
particular case is a matter of taste or convenience.

Since $I_{3}$ is not affected by coupling, the common formula for the energy
spread is reproduced using the longitudinal damping partition from eq. 62:

$\sigma_{e}^{2}=C_{q}\gamma^{2}\frac{I_{3}}{I_{2}{\cal
J}_{e}}=C_{q}\gamma^{2}\frac{I_{3}}{2I_{2}+I_{4a}+I_{4b}}=C_{q}\gamma^{2}\frac{I_{3}}{2I_{2}+I_{4}}$
(66)

Finally the r.m.s. bunch length follows with $I_{1}$ from the energy spread
and the synchrotron tune $Q_{s}$ (for $Q_{s}<<1$):

$\sigma_{s}=\frac{I_{1}}{2\pi Q_{s}}\sigma_{e}.$ (67)

Eqs. 65–67 give the 6-dimensional equilibrium beam parameters at any location
in a coupled lattice. They are determined by the generalized radiation
integrals together with transverse (magnet gradients) and longitudinal (RF
voltage) focusing.

##  VI Orbit distortion

Magnet misalignments, active steerers or a beam energy offset cause a
deviation of the beam centroid (orbit) from the reference path, which is
defined by the multipole centers. As a consequence the beam experiences
transverse fields and gradients in all multipoles, i.e. all magnets act as
small bending magnets and thus contribute to the radiation integrals via feed-
down of fields and gradients at the orbit position. We do not explicitely
consider skew multipoles instead they are realized by rotations of regular
multipoles.

Defining $X(s)$, $Y(s)$ as the local orbit relative to the reference
trajectory, the coordinates of a particle oscillating around the orbit are
$(X+x)$ and $(Y+y)$. Keeping only linear, regular multipoles and using again
curvatures (inverse bending radii) we define local curvatures

$h_{X}=\left.\frac{B_{y}}{(B\rho)}\right|_{x,y=0}=h_{0}+kX\qquad
h_{Y}=-\left.\frac{B_{x}}{(B\rho)}\right|_{x,y=0}=-kY,$ (68)

with $h_{0}$ the curvature of the reference trajectory and $k$ the gradient.

Generalizing the preceding treatments to an arbitrary orbit, we first realize
that the radiation integrals $I_{2}$, $I_{3}$, $I_{5}$ contain only the
quantum emission term which is given by the absolute curvature. At the orbit
it is given by

$|h|=\sqrt{h_{X}^{2}+h_{Y}^{2}}=\sqrt{h_{0}^{2}+2h_{0}kX+k^{2}(X^{2}+Y^{2})}.$
(69)

For the damping integral $I_{4}$ the variation of path length with offsets
$x,y$ to the displaced orbit at $X,Y$ has to be considered. Using local orbit
curvatures $h_{X}$, $h_{Y}$ the path length is approximated by

$dl\approx ds(1+h_{X}x+h_{Y}y),$ (70)

if we keep only linear terms in coordinates. We use $dl$ as a generic symbol
for a path length element and not as global variable.

As derived in sec. IV, $I_{4}$ depends on the radiated power as a function of
transverse position, which is proportional to the local field squared, see eq.
24: the radiation on the orbit depends on the field at the orbit and the
energy of the orbit. So the average relative energy of the emitted photons
from eq. 40 is given by

$\Delta\bar{u}=\frac{P_{\gamma}}{E}\frac{dl}{c}=\hat{C}\left[(h_{X}+kx)^{2}+(h_{Y}-ky)^{2}\right]\cdot(1+h_{X}x+h_{Y}y)\,ds$
(71)

with $\hat{C}=C_{\gamma}E^{3}/(2\pi)$ and the path length element from eq. 70,
now also including the relative offset to the orbit. The first bracket
corresponds to the squared magnetic field at the location of the particle.
Keeping only linear terms in the coordinates $x,y$ we get

$\displaystyle\Delta\bar{u}$ $\displaystyle=$
$\displaystyle\hat{C}\left(C_{0}+C_{x}x+C_{y}y\right)\,ds,\quad\mbox{with}$
(72) $\displaystyle C_{0}$ $\displaystyle=$ $\displaystyle|h|^{2}$ (73)
$\displaystyle C_{x}$ $\displaystyle=$ $\displaystyle h_{X}\,(|h|^{2}+2k)$
$\displaystyle C_{y}$ $\displaystyle=$ $\displaystyle h_{Y}\,(|h|^{2}-2k)$

where we used eq. 68 and eq. 69. We insert 72 into eq. 39 describing the
change of betatron amplitudes due to radiation emission. Further we insert eq.
42 expressing $x$ in normalized coordinates $\tilde{x}_{i}$, and the
corresponding expression for $y$:

$\Delta
J_{a}=\left(\tilde{x}_{1}\tilde{D}_{1}+\tilde{x}_{2}\tilde{D}_{2}\right)\hat{C}\left[C_{0}+C_{x}\sum_{i=1}^{4}S_{1i}\tilde{x}_{i}+C_{y}\sum_{i=1}^{4}S_{3i}\tilde{x}_{i}\right]ds\,.$
(74)

Averaging over betatron phases (cf. treatments following eq. 43) and
performing the corresponding procedure for $J_{b}$ results in

$\displaystyle\langle\Delta J_{a}\rangle$ $\displaystyle=$
$\displaystyle\hat{C}\left[C_{x}\left(S_{11}\tilde{D}_{1}+S_{12}\tilde{D}_{2}\right)+C_{y}\left(S_{31}\tilde{D}_{1}+S_{32}\tilde{D}_{2}\right)\right]ds\,\langle
J_{a}\rangle$ $\displaystyle\langle\Delta J_{b}\rangle$ $\displaystyle=$
$\displaystyle\hat{C}\left[C_{x}\left(S_{13}\tilde{D}_{3}+S_{14}\tilde{D}_{4}\right)+C_{y}\left(S_{33}\tilde{D}_{3}+S_{34}\tilde{D}_{4}\right)\right]ds\,\langle
J_{b}\rangle$ (75)

With the relations between normalized, decoupled and real space dispersions
(see eqs. 46, 47 and the corresponding ones for $D_{y}$, ${\cal D}_{b}$) we
get

$\displaystyle\frac{\langle\Delta J_{a}\rangle}{\langle J_{a}\rangle}$
$\displaystyle=$ $\displaystyle\hat{C}\left[C_{x}g{\cal
D}_{a}+C_{y}\left(D_{y}-g{\cal D}_{b}\right)\right]ds$
$\displaystyle\frac{\langle\Delta J_{b}\rangle}{\langle J_{b}\rangle}$
$\displaystyle=$ $\displaystyle\hat{C}\left[C_{y}g{\cal
D}_{b}+C_{x}\left(D_{x}-g{\cal D}_{a}\right)\right]ds$ (76)

On the reference orbit ($X=Y=0$) we get $C_{x}=h_{0}(2k+h_{0}^{2})$, $C_{y}=0$
and recover eq. 48. For a pure quadrupole we simply set $h_{0}=0$.

If the magnet has a rotated edge of angle $\theta$, then the path is shortened
by $\Delta L=-x\tan\theta$, see Fig. 1 – now the orbit at $X$ is the
reference. Using $\Delta L$ instead of $dl=ds(1+h_{X}x+h_{Y}y)$ we repeat the
calculations and arrive at an expressions that corresponds to eq. 73 for the
edge angles:

$\bar{C}_{0}=0\qquad\bar{C}_{x}=-\tan\theta\,|h|^{2}\qquad\bar{C}_{y}=0$ (77)

The expressions are the same like on-axis, cf. 54.

The radiation integrals also include $I_{1}$, which actually is not related to
radiation but describes the energy dependent path length. Including an energy
offset $u$ the path length element from eq. 70 becomes

$dl=ds(1+(h_{X}D_{x}+h_{Y}D_{y})u)$

Any higher orders of path length (i.e. non-linear momentum compaction) are
contained in the local dispersion $D_{x}$.

Collecting our results we obtain Table 1 giving a summary of the contributions
to the radiation integrals from a general (combined function) bending magnet
with (constant) curvature $h_{0}$, quadrupole strength $k$ and edge angles
$\theta_{1,2}$ including coupling and orbit distortions. Our results agree
with prior work by Sagan as documented in Sagan .

Table 1: Radiation integrals for a general bending magnet with curvature
$h_{0}$ and focusing strength $k$ including coupling and orbit distortions.

$\displaystyle I_{1}$ $\displaystyle=$
$\displaystyle\int(h_{X}D_{x}+h_{Y}D_{y})\,ds\qquad
I_{2}=\int|h(s)|^{2}\,ds\qquad I_{3}=\int|h(s)|^{3}\,ds$ $\displaystyle
I_{4a}$ $\displaystyle=$ $\displaystyle\int\left[C_{x}(s)g{\cal
D}_{a}(s)+C_{y}(s)\left(D_{y}(s)-g{\cal
D}_{b}(s)\right)\right]\,ds-\sum_{i=1,2}\tan\theta_{i}\bar{C}_{xi}g{\cal
D}_{ai}$ $\displaystyle I_{4b}$ $\displaystyle=$
$\displaystyle\int\left[C_{y}(s)g{\cal D}_{b}(s)+C_{x}(s)\left(D_{x}(s)-g{\cal
D}_{a}(s)\right)\right]\,ds-\sum_{i=1,2}\tan\theta_{i}\bar{C}_{xi}\left(D_{xi}-g{\cal
D}_{ai}\right)$ $\displaystyle I_{5a}$ $\displaystyle=$
$\displaystyle\int\vec{\tilde{D}}_{a}^{T}(s)\cdot\vec{\tilde{D}}_{a}(s)\,|h(s)|^{3}\,ds\qquad
I_{5b}=\int\vec{\tilde{D}}_{b}^{T}(s)\cdot\vec{\tilde{D}}_{b}(s)\,|h(s)|^{3}\,ds$

${\rm with}\qquad h_{X}(s)=h_{0}+kX(s)\qquad
h_{Y}(s)=-kY(s)\qquad|h(s)|=\sqrt{(h_{X}(s))^{2}+(h_{Y}(s))^{2}}$

$C_{x}(s)=h_{X}(s)\,(|h(s)|^{2}+2k(s))\qquad
C_{y}(s)=h_{Y}(s)\,(|h(s)|^{2}-2k(s))\qquad\bar{C}_{x}=|h_{i}|^{2}$

In the calculations we made the following assumptions and simplifications:

(a) The integrals extend along the displaced orbit. In a magnet of cylindrical
symmetry with design curvature $h_{0}$, the integration path element thus
becomes

$dl=ds\,(1+h_{0}X(s)).$

If the magnet further has a rotated edge of angle $\theta$, then the path is
shortened by

$\Delta L=-X\tan\theta,$

see Fig. 1. An integral over a function $f(s)$ thus becomes

$\int_{\rm path}f(s)\,dl=\int_{\Delta L_{1}}^{L-\Delta
L_{2}}f(s)\,(1+h_{0}X)\,ds\approx\int_{0}^{L}f(s)\,(1+h_{0}X)\,ds-\sum_{i=1,2}f(s_{i})\tan\theta_{i}\,X_{i}.$

Basically this path length correction has to be applied to all radiation
integrals. However, for storage rings operating at several GeV beam energy,
the radius of curvature $\rho=1/h_{0}$ is on the orders of meters, while the
orbit excursion $X$ is typically on the order of some millimeter. As a
consequence, the path length modifications are a per-mill effect and may be
neglected. For rectangular bends, as commonly used, where
$\theta_{1}=\theta_{2}=h_{0}L/2\ll 1$ the two effects, lengthening of the path
in the arc and shortening due to cut off by rotated edges, cancel anyway.

(b) It would be straightforward to extend the treatment to non-linear
multipoles too, however the contribution to the integrals is rather small: the
field increases quadratically (sextupole) or at even higher power with the
distance from the axis, and the radiation integrals scale with the second or
third power of the field. Even if the sextupoles are rather strong, like in a
low-emittance storage ring, their contribution to the integrals is on a per-
mill level compared to the quadrupole contribution – which in itself is a
fraction of the bending magnet contribution – and thus may well be neglected.

(c) We only consider the design edge rotation (around the $y$-axis) of a
bending magnet and neglect additional edge effects due to orbit slopes
$X^{\prime}$, $Y^{\prime}$, because these angles are usually small compared to
bending and edge angles. Thus the edge related path shortening does not depend
on the vertical coordinates. Note that we only considered horizontal bending
magnets, which may be sandwiched between rotations to model vertical or
slanted deflections.

(d) During propagation the effect known as “mode flip” may occur: the normal
mode decomposition locally may result in one or two solutions for the matrix
$T$ of Eq. 5. In strongly coupled lattices it can happen, that the solution in
use does not exist anymore at the next element, and the other solution has to
be used to propagate further, as described in Sagan and Rubin (1999). In this
case we proceed by adding the $a$-integrals to the $b$-sum integral of the
ring and vice versa, until the originally used solution exists again. However,
since we restricted ourselves to regular magnets all changes of coupling (as
expressed by the parameter $g$) occur at rotations, which are inserted before
and after the magnets if needed, and no mode flip can occur during integration
inside the magnet.

##  VII Applications

Several simulation codes are in use, e.g. MAD-X MADX , elegant ELEGANT , AT AT
, BMAD BMAD , Tracy TRACY , which incorporate different coupling formalisms
Chao (1979a)–Sagan and Rubin (1999) and calculate the equilibrium parameters
for a given storage ring lattice. The code OPA OPA  emerged from the former
code OPTIK, which was written by K. Wille and colleagues at Dortmund
university in the 1980’s. The philosophy of the code was different providing
an interactive, visual tool to build up a lattice from scratch like playing
with LEGOTM blocks, while restricting the models to elementary magnet types
and functionalities.

For interactive lattice design from scratch it is convenient, although
physically meaningless, to calculate radiation integrals even if it is not yet
known if and how the lattice will be terminated or closed, i.e. irrespective
of the existence of a periodic solution. For this purpose a forward
propagation using analytical expressions for the integrals is required, which
starts from initial beam parameters defining orbit, normal mode betas,
coupling and dispersion. If a periodic solution exists, it will deliver the
initial parameters and the same calculation can be applied. The OPA-code is
based on this concept. Extending it to coupled lattices partially motivated
the work described here.

The presented method to calculate the radiation integrals for the normal modes
permits a straightforward implementation of the formulae in Table 1: We assume
that curvature $h$ and gradient $k$ are piecewise constant, then the transfer
matrix $R$ and the vector $\vec{V}$ of dispersion production in eq. 1 contain
trigonometric or hyperbolic functions (to be found in any text book on
accelerators, e.g.Edwards and Sypers (1993), Zimmermann et al. (2013), Lee
(2004)), whereas otherwise analytic solutions can only be found for a few
special cases of longitudinal field variation.

Constant $h$, $k$ parameters can be extracted from the radiation integrals,
and we are left with sums of definite integrals over multiple products of
trigonometric or hyperbolic functions, which have elementary solutions. These
solutions become rather lengthy and are not instructive at all, therefore we
do not show them here. But they can be evaluated and exported as program code
by a symbolic code, which was the first and most efficient implementation.
Solving the integrals numerically by dividing the magnets into sufficiently
small slices instead of evaluating the analytical solutions is less elegant
and computationally more expensive, but easy and straightforward to implement
and to extend further, in order to also include orbit distortions,
longitudinal field variations etc.

Linear beam optics is calculated from the local transfer-matrix, which is the
Jacobian of the non-linear map at the orbit position. The orbit itself is the
fixpoint of the one-turn map in a periodic system, or just the propagation of
initial conditions including all multipoles in a single-pass system. Off-axis
down-feeds have to be included in linear optics, like the additional
dispersion production on the displaced orbit,
$\Delta\vec{D}(s)=d\vec{X}(s)/du|_{u=0}$, and local quadrupole and skew
quadrupole gradients from non-linear multipoles.

Implementation details and practical issues will be described elsewhere OPA .
In the following we will show four example applications:

Figure 2: Schematic layout of the Spiral-COSAMI lattice (from Wrulich et al.
(2021)).

Example 1: The Spiral-COSAMI storage ring is a machine for industrial
application of extreme ultraviolet light Wrulich et al. (2021). It is build
from three vertically stacked turns as sketched in Fig. 2. The tilting of the
arcs introduces coupling and excites a wave of vertical dispersion as shown in
Fig. 3. The figure shows the normal mode beta functions which almost coincide
with the projections to physical $x,y$ space since coupling is low. The
physical dispersions also shown in the figure are defined in the local
coordinate system and thus show discontinuities at locations of coordinate
rotation. The normal mode equilibrium emittances of this lattice amount to
3.43/0.23 nm at 430 MeV. This result was confirmed by the code TRACY TRACY ,
which employs a different method of emittance calculation based on the
6-dimensional periodic sigma-matrix Chao (1979b).

Figure 3: Optics of the Spiral COSAMI lattice showing normal mode beta
functions (left axis) and physical dispersions (right axis). The triangles at
the bottom indicate points of rotations by angles
$\\{-\beta,\alpha,\alpha,\alpha,\alpha,-\beta\\}$ with $\alpha=3.86^{\circ}$,
$\beta=7.69^{\circ}$

Example 2: The lattice for the upgrade of the Swiss Light Source, SLS 2.0
Braun et al.  has a natural horizontal emittance of 158 pm at 2.7 GeV. 264
skew quadrupoles are set to generate closed bumps of vertical dispersion in
the arcs in order to create 10 pm of $b$-mode emittance, which is a compromise
between beam life time due to Touschek scattering and photon beam brightness.
Due to the rather low coupling the vertical emittance basically is given by
the $b$-mode emittance. Closing all undulators to minimum gaps increases the
radiated power as given by the $I_{2}$-integral while little affecting the
other integrals, and thus reduces the normal mode emittances to 134 pm and 7.7
pm. Scaling all skew quadrupoles by the same factor proportional to
$\sqrt{I_{2}}$ keeps the $b$-mode emittance at 10 pm for the fully loaded
lattice. Fig. 4 shows the normal mode beta functions $\beta_{a}$, $\beta_{b}$
and the physical dispersions $D_{x}$, $D_{y}$ for one super-period of the
fully loaded, ideal lattice. Fig. 5 demonstrates for the unloaded lattice how
the vertical emittance increases due to orbit distortion when exciting one
vertical corrector. We tracked the residual discrepancy between the results
from TRACY and OPA to slightly different methods to calculate the off-axis
dispersion.

Figure 4:  Normal normal mode beta functions (blue $a$-mode, red $b$-mode) and
dispersions (green horizontal, magenta vertical) for setors 6 to 9 of the SLS
2.0 lattice loaded with insertion devices. Figure 5:  Comparison of OPA- and
TRACY calculations of vertical emittance due to orbit distortion, if one
vertical corrector is excited (located at $\beta_{y}=7.6~{}m$).

Example 3: As an example of an extremely coupled machine we designed a ring in
which we can continuously adjust the coupling between an uncoupled
configuration and a Möbius configuration Talman (1995), where the transverse
planes are exchanged after one turn. This racetrack-shaped ring Ziemann  is
based on 10 m long FODO cells with $90\mbox{${}^{\circ}$}$ phase advance in
both planes. In the two arcs the phase advance in the vertical plane is
slightly reduced in order to split the integer tunes, while the phase advance
in the horizontal plane remains at $90\mbox{${}^{\circ}$}$, which allows us to
implement a dispersion suppressor with four half-length dipoles. The two
straight sections consist of six $90\mbox{${}^{\circ}$}$ FODO cells. One is
used to adjust tunes and the other implements a Möbius tuner with three skew
quadrupoles placed in the middle of corresponding drifts spaces in consecutive
cells Aiba et al. (2015). This ensures that the phase advance between
consecutive skew quadrupoles is $90\mbox{${}^{\circ}$}$ in both planes. If we
now choose the focal length of the skew quadrupole to be equal to the beta
function in the middle of the straight, the transfer matrix for this section
has $2\times 2$ blocks with zeros on the diagonal and non-zero entries in the
off-diagonal blocks. In other words, it exchanges the transverse planes.
Adjusting the excitations of the three skew quadrupoles by the same factor
allows us to continuously vary the coupling between an uncoupled ring and one
in the Möbius configuration.

Figure 6: The fractional tunes (top) and the emittances (middle) of the Möbius
ring (top) as a function of the excitation of the skew quadrupoles. The bottom
panel shows the sum and difference of the tunes.

The upper panel in Figure 6 shows the fractional tunes for the ring as the
excitation of the skew quadrupoles is increased from zero to the Möbius
configuration and up to 1.5 times that excitation. The integer part of the
tune is indicated on the left-hand side of the plot. We see that increasing
the excitation of the skew quadrupoles “pushes the fractional tunes apart”,
which makes it beneficial to place one tune above the half-integer and the
other below. This prevents the crossing of the half integer with the ensuing
instability for one of the tunes. The middle panel shows the emittances
calculated with Eq. 63. We observe that increasing the excitation of the skew
quadrupoles to the Möbius configuration reduces one of the emittances and
increases the other one such that they become equal. The lower panel shows
that the difference between the fractional part of the tunes becomes half-
integer at this point. Increasing their excitation further causes both
emittances to increase dramatically. The plot on the lower panel provides an
explanation; the sum of the tunes $Q_{a}+Q_{b}$ becomes an integer and the
system crosses a sum resonance shown by the asterisks in the lower panel from
Figure 6. And on a sum resonance, the emittances can become arbitrarily large,
because only their difference is bounded (Sec.2.1.3 in Zimmermann et al.
(2013)).

Example 4: A Möbius insertion in a light source is an option to lower the
bunch density in order to minimize emittance growth and particle losses
(Touschek effect) due to intra-beam scattering. Furthermore, some beam lines
prefer photon beams of almost round cross section obtained in this way. Our
last example demonstrates application to the SLS 2.0 lattice: The Möbius
insertion could be realized in one of the of the short straight sections. It
is made from five skew and two regular quadrupoles Aiba et al. (2015). Figure
7 shows the optical functions in magnification: normal mode beta functions
jump at the location of the (thin) skew quadrupoles, whereas their projections
to $x$\- and $y$-axes, which correspond to the physical beam envelopes, vary
continuously. Outside the Möbius insertion the projections from the $a$\- and
$b$-modes to $x$, resp. to $y$ coincide, indicating full coupling at
$g=\sqrt{1/2}$. As a consequence the radiation integrals for both modes are
identical, and so are the damping times and the emittances. The uncoupled
lattice without insertion devices has a natural emittance of
$\epsilon_{xo}=158$ pm and damping partitions of ${\cal J}_{xo}=1.828$, ${\cal
J}_{yo}=1.000$. The Möbius lattice has ${\cal J}_{a}={\cal J}_{b}=({\cal
J}_{xo}+{\cal J}_{yo})/2=1.414$ and emittances of
$\epsilon_{a}=\epsilon_{b}=\epsilon_{xo}/(1+{\cal J}_{yo}/{\cal J}_{xo})=102$
pm. These result in projected emittances of same amount,
$\epsilon_{x}=\epsilon_{y}=102$ pm outside the insertion, whereas projected
emittances are larger inside the insertion.

Figure 7:  Möbius insertion in straight 4 of the SLS 2.0 storage ring. Solid
lines show the normal mode betafunctions (blue $a$-mode, red $b$-mode), dashed
and dotted lines the projected beta functions $\beta_{xa},\beta_{yb}$ and
$\beta_{xb},\beta_{ya}$ (blue $\beta_{x}$, red $\beta_{y}$). The square line
shows the determinant of the coupling matrix $\det C=1-g^{2}$ (axis at right).

##  VIII Conclusion

We found a generalization of the synchrotron radiation integrals that are
responsible for the emittance growth in coupled lattices. It turned out that
the resulting methodology is rather intuitive. It consist of projecting the
four dimensional coupled dispersion from eq. 8 onto the normal modes that are
constructed from the Edwards-Teng formalism. Then the emittance growth
integral is given through the normal-mode beta functions and dispersion
$\vec{\tilde{D}}$. The same method also worked to project the effect of
damping due to the emission of synchrotron radiation on to the normal modes.
The formalism allows off-axis contribution to be included, and implementation
in a beam dynamics code is straightforward. In this way the formalism allows
to introduce the coupling of the normal-mode emittances in a natural way.

## Acknowledgements

We would like to thank Masamitsu Aiba (PSI) for setting up the Möbius
insertion in the latest SLS 2.0 lattice, and for useful discussions, and
Michael Böge (PSI) and Bernard Riemann (PSI) to perform cross-checks of our
example applications using the TRACY code. We further gratefully acknowledge
helpful comments from David C. Sagan (Cornell University).

## References

  * Sands (1969) M. Sands, Conf. Proc. C6906161, 257 (1969).
  * Helm et al. (1973) R. H. Helm, M. J. Lee, P. L. Morton, and M. Sands, IEEE Trans. Nucl. Sci. 20, 900 (1973).
  * Edwards and Sypers (1993) D. Edwards and M. Sypers, _An Introduction to the Physics of High Energy Accelerators_ (John Wiley & Sons, Inc., New York, 1993).
  * Chao (1979a) A. W. Chao, J. Appl. Phys. 50, 595 (1979a).
  * Ohmi et al. (1994) K. Ohmi, K. Hirata, and K. Oide, Phys. Rev. E49, 751 (1994).
  * Wolski (2006) A. Wolski, Phys. Rev. ST Accel. Beams 9, 024001 (2006).
  * (7) OPA, _Project web site:http://ados.web.psi.ch/opa/_.
  * Wrulich et al. (2021) A. Wrulich, A. Streun, and L. Rivkin, Nucl. Instrum. Meth. A 1014, 165731 (2021).
  * Edwards and Teng (1973) D. Edwards and L. Teng, IEEE Transactions on Nuclear Science 20, 885 (1973).
  * Sagan and Rubin (1999) D. Sagan and D. Rubin, Phys. Rev. ST Accel. Beams 2, 074001 (1999), [Phys. Rev. ST Accel. Beams3,059901(2000)].
  * Ziemann (2019) V. Ziemann, _Hands-On Accelerator Physics Using MATLAB ®_ (CRC Press, 2019), 1st ed.
  * Zimmermann et al. (2013) F. Zimmermann, K. Mess, and M. Tigner, _Handbook of Accelerator Physics and Engineering_ (World Scientific, 2013), 2nd ed.
  * Lee (2004) S. Y. Lee, _Accelerator Physics_ (World Scientific, 2004), 2nd ed.
  * Robinson (1958) K. W. Robinson, Phys. Rev. 111, 373 (1958), [205(1958)].
  * (15) D. Sagan, _The BMAD reference manualBMAD _.
  * (16) MADX, _Project web site:http://madx.web.cern.ch/madx/_.
  * (17) ELEGANT, _Project web site:https://ops.aps.anl.gov/elegant.html_.
  * (18) AT, _Project web site:http://atcollab.sourceforge.net/_.
  * (19) BMAD, _Project web site:https://www.classe.cornell.edu/bmad/_.
  * (20) TRACY, _Project web site:https://github.com/jbengtsson/tracy-3.5/_.
  * Chao (1979b) A. W. Chao, Journal of Applied Physics 50, 595 (1979b).
  * (22) H. Braun, T. Garvey, M. Jörg, A. Ashton, P. Wilmott, and R. Kobler, _SLS 2.0 storage ring technical design report, PSI 2021,https://www.dora.lib4ri.ch/psi/islandora/object/psi:39635_.
  * Talman (1995) R. Talman, Phys. Rev. Lett. 74, 1590 (1995).
  * (24) V. Ziemann, _Beam parameters in a Möbius ring, FREIA Report 2022/01, Jan. 2022,https://arxiv.org/abs/2201.01083_.
  * Aiba et al. (2015) M. Aiba, M. Ehrlichman, and A. Streun, in _Proc. IPAC 2015, Richmond VA, USA_ (2015), pp. 1716–1719.


# Super heavy dark matter from inflationary Schwinger production

Mar Bastero-Gil<EMAIL_ADDRESS>Departamento de Física Teórica y del Cosmos and
CAFPE, Universidad de Granada, Campus de Fuentenueva, E-18071 Granada, Spain
Paulo B. Ferraz<EMAIL_ADDRESS>Departamento de Física Teórica y
del Cosmos and CAFPE, Universidad de Granada, Campus de Fuentenueva, E-18071
Granada, Spain Universidad de Coimbra, Faculdade de Ciências e Tecnologia da
Universidade de Coimbra and CFisUC, Rua Larga, 3004-516 Coimbra, Portugal
Lorenzo Ubaldi<EMAIL_ADDRESS>Jožef Stefan Institute, Jamova 39, 1000
Ljubljana, Slovenia Institute for Fundamental Physics of the Universe (IFPU),
Via Beirut 2, 34014 Trieste, Italy Roberto Vega-Morales<EMAIL_ADDRESS>
Departamento de Física Teórica y del Cosmos and CAFPE, Universidad de Granada,
Campus de Fuentenueva, E-18071 Granada, Spain

###### Abstract

We consider a simple setup with a dark sector containing dark electrons
charged under an abelian $U(1)_{D}$ gauge symmetry. We show that if the
massless dark photon associated to the $U(1)_{D}$ is produced during inflation
in such a way as to form a classical dark electric field, then dark electron-
positron pairs are also produced close to the end of inflation via the
Schwinger mechanism. Their mass is larger than the Hubble scale so they are
non-relativistic at production and remain so for the rest of their cosmic
evolution following reheating. They can account for the dark matter abundance
today in scenarios with a high scale of inflation and efficient reheating.
This points to a very heavy mass for the dark electrons of order $10^{14}$
GeV. Even with a dark gauge coupling comparable to the visible one of
electromagnetism, the dark electrons are so heavy that they do not thermalize
with the dark photons throughout their cosmic history. Assuming negligible
kinetic mixing with the visible $U(1)$ the dark electrons remain decoupled
from the Standard Model thermal bath as well.

††preprint: UG-FT 328-23, CAFPE 198-23

##  I Introduction

The nature and production mechanism of dark matter is still unknown.
Significant research effort has been put into candidates like weakly
interacting massive particles, axions, sterile neutrinos, but we have seen no
experimental evidence for them yet Bertone and Tait (2018). This motivates
studying alternative candidates as well. One class which has received
considerable attention in the recent past is that of dark photon dark matter,
produced via non-thermal mechanisms Graham et al. (2016); Agrawal et al.
(2020); Dror et al. (2019); Co et al. (2019); Bastero-Gil et al. (2019); Long
and Wang (2019); McDermott and Witte (2020); Nakai et al. (2020); Ahmed et al.
(2020); Kolb and Long (2021); Salehian et al. (2021); Bastero-Gil et al.
(2022); Gorghetto et al. (2022); Sato et al. (2022); Redi and Tesi (2022);
Barrie (2022) where the dark photon is associated with a dark $U(1)_{D}$
abelian gauge theory. At some point before matter radiation equality the dark
photon must become massive in order to provide a cold (non relativistic) dark
matter candidate. However, it has been argued East and Huang (2022) that the
formation of vortices in these scenarios spoils the possibility of the dark
photon being dark matter.

One can consider adding scalar and/or fermion fields charged under the
$U(1)_{D}$ to the dark sector, of which the dark photon is the gauge boson.
Such a ‘dark QED’ setup was studied in Arvanitaki et al. (2021) where the
authors considered a massive dark photon whose longitudinal component is
produced via inflationary fluctuations as in Graham et al. (2016) which is
coupled to dark electrons. They showed how via various processes, such as dark
electromagnetic cascades and plasma dynamics, the dark electrons are produced
and eventually thermalize with the dark photons while remaining decoupled from
the visible thermal bath. In the end it was found that the dark electrons can
constitute the dark matter in a mass range between 50 MeV and 30 TeV. An
important point Graham et al. (2016); Arvanitaki et al. (2021) is that the
dark photons produced during inflation form a ‘condensate’ which can only be
treated as a classical dark electric field in the radiation dominated era once
their momentum modes re-enter the horizon. It is only then that the processes
leading to the production of dark electrons and to thermalization can begin
taking place.

In this letter we consider a ‘dark QED’ scenario as well, but utilize a
mechanism where the classical dark electric field is present already by the
end of inflation Bastero-Gil et al. (2019, 2022). We show that in such a case
it is possible to generate by the end of inflation a relic of dark electrons
which today matches the observed dark matter abundance. The dark photon can be
massless as opposed to the scenario in Graham et al. (2016); Arvanitaki et al.
(2021) where it is necessarily massive. Furthermore, the dark electrons are
produced entirely via the Schwinger effect as opposed to via thermalization
within the dark sector as in Graham et al. (2016); Arvanitaki et al. (2021).
Requiring the dark electrons to be non-relativistic by the time of matter
radiation equality pushes their mass to be heavier than the Hubble scale at
the end of inflation which itself must be $\sim 10^{14}\,GeV$.

##  II Generating a dark electric field

We consider an abelian $U(1)_{D}$ group whose gauge boson is the dark photon.
A simple way of generating a classical dark electric field during inflation is
by considering the coupling $\frac{\alpha}{4f}\phi
F_{\mu\nu}\tilde{F}^{\mu\nu}$ where $\phi$ is the inflaton, $F_{\mu\nu}$ is
the dark photon field strength Bastero-Gil et al. (2019, 2022). The coupling
leads to the exponential production of one of the two transverse polarizations
of the dark photon. The typical wavelength of the produced dark photons is
roughly the size of the Hubble horizon. This produces a classical dark
electric field Anber and Sorbo (2010); Bastero-Gil et al. (2022) that at the
end of inflation has magnitude approximately given by,

$E\approx 10^{-2}\frac{e^{\pi\xi_{\rm end}}}{\xi_{\rm end}^{3/2}}H_{\rm
end}^{2}\,,$ (1)

which is coherent within the horizon with energy density,

$\rho_{E}=\frac{1}{2}E^{2}.$ (2)

Here $\xi\equiv\frac{\alpha\dot{\phi}}{2Hf}$ with $H$ the Hubble parameter.
The subscript ‘end’ indicates that the quantities are evaluated at the end of
inflation. We can consider values of $\xi_{\rm end}$ roughly in the range
$1<\xi_{\rm end}<10$ in which case we easily evade constraints from the Cosmic
Microwave Background (CMB) and we can neglect the backreaction of the produced
dark electric field on the inflaton dynamics Bastero-Gil et al. (2022). Note
that in such a range we have a dark electric field with magnitude larger than
$H_{\rm end}^{2}$. A dark magnetic field is produced as well, but its
magnitude is suppressed by $\xi_{\rm end}^{2}$ compared to $E$ so we neglect
it. In this mechanism for generating the dark electric field, the dark photon
can be massless as the production involves the transverse polarizations. This
is in contrast to mechanisms explored in Graham et al. (2016); Arvanitaki et
al. (2021) where the dark photon is inevitably massive.

##  III Dark Schwinger production

We now introduce the dark electrons which we label $\chi$. We parametrize the
background de-Sitter spacetime by the Friedmann-Robertson-Walker metric with
$ds^{2}=-dt^{2}+a^{2}(t)d\vec{x}^{2}$ and consider the action,

$\displaystyle S$ $\displaystyle=$ $\displaystyle\int
d^{4}x\sqrt{-g}\left(-\frac{1}{4}F^{\mu\nu}F_{\mu\nu}+\,\mathcal{L}_{\rm{charged}}(A_{\mu},\chi)\right)\,.~{}~{}$
(3)

For simplicity and for concreteness we focus on scalar dark electrons in this
work, but our results also apply to fermions up to numerical factors of order
one. We thus consider the lagrangian,

$\displaystyle\mathcal{L}_{\rm charged}^{\rm
scalar}=-|D_{\mu}\chi|^{2}-m_{\chi}^{2}|\chi|^{2}\,,$ (4)

where we have the covariant derivative $D_{\mu}=\partial_{\mu}+ig_{D}A_{\mu}$,
with $g_{D}$ the $U(1)_{D}$ dark gauge coupling.

In the presence of a constant background dark electric field $E$ during
inflation, such as the one in Eq. (1), pairs of dark electrons and positrons
are produced via the Schwinger effect. To see this we rely on previous studies
of magnetogenesis Kobayashi and Afshordi (2014); Hayashinaka et al. (2016);
Banyeres et al. (2018) which computed the Schwinger current $J$ of the
generated $\chi$ charged particles. To estimate the dark electron energy
density $\rho_{\chi}$ we consider Sobol et al. (2019) the evolution equation,

$\displaystyle\dot{\rho}_{\chi}+3H\rho_{\chi}=E\cdot J=\sigma E^{2}\,.$ (5)

Here we are treating the produced $\chi$’s as non relativistic (this will be
justified in the next section) and we have introduced the conductivity
$\sigma=J/E$. This evolution equation can be derived from the energy momentum
tensor and expresses energy conservation where $E\cdot J$ acts as a source
term for the $\chi$ energy density. A corresponding term with the opposite
sign appears in the evolution equation for the energy density of the dark
electric field which is sourced by the inflaton, but depleted by the $-\sigma
E^{2}$ term. The latter is negligible because $\sigma/H$ is small, as we will
discuss below. To a good approximation then the dark electric field is
unaffected by the produced dark electrons.

To ensure Schwinger production dominates over purely gravitational particle
production, we consider the strong dark electric force regime where
$g_{D}E>H^{2}$. The dimensionless conductivity in this regime is given by
Kobayashi and Afshordi (2014),

$\bar{\sigma}\equiv\frac{\sigma}{H}=\frac{g_{D}^{2}}{12\pi^{3}}\frac{g_{D}E}{H^{2}}{\rm
Exp}\left[-\pi\frac{m_{\chi}^{2}}{g_{D}E}\right]\,.$ (6)

As we will see in the next section, we have $m_{\chi}^{2}>g_{D}E$ so the
conductivity is exponentially suppressed, $\bar{\sigma}\ll 1$. As we consider
the Schwinger production during inflation, it is convenient to rewrite Eq. (5)
switching to conformal time, $d\tau=dt/a$, and use $\tau=-1/(aH)$ to obtain,

$\displaystyle\partial_{\tau}\rho_{\chi}-\frac{3}{\tau}\rho_{\chi}=-\bar{\sigma}E^{2}\frac{1}{\tau}.$
(7)

Evolving from $\tau_{0}$, taking $\rho_{\chi}(\tau=\tau_{0})=0$ as the initial
condition, to the end of inflation $\tau_{\rm end}$ one finds,

$\displaystyle\rho_{\chi}^{\rm end}\equiv\rho_{\chi}(\tau_{\rm
end})=\frac{\bar{\sigma}}{3}E^{2}\Big{[}1-\Big{(}\frac{\tau_{0}}{\tau_{\rm
end}}\Big{)}^{-3}\Big{]}\simeq\frac{\bar{\sigma}_{\rm end}}{3}E^{2}\,.~{}~{}$
(8)

We have used $\tau_{0}/\tau_{\rm end}=e^{N_{e}}$, with $N_{e}$ the number of
e-folds between $\tau_{0}$ and $\tau_{\rm end}$, to reach the last approximate
equality. Even if we have a dark electric field present only for a few e-folds
close to the end of inflation, the solution in Eq. (8) holds. Since
$\bar{\sigma}\ll 1$, we also have that $\rho_{\chi}^{\rm end}$ is much smaller
than the energy density of the electric field which justifies neglecting the
backreaction of $\chi$ on $E$.

##  IV Cosmic evolution

For simplicity, we approximate the reheating process that takes us from the
end of inflation to the radiation dominated era as instantaneous. At the end
of inflation the inflaton dominates the energy density of the Universe,

$\displaystyle\rho^{\rm{end}}_{I}=3H_{\rm end}^{2}M_{\rm Pl}^{2}.$ (9)

It must transfer most of this energy density into visible radiation whose
energy density we parametrize as,

$\displaystyle\rho_{R}(T_{\rm RH})=\frac{\pi^{2}}{30}g_{*}(T_{\rm RH})T_{\rm
RH}^{4}\,\equiv\epsilon_{R}^{4}\rho^{\rm{end}}_{I}=3H_{\rm RH}^{2}M_{\rm
Pl}^{2}\,,~{}~{}$ (10)

where $\epsilon_{R}<1$ parametrizes the efficiency of reheating. Here $T_{\rm
RH}$ denotes the reheating temperature while $g_{*}(T_{\rm RH})$ is the number
of relativistic degrees of freedom which we fix to $g_{*}(T_{\rm RH})\sim 100$
and restrict ourselves to reheating temperatures above the electroweak scale.
Given the inflaton energy density at the end of inflation in Eq. (9) we see
Eq. (10) also defines the relation $H_{\rm RH}=\epsilon_{R}^{2}H_{\rm end}$
which accounts for the fact that the Hubble scale at the end of reheating will
be lower than at the end of inflation. At reheating we take the energy
densities of dark electric field and dark electrons as equal to those at the
end of inflation, $\rho_{E}^{\rm RH}=\frac{1}{2}E^{2}$ and $\rho_{\chi}^{\rm
RH}=\frac{\bar{\sigma}_{\rm end}}{3}E^{2}$ with $E$ given by Eq. (1). So we
have the hierarchy $\rho_{R}(T_{\rm RH})\gg\rho_{E}^{\rm
RH}\gg\rho_{\chi}^{\rm RH}$.

The expression for the conductivity in Eq. (6) was derived in Kobayashi and
Afshordi (2014) assuming a de-Sitter background so it does not hold in the
epoch of radiation domination (RD) for which, as far as we know, there are no
studies of the conductivity. However, we presume it will maintain the same
exponential suppression during RD of ${\rm Exp}\left[-\pi
m_{\chi}^{2}/g_{D}E\right]$. After reheating the dark electric field (energy
density) redshifts like radiation so the exponential suppression quickly
becomes more severe. We thus expect Schwinger pair production to be negligible
during RD. In this case the energy density of $\chi$ simply redshifts like
non-relativistic matter from reheating until today,

$\rho_{\chi}(T)=\frac{\bar{\sigma}_{\rm end}}{3}E^{2}\left(\frac{T}{T_{\rm
RH}}\right)^{3}\,,$ (11)

while $\rho_{E}$, which redshifts like $T^{4}$, remains negligible compared to
$\rho_{R}$ throughout the history of the universe.

To justify the non-relativistic treatment of $\chi$ we can examine how the
physical momentum $q_{\chi}$ evolves. It will be subject to Hubble expansion
which wants to decrease it and to the electric force which wants to increase
it. This is encapsulated by the evolution equation for $q_{\chi}$,

$\frac{dq_{\chi}}{dt}+Hq_{\chi}=g_{D}E(t)\,.$ (12)

Here $E(t)$ is the physical dark electric field which depends on time in RD.
It is convenient to switch variables from time $t$ to temperature $T$ in terms
of which we have $E(T)=E(T/T_{\rm RH})^{2}$ and the evolution equation
becomes,

$\frac{dq_{\chi}}{dT}-\frac{q_{\chi}}{T}=-\frac{1}{T}\frac{g_{D}E}{H_{\rm
RH}}\,.$ (13)

Here $E$ is a constant which is given by Eq. (1) for the mechanism utilized in
this study. We take as the initial condition $q_{\chi}(T_{\rm
RH})=q_{\chi}^{\rm end}\approx H_{\rm end}$. This is based on the expectation
that the typical momentum of $\chi$ at production will be of the same order as
for the dark photon which generates $E$ and which was shown Bastero-Gil et al.
(2022) to be of order the Hubble scale. The solution to Eq. (13) is,

$\displaystyle q_{\chi}(T)=q_{\chi}^{\rm end}\frac{T}{T_{\rm
RH}}+\frac{g_{D}E}{\epsilon_{R}^{2}H_{\rm end}}\left(1-\frac{T}{T_{\rm
RH}}\right)\,.$ (14)

As the universe expands and cools to temperatures much below the reheat
temperature the momentum tends to,

$q_{\chi}(T\ll T_{\rm RH})\approx H_{\rm
end}\frac{1}{\epsilon_{R}^{2}}\frac{g_{D}E}{H_{\rm end}^{2}}\,,$ (15)

erasing any information about the initial condition $q_{\chi}^{\rm end}$. To
make a suitable cold dark matter candidate, we must require that $\chi$ become
non-relativistic before matter radiation equality. This amounts to the
requirement $q_{\chi}(T=T_{MRE}\ll T_{\rm RH})<m_{\chi}$ which gives the
condition,

$1<\frac{g_{D}E}{H_{\rm end}^{2}}<\epsilon_{R}^{2}\frac{m_{\chi}}{H_{\rm
end}}\,.$ (16)

The first inequality expresses the strong dark electric force condition. We
see that $m_{\chi}>H_{\rm end}\approx q_{\chi}^{\rm end}$ so the dark
electrons are non-relativistic already at production during inflation and
remain so after reheating. This also implies that purely gravitational
production is exponentially suppressed though this doesn’t hold in some
specific models of inflation Chung et al. (2019); Ema et al. (2018). From Eq.
(16) we also obtain,

$\frac{m_{\chi}^{2}}{g_{D}E}>\frac{1}{\epsilon_{R}^{4}}>1\,,$ (17)

which justifies assuming that the conductivity given in Eq. (6) is
exponentially suppressed.

We also briefly argue that the dark electrons do not thermalize with the dark
photons. Since $\chi$ is always non-relativistic, the cross section for
$2\leftrightarrow 2$ processes in the dark sector scales as
$\sigma_{2\leftrightarrow 2}\sim g_{D}^{4}/m_{\chi}^{2}$. The number density
is $n_{\chi}(T)\simeq\rho_{\chi}(T)/m_{\chi}$ with $\rho_{\chi}(T)$ given in
Eq. (11). The interaction rate is
$\Gamma(T)=n_{\chi}(T)\sigma_{2\leftrightarrow 2}$ and inherits the
exponential suppression of the conductivity $\bar{\sigma}_{\rm end}$. It is
then easy to check that $\Gamma(T)/H(T)\ll 1$ for any $T<T_{\rm RH}$ so there
is no thermalization and all of the dark electron production comes from the
inflationary Schwinger effect in contrast to Arvanitaki et al. (2021).

The scenario proposed could also have interesting implications for structure
formation. The dark coulomb interaction must be screened at length scales of
roughly the size of cluster of galaxies or it would be stronger than the
gravitational interaction. Some screening occurs due to the fact that the dark
photon will inevitably get an effective mass (plasma frequency) through the
interaction with the dark electron plasma. We can estimate the effective mass
at redshift around 100 as
$\omega_{p}(T_{*})=g_{D}\sqrt{\rho_{\chi}(T_{*})}/m_{\chi}$ taking
$T_{*}\simeq 100\ T_{0}$. This gives $\omega_{p}(T_{*})\sim g_{D}10^{4}\ {\rm
Mpc}^{-1}$ for $m_{\chi}\sim 10^{14}$ GeV which implies that for a dark gauge
coupling $\gtrsim 10^{-4}$ there should be enough screening. We leave a
detailed study of the impact of our scenario on structure formation to future
work. One should also worry about constraints from the observed triaxality of
dark matter halos Agrawal et al. (2017) and from plasma instabilities Lasenby
(2020), but because the dark electrons in the mechanism considered here are so
heavy we are safe from these bounds.

##  V Dark matter parameter space

Figure 1: Contours of $\Omega_{\chi}/\Omega_{\rm CDM}=1$ with $g_{D}$ the dark
gauge coupling, $E$ the magnitude of the dark electric field at the end of
inflation, and $m_{\chi}$ the mass of the dark (scalar) electron. The Hubble
parameter at the end of inflation is labeled $H_{\rm end}$ while
$\epsilon_{R}$ parametrizes the fraction of the energy density transferred
from the inflaton to the visible radiation at reheating (see Eq. (10)). The
dotted, dashed, and solid black lines show contours where the dark electrons
relic abundance matches the one observed for dark matter,
$\Omega_{\chi}=\Omega_{\rm CDM}$. The light blue region is excluded by the
constraint in Eq. (17).

We have learned that once the dark electrons are produced at the end of
inflation they simply redshift like non-relativistic matter until today. Thus
to see if they can account for the observed dark matter relic abundance we
only have to use Eq. (11) and set $\rho_{\chi}(T_{0})=\rho_{\rm CDM}=9.6\cdot
10^{-48}\ {\rm GeV}^{4}$ with $T_{0}\simeq 2.3\cdot 10^{-13}$ GeV the current
temperature of the universe. We can then write the ratio
$\rho_{\chi}(T_{0})/\rho_{\rm CDM}=\Omega_{\chi}/\Omega_{\rm CDM}$ as,

$\frac{\Omega_{\chi}}{\Omega_{\rm CDM}}=\frac{8.2\cdot
10^{13}}{\epsilon_{R}^{3}}\left(\frac{g_{D}E}{H_{\rm
end}^{2}}\right)^{3}e^{-\pi\frac{m_{\chi}^{2}}{g_{D}E}}\left(\frac{H_{\rm
end}}{10^{14}\ {\rm GeV}}\right)^{5/2}.$ (18)

In Fig. 1 we present contours of $\Omega_{\chi}/\Omega_{\rm CDM}=1$ showing
the parameter space where the dark electrons can constitute the entirety of
dark matter. As can be seen, due to the exponential suppression in Eq. (18),
the final dark matter parameter space available is limited to values of
$m_{\chi}\gtrsim H_{\rm end}\simeq(10^{12}-10^{14})$ GeV. Furthermore, the
constraints of Eq. (16) and Eq. (17), combined with the dark matter abundance
requirement of $\Omega_{\chi}/\Omega_{\rm CDM}=1$, imply $\epsilon_{R}$ is
relatively close to 1. Thus inflationary Schwinger production points to a high
scale of inflation with a dark electron dark matter mass a few times larger
and very efficient reheating. This holds generically for any mechanism which
can generate the necessary background dark electric field though the viable
parameter space can be somewhat larger in specific cases such as the one
explored in Eq. (1) which leads to Bastero-Gil et al. (2019, 2022).

##  VI Summary and Discussion

We have considered a minimal setup to explain the origin of dark matter
consisting of a dark $U(1)_{D}$ gauge group containing a dark photon and dark
(scalar) electrons, $\chi$. If the dark photons form a strong constant
classical electric field during inflation, dark electrons are inevitably
produced via the Schwinger effect. The production occurs in the ‘strong
electric force’ regime where $g_{D}E>H^{2}$. The electric force acts to
accelerate the dark electrons once they are produced. Thus, to ensure a viable
cold dark matter candidate, the dark electrons must be very heavy with
$m_{\chi}\gtrsim H$ implying they are non-relativistic at production. Their
energy density is exponentially suppressed compared to the dark electric field
which in turn is small compared to that stored in the inflaton. However,
because the dark electrons redshift as non-relativistic matter eventually they
come to dominate the energy budget of the universe and can provide a viable
cold dark matter candidate. This requires a high Hubble scale of inflation of
order $H\sim 10^{14}$ GeV. Furthermore, reheating must be efficient meaning
that the inflaton has to transfer roughly a third of its energy density to
radiation in the visible sector. We did not explore mechanisms of reheating,
but it would be interesting to explore those that can achieve a high
efficiency.

The dark matter production mechanism discussed here is non-thermal. The dark
sector does not thermalize with itself and remains decoupled from the thermal
bath of the visible sector as we assumed no coupling between the two sectors.
One could contemplate introducing a kinetic mixing between the dark and
visible (electromagnetic) $U(1)$’s. Since the dark electrons are super heavy,
$m_{\chi}\sim 10^{14}$ GeV one could in principle allow for a large kinetic
mixing and still avoid thermalization as the cross sections and reaction rates
are suppressed by the large $m_{\chi}$. A large kinetic mixing could give some
leverage for possible experimental searches of this type of dark matter. As
discussed above, there could also be interesting implications for structure
formation due to the dark coulomb interaction, but we leave a dedicated study
of these interesting possibilities for future work.

Acknowledgments: The authors thank Antonio Torres Manso, Manel Masip, Jose
Santiago and Takeshi Kobayashi for useful comments and discussions. This work
has been partially supported by Junta de Andalucá Project A-FQM-472-UGR20
(fondos FEDER) (R.V.M.) and by SRA (10.13039/501100011033) and ERDF under
grant PID2022-139466NB-C21(R.V.M.) as well as by MICINN
(PID2019-105943GB-I00/AEI/10.13039/501100011033) (M.B.G.), FCT-CERN grant No.
CERN/FIS-PAR/0027/2021 (M.B.G., P.F.), FCT Grant No.SFRH/BD/151475/2021
(P.F.). This article is based upon work from COST Action COSMIC WISPers
CA21106, supported by COST (European Cooperation in Science and Technology).

## References

  * Bertone and Tait (2018) G. Bertone and M. P. Tait, Tim, Nature 562, 51 (2018), eprint 1810.01668.
  * Graham et al. (2016) P. W. Graham, J. Mardon, and S. Rajendran, Phys. Rev. D93, 103520 (2016), eprint 1504.02102.
  * Agrawal et al. (2020) P. Agrawal, N. Kitajima, M. Reece, T. Sekiguchi, and F. Takahashi, Phys. Lett. B801, 135136 (2020), eprint 1810.07188.
  * Dror et al. (2019) J. A. Dror, K. Harigaya, and V. Narayan, Phys. Rev. D99, 035036 (2019), eprint 1810.07195.
  * Co et al. (2019) R. T. Co, A. Pierce, Z. Zhang, and Y. Zhao, Phys. Rev. D99, 075002 (2019), eprint 1810.07196.
  * Bastero-Gil et al. (2019) M. Bastero-Gil, J. Santiago, L. Ubaldi, and R. Vega-Morales, JCAP 1904, 015 (2019), eprint 1810.07208.
  * Long and Wang (2019) A. J. Long and L.-T. Wang, Phys. Rev. D 99, 063529 (2019), eprint 1901.03312.
  * McDermott and Witte (2020) S. D. McDermott and S. J. Witte, Phys. Rev. D 101, 063030 (2020), eprint 1911.05086.
  * Nakai et al. (2020) Y. Nakai, R. Namba, and Z. Wang (2020), eprint 2004.10743.
  * Ahmed et al. (2020) A. Ahmed, B. Grzadkowski, and A. Socha, JHEP 08, 059 (2020), eprint 2005.01766.
  * Kolb and Long (2021) E. W. Kolb and A. J. Long, JHEP 03, 283 (2021), eprint 2009.03828.
  * Salehian et al. (2021) B. Salehian, M. A. Gorji, H. Firouzjahi, and S. Mukohyama, Phys. Rev. D 103, 063526 (2021), eprint 2010.04491.
  * Bastero-Gil et al. (2022) M. Bastero-Gil, J. Santiago, R. Vega-Morales, and L. Ubaldi, JCAP 02, 015 (2022), eprint 2103.12145.
  * Gorghetto et al. (2022) M. Gorghetto, E. Hardy, J. March-Russell, N. Song, and S. M. West, JCAP 08, 018 (2022), eprint 2203.10100.
  * Sato et al. (2022) T. Sato, F. Takahashi, and M. Yamada, JCAP 08, 022 (2022), eprint 2204.11896.
  * Redi and Tesi (2022) M. Redi and A. Tesi, JHEP 10, 167 (2022), eprint 2204.14274.
  * Barrie (2022) N. D. Barrie (2022), eprint 2211.03902.
  * East and Huang (2022) W. E. East and J. Huang, JHEP 12, 089 (2022), eprint 2206.12432.
  * Arvanitaki et al. (2021) A. Arvanitaki, S. Dimopoulos, M. Galanis, D. Racco, O. Simon, and J. O. Thompson, JHEP 11, 106 (2021), eprint 2108.04823.
  * Anber and Sorbo (2010) M. M. Anber and L. Sorbo, Phys. Rev. D 81, 043534 (2010), eprint 0908.4089.
  * Kobayashi and Afshordi (2014) T. Kobayashi and N. Afshordi, JHEP 10, 166 (2014), eprint 1408.4141.
  * Hayashinaka et al. (2016) T. Hayashinaka, T. Fujita, and J. Yokoyama, JCAP 1607, 010 (2016), eprint 1603.04165.
  * Banyeres et al. (2018) M. Banyeres, G. Domènech, and J. Garriga, JCAP 1810, 023 (2018), eprint 1809.08977.
  * Sobol et al. (2019) O. O. Sobol, E. V. Gorbar, and S. I. Vilchinskii, Phys. Rev. D100, 063523 (2019), eprint 1907.10443.
  * Chung et al. (2019) D. J. H. Chung, E. W. Kolb, and A. J. Long, JHEP 01, 189 (2019), eprint 1812.00211.
  * Ema et al. (2018) Y. Ema, K. Nakayama, and Y. Tang, JHEP 09, 135 (2018), eprint 1804.07471.
  * Agrawal et al. (2017) P. Agrawal, F.-Y. Cyr-Racine, L. Randall, and J. Scholtz, JCAP 05, 022 (2017), eprint 1610.04611.
  * Lasenby (2020) R. Lasenby, JCAP 11, 034 (2020), eprint 2007.00667.