#  Electron beams traversing spherical nanoparticles: analytic and numerical
treatment

P. Elli Stamatopoulou  0000-0001-9121-911X<EMAIL_ADDRESS>POLIMA—Center for
Polariton-driven Light-Matter Interactions, University of Southern Denmark,
Campusvej 55, DK-5230 Odense M, Denmark  Wenhua Zhao∗ 0009-0004-5721-607X
Humboldt-Universit ät zu Berlin, Institut f ür Physik, AG Theoretische Optik
and Photonik, 12489 Berlin, Germany  Max-Born-Institut, 12489 Berlin, Germany
Álvaro Rodríguez Echarri  0000-0003-4634-985X Max-Born-Institut, 12489 Berlin,
Germany  N. Asger Mortensen  0000-0001-7936-6264 POLIMA—Center for Polariton-
driven Light-Matter Interactions, University of Southern Denmark, Campusvej
55, DK-5230 Odense M, Denmark  Danish Institute for Advanced Study, University
of Southern Denmark, Campusvej 55, DK-5230 Odense M, Denmark  Kurt Busch
0000-0003-0076-8522<EMAIL_ADDRESS>Humboldt-Universität zu
Berlin, Institut für Physik, AG Theoretische Optik and Photonik, 12489 Berlin,
Germany  Max-Born-Institut, 12489 Berlin, Germany  Christos Tserkezis
0000-0002-2075-9036 POLIMA—Center for Polariton-driven Light-Matter
Interactions, University of Southern Denmark, Campusvej 55, DK-5230 Odense M,
Denmark  Christian Wolff  0000-0002-5759-6779<EMAIL_ADDRESS>POLIMA—Center for
Polariton-driven Light-Matter Interactions, University of Southern Denmark,
Campusvej 55, DK-5230 Odense M, Denmark

###### Abstract

We present an analytic, Mie-theory based solution for the energy-loss and the
photon-emission probabilities in the interaction of spherical nanoparticles
with electrons passing nearby and through them, in both cathodoluminescence
(CL) and electron energy-loss spectroscopies (EELS). In particular, we focus
on the case of penetrating electron trajectories, for which the complete fully
electrodynamic and relativistic formalism has not been reported as yet. We
exhibit the efficiency of this method in describing collective excitations in
matter through calculations for a dispersive and lossy system, namely a sphere
described by a Drude permittivity, and discuss possible complications when
computing contributions from higher-order modes. Subsequently, we use the
analytic solution to corroborate the implementation of electron-beam sources
in a state-of-the-art numerical method for problems in electrodynamics, the
discontinuous Galerkin time-domain (DGTD) method. We show that the two
approaches produce spectra in good mutual agreement, and demonstrate the
versatility of DGTD via simulations of spherical nanoparticles characterized
by surface roughness. The possibility of simultaneously employing both kinds
of calculations (analytic and numerical) facilitates a better understanding of
the rich optical response of nanophotonic architectures excited by fast
electron beams.

##  I Introduction

In recent decades, electron-beam spectroscopy has emerged as a revolutionary
tool for the optical characterization of materials. Swift electrons passing in
close proximity or through a specimen, undergo energy loss owing to energy
transfer to the optical modes sustained in the material García de Abajo
(2010). From localized and propagating surface plasmons in metallic structures
Watanabe (1956); Ritchie (1957); Yamamoto  _et al._ (2001); Vesseur  _et al._
(2007), to Mie resonances in dielectric resonators Coenen  _et al._ (2013);
Matsukata  _et al._ (2019); Fiedler  _et al._ (2022) and phonon polaritons in
polar crystals Lagos  _et al._ (2017); Hage  _et al._ (2019); Maciel-Escudero
_et al._ (2020), electron-beam spectroscopy has proven quintessential for
mapping collective excitations in a broad spectral range that spans from
ultraviolet to far-infrared frequencies.

With the diffraction limit ultimately being controlled by the de Broglie
wavelength, highly energetic electrons are excellent probes to study the
optical properties of truly nanoscale structures, with atomic spatial
resolution and sub-meV energy resolution Polman  _et al._ (2019); García de
Abajo and Di Giulio (2021). In electron energy-loss spectroscopy (EELS), the
sample is excited by a high-energy ($30-300$ keV) electron beam and the energy
lost to the interaction is measured in a transmission electron microscope
(TEM) setup Egerton (2011); Wu  _et al._ (2018). EELS allows, thereby, the
detection of both radiative and dark modes, including longitudinal bulk
plasmons (BPs) Batson and Silcox (1983), breathing modes Schmidt  _et al._
(2012), or antibonding modes in nanoparticle (NP) dimers Koh  _et al._ (2009).
Optical excitations in thick samples can be imaged in cathodoluminescence (CL)
spectroscopy, performed in scanning electron microscopes (SEMs) at
intermediate beam energies ($1-50$ keV) Polman  _et al._ (2019); García de
Abajo (2010). In CL measurements, the signal collected is the result of far-
field photon emission from the sample, originating from the radiative decay of
the excited modes. Recent advances in instrumentation have even added temporal
resolution in EEL and CL spectra, introducing the field of ultrafast electron
microscopy (UEM) Barwick  _et al._ (2008, 2009).

Considering the recent progress in electron spectroscopy techniques, robust
analytic and computational tools are evidently required to interpret the
plethora of experimental data. While first theoretical efforts were performed
within the non-retarded approximation for the description of plasmons in thin
films Ritchie (1957), the theory was gradually generalized to account for
collective excitations in diverse media and geometries Ferrell and Echenique
(1985); Garcia-Molina  _et al._ (1985), also considering retardation effects
Rivacoba  _et al._ (1992); Zabala  _et al._ (1989). Quantum approaches
Echenique  _et al._ (1987); Ritchie and Howie (1988) and analytic solutions
including relativistic effects for simple geometries were later developed
Kröger (1970); García de Abajo (1999), allowing the combination of high-
velocity electron beams with both common and less conventional materials,
including dielectric media, polar crystals, graphene and other two-dimensional
(2D) materials Fiedler  _et al._ (2022); Maciel-Escudero  _et al._ (2020);
Hage  _et al._ (2020); Yankovich  _et al._ (2019). Most theoretical
relativistic descriptions have focused on aloof electron trajectories, that do
not penetrate the specimen, in contrast to experimental practices, where the
electron beam is typically scanned over the entire sample area Egerton (2011).
However, aloof electron trajectories oftentimes fail to capture intriguing
phenomena associated with bulk properties, such as BPs and bulk phonons Koh
_et al._ (2009); Lagos  _et al._ (2017), and other sources of electron-induced
photon emission, like Cherenkov or transition radiation Kröger (1970);
Pogorzelski and Yeh (1973); Ginzburg (1996); Yamamoto  _et al._ (1996).

Despite the undeniable advantage of analytic solutions in data analysis, their
applicability is limited to a handful of highly symmetric geometries. Over the
years, different numerical schemes have been consolidated to complement
analytic approaches in simulating the electromagnetic properties of
nanophotonic systems of diverse shapes and forms Gallinet  _et al._ (2015),
such as the boundary element method (BEM) García de Abajo and Howie (2002);
Hohenester (2014), the finite element method (FEM) Pomplun  _et al._ (2007);
Burger  _et al._ (2012), or the finite-difference frequency-domain (FDFD) and
finite-difference time-domain (FDTD) methods Das  _et al._ (2012); Cao  _et
al._ (2015). An alternative route is offered by the discontinuous Galerkin
time-domain (DGTD) method, which employs the Galerkin scheme to solve
Maxwell’s equations in the time domain Hesthaven and Warburton (2007);
Matyssek  _et al._ (2011); Busch  _et al._ (2011); Husnik  _et al._ (2013);
Prämassing  _et al._ (2021). This method combines the flexible space
discretization of finite elements, with the memory efficiency and the ability
to include nonlinearities, characteristic of time-domain methods. As a
consequence, DGTD offers great versatility in simulating objects of complex
geometry and nonlinear response.

In this work, we present and compare an analytic approach and the DGTD method
for the study of spherical nanostructures excited by aloof and penetrating
electron beams. Following the work of García de Abajo for aloof electron beams
García de Abajo (1999), we derive analytic formulas for the energy loss and
photon emission probability, generalized here to account for penetrating
trajectories. We then validate the implementation of electron-beam excitation
of nanostructures in DGTD Matyssek  _et al._ (2011) by comparing the EEL and
CL spectra produced by the two methods for a perfectly spherical plasmonic NP
featuring localized surface plasmons (LSPs) and BPs. Finally, we apply the
numerical method to study the optical response of a NP with surface roughness,
showcasing the ability of the DGTD method to emulate scenarios aligned with
realistic experimental conditions that involve imperfect structures Maradudin
and Mills (1975); Maradudin and Méndez (2007).

##  II Methods

###  II.1 Analytic approach

As a first step to examine the agreement and complementarity of analytic and
numerical tools, we outline the modeling of the physical system and the
assumptions made in each method. As a testbed, we consider a perfectly
spherical metal NP of radius $R$ suspended in air. The NP is characterized by
a unity relative permeability ($\mu=1$), while its relative permittivity
depends on the angular frequency $\omega$ as described by the Drude model

$\varepsilon(\omega)=1-\frac{\omega^{2}_{\mathrm{p}}}{\omega(\omega+i\tau^{-1})},$
(1)

with plasma frequency $\omega_{\mathrm{p}}$ and damping rate $\tau^{-1}$.

In the analytic calculation, the electron beam is modeled as a single point
particle that carries the elementary charge $-e$. We assume that it travels
with constant velocity ${\mathbf{v}}$, say along the $z$-axis, following thus
a trajectory ${\mathbf{r}}_{e}={\mathbf{r}}_{0}+{\mathbf{v}}t$, where
${\mathbf{r}}_{0}=(b,\phi_{0},z=-\infty)$ is its initial position in
cylindrical coordinates. The impact parameter $b$ indicates the distance
between the electron trajectory and the center of the sphere, as shown in the
schematics of Fig. 1. In passing, we mention that the angular coordinate
$\phi_{0}$ need not be specified, as it does not enter the calculation, due to
the symmetry of the problem.

Figure 1: Schematic illustration of the geometry under study; a metallic NP of
radius $R=75$ nm, and permittivity described by the Drude model of Eq. (1), is
excited by an electron beam passing with velocity ${\mathbf{v}}$ at impact
parameter $b$ with respect to its center. The color dots correspond to the
selection of impact parameters examined in our study. Labels $\mathrm{I}$ and
$\mathrm{II}$ indicate the region inside the NP and the surrounding medium
(air), respectively.

The fast electron drives plasmon oscillations in the metal, generating an
induced electric and magnetic field, ${\mathbf{E}}_{\mathrm{ind}}$ and
${\mathbf{H}}_{\mathrm{ind}}$ respectively, eventually resulting in the
emission of radiation to the environment (labeled as region $\mathrm{II}$ in
Fig. 1). Following the basic steps of Mie theory Mie (1908), one can decompose
the induced electromagnetic field into transverse electric (TE) and transverse
magnetic (TM) components. In region $\mathrm{II}$ the induced fields take the
general form of outgoing spherical waves Jackson (1998)

$\displaystyle{\mathbf{E}}_{\mathrm{ind}}^{\mathrm{II}}({\mathbf{r}})$
$\displaystyle=\sum_{\ell=1}^{\infty}\sum_{m=-\ell}^{+\ell}\Big{\\{}b_{\ell
m}^{\mathrm{II}}h^{+}_{\ell}(k_{0}r){\mathbf{X}}_{\ell
m}({\bm{\hat{\mathbf{r}}}})$ $\displaystyle\qquad\qquad+\frac{i}{k_{0}}a_{\ell
m}^{\mathrm{II}}\bm{\nabla}\times h^{+}_{\ell}(k_{0}r){\mathbf{X}}_{\ell
m}({\bm{\hat{\mathbf{r}}}})\Big{\\}}$ (2a)
$\displaystyle{\mathbf{H}}_{\mathrm{ind}}^{\mathrm{II}}({\mathbf{r}})$
$\displaystyle=\frac{1}{Z_{0}}\sum_{\ell=1}^{\infty}\sum_{m=-\ell}^{+\ell}\Big{\\{}a_{\ell
m}^{\mathrm{II}}h_{\ell}^{+}(k_{0}r){\mathbf{X}}_{\ell
m}({\bm{\hat{\mathbf{r}}}})$ $\displaystyle\qquad\qquad-\frac{i}{k_{0}}b_{\ell
m}^{\mathrm{II}}\bm{\nabla}\times h_{\ell}^{+}(k_{0}r){\mathbf{X}}_{\ell
m}({\bm{\hat{\mathbf{r}}}})\Big{\\}},$ (2b)

where $\ell$, $m$ are the angular momentum quantum numbers, and $a/b_{\ell
m}^{\mathrm{II}}$ denote the expansion coefficients of the TE/TM components
corresponding to modes of electric/magnetic multipole character (see the
appendix for the analytic expressions). Furthermore, in Eqs. (2)
$h_{\ell}^{+}$ is the spherical Hankel function of the first kind,
${\mathbf{X}}_{\ell m}$ are the vector spherical harmonics, $k_{0}=\omega/c$,
is the wave number in free space, $c=1/\sqrt{\varepsilon_{0}\mu_{0}}$ the
speed of light in vacuum, where $\varepsilon_{0}$ and $\mu_{0}$ denote the
vacuum permittivity and permeability, respectively, and
$Z_{0}=\sqrt{\mu_{0}/\varepsilon_{0}}$ is the impedance in free space. Details
about the field expansions can be found in section S.I of the Supporting
Information (SI).

The energy radiated in the far field can be found by integrating the Poynting
flux at a spherical surface of radius $r\to\infty$, in the normal direction
${\bm{\hat{\mathbf{r}}}}$. Then the (CL) probability of collecting a photon of
energy $\hbar\omega$ is given by

$\Gamma_{\mathrm{CL}}(\omega)=\frac{r^{2}}{\pi\hbar\omega}\int
d\Omega\,\mathrm{Re}\big{\\{}{\mathbf{E}}_{\mathrm{ind}}^{\mathrm{II}}({\mathbf{r}},\omega)\times{\mathbf{H}}_{\mathrm{ind}}^{\mathrm{II}^{\ast}}({\mathbf{r}},\omega)\big{\\}}\cdot{\bm{\hat{\mathbf{r}}}},$
(3)

where $d\Omega$ denotes the infinitesimal solid angle. By inserting Eqs. (2)
into Eq. (3), and evaluating the result in the far field ($k_{0}r\to\infty$),
we find

$\Gamma_{\textrm{CL}}(\omega)=\frac{1}{\pi\hbar\omega
Z_{0}k^{2}_{0}}\sum_{\ell=1}^{\infty}\sum_{m=-\ell}^{+\ell}\Big{\\{}\big{|}b_{\ell
m}^{\textrm{II}}\big{|}^{2}+\big{|}a_{\ell
m}^{\textrm{II}}\big{|}^{2}\Big{\\}}.$ (4)

Apart from the emission of radiation, part of the energy transferred to the
optical modes of the NP dissipates non-radiatively, owing to the intrinsic
losses within the material. The total energy lost can be calculated by the
work done by the electron against the induced field along the entire electron
trajectory. Then the (EEL) probability of the electron losing energy
$\hbar\omega$ is given by

$\Gamma_{\mathrm{EELS}}(\omega)=\frac{e}{\pi\hbar\omega}\int
dt\,\mathrm{Re}\big{\\{}\exp{(-i\omega
t)}\,{\mathbf{v}}\cdot{\mathbf{E}}_{\mathrm{ind}}({\mathbf{r}}_{e},\omega)\big{\\}}.$
(5)

The integral in Eq. (5) can be decomposed into three terms (see details in
section S.I of the SI)

$\Gamma_{\mathrm{EELS}}(\omega)=\Gamma_{\mathrm{bulk}}(\omega)+\Gamma_{\mathrm{surf}}(\omega)+\Gamma_{\mathrm{Begr}}(\omega).$
(6)

Here, $\Gamma_{\mathrm{bulk}}$ is related to the bulk modes of the unbound
medium, reduced by the Begrenzung term $\Gamma_{\mathrm{Begr}}$ that accounts
for the presence of a boundary Lucas and Šunjić (1972). The
$\Gamma_{\mathrm{surf}}$ term contains the contribution from modes excited by
the part of the electron trajectory lying externally to the NP, and is, thus,
associated with the excitation of LSPs. The terms entering Eq. (6) are given
by the following formulas

$\Gamma_{\mathrm{bulk}}(\omega)=\frac{e^{2}z_{e}}{2\pi^{2}\varepsilon_{0}\hbar
v^{2}}\mathrm{Im}\Bigg{\\{}\frac{1}{\gamma_{0}^{2}}\mathrm{ln}\left(\left[\frac{q_{\mathrm{c}}\gamma_{0}v}{\omega}\right]^{2}+1\right)-\frac{1}{\gamma^{2}\varepsilon}\mathrm{ln}\left(\left[\frac{q_{\mathrm{c}}\gamma
v}{\omega}\right]^{2}+1\right)\Bigg{\\}},$ (7a)
$\displaystyle\Gamma_{\mathrm{surf}}(\omega)=\frac{e}{\pi\hbar\omega}\mathrm{Re}\sum_{\ell=1}^{\infty}\sum_{m=-\ell}^{+\ell}$
$\displaystyle\Bigg{\\{}\frac{K_{m}\left(\omega
b/[v\gamma_{0}]\right)}{ik_{0}\sqrt{\ell(\ell+1)}}\left[mb_{\ell
m}^{\mathrm{II}}\mathcal{M}_{\ell m}^{*}-a_{\ell
m}^{\mathrm{II}}\frac{\mathcal{N}_{\ell m}^{*}}{\beta\gamma_{0}}\right]$
$\displaystyle-$
$\displaystyle\int_{-z_{\mathrm{e}}}^{z_{\mathrm{e}}}dz\,\frac{\exp{(-i\omega
z/v)}}{\sqrt{\ell(\ell+1)}}\left[mb_{\ell
m}^{\mathrm{II}}h_{\ell}^{+}(k_{0}r){Y^{m}_{\ell}}\left(\theta,0\right)-\frac{a_{\ell
m}^{\mathrm{II}}}{k_{0}b}\big{\\{}\mathcal{H}_{\ell
m}^{+}(k_{0}z)+\mathcal{H}_{\ell m}^{-}(k_{0}z)\big{\\}}\right]\Bigg{\\}},$
(7b) and
$\Gamma_{\mathrm{Begr}}(\omega)=\frac{e}{\pi\hbar\omega}\mathrm{Re}\sum_{\ell=1}^{\infty}\sum_{m=-\ell}^{+\ell}\int_{-z_{\mathrm{e}}}^{z_{\mathrm{e}}}dz\frac{\exp{(-i\omega
z/v)}}{\sqrt{\ell(\ell+1)}}\left[mb_{\ell
m}^{\mathrm{I}}j_{\ell}(kr){Y^{m}_{\ell}}\big{(}\theta,0\big{)}-\frac{a_{\ell
m}^{\mathrm{I}}}{kb}\big{\\{}\mathcal{J}_{\ell m}^{+}(kz)+\mathcal{J}_{\ell
m}^{-}(kz)\big{\\}}\right].$ (7c)

Here, $2z_{\mathrm{e}}=2\sqrt{R^{2}-b^{2}}$ is the length of the electron path
inside the NP, and $\gamma=1/[1-\varepsilon\beta^{2}]^{1/2}$ with $\beta=v/c$
are the Lorentz kinematic factors ($\gamma_{0}$ is evaluated in free space).
In Eqs. (7) and (7c) $K_{m}$ is the modified Bessel function of the second
kind and $Y_{\ell}^{m}$ are the spherical harmonics. In addition, we have set
$r=\sqrt{b^{2}+z^{2}}$, and $\theta=\arccos(z/r)$, while analytic expressions
for coefficients $a/b_{\ell m}^{\mathrm{I}}$, $\mathcal{M}_{\ell m}$,
$\mathcal{N}_{\ell m}$, $\mathcal{H}_{\ell m}^{\pm}$, and $\mathcal{J}_{\ell
m}^{\pm}$ can be found in the appendix.

It is important to note here that the aforementioned decomposition of the EEL
spectra introduces a free parameter; assuming that upon losing energy
$\hbar\omega$ the electron transfers a transverse (with respect to the
electron trajectory) momentum $q$ to excite an optical mode, $q_{\mathrm{c}}$
is the maximum transverse momentum collected. In an experiment, the momentum
cutoff is determined by the half-aperture collection angle $\varphi$ of the
microscope spectrometer, as

$\hbar q_{\mathrm{c}}\approx\sqrt{(m_{e}v\varphi)^{2}+(\hbar\omega/v)^{2}},$
(8)

where $m_{e}$ is the electron mass. We may freely choose the value for this
momentum cutoff, making sure that it aligns with the typical values for the
collection angle in scanning TEM (STEM) setups, which are in the order of a
few mrad. Naturally, this introduces a level of arbitrariness in the EEL
spectra, as various values of $q_{\mathrm{c}}$ lead to different peak
intensities at the BP energy.

###  II.2 DGTD simulation

We complement our analytical work with numerical simulations of the
electromagnetic problem of a NP excited by a moving Gaussian charge
distribution. To this end, we employ the DGTD method (see details in section
S.II of the SI), which combines a piecewise polynomial spatial interpolation
on an unstructured tetrahedral mesh with a Runge–Kutta time integrator to
obtain a high–order accurate explicit solver for Maxwell’s equations in time
domain

$\displaystyle\partial_{t}{\mathbf{H}}({\mathbf{r}},t)=$
$\displaystyle-\mu^{-1}_{0}\mu^{-1}({\mathbf{r}})\,\bm{\nabla}\times{\mathbf{E}}({\mathbf{r}},t),$
(9a) $\displaystyle\partial_{t}{\mathbf{E}}({\mathbf{r}},t)=$
$\displaystyle\,\varepsilon^{-1}_{0}\varepsilon^{-1}({\mathbf{r}})\,\left[\bm{\nabla}\times{\mathbf{H}}({\mathbf{r}},t)-{\mathbf{j}}({\mathbf{r}},t)\right].$
(9b)

Here, ${\mathbf{j}}$ is the total current density that encompasses both any
current associated with the excitation source, as well as dispersive
polarization currents. The resulting method is memory–efficient compared to
traditional finite elements and especially well–suited for the calculation of
wide–band spectra.

One key difference between the Mie-based theory and the DGTD simulations,
which can potentially lead to deviations between the two approaches, is the
implementation of the excitation source. In the numerical treatment, we model
the electron beam with a Gaussian charge distribution of the form

$\rho({\mathbf{r}})=-\frac{e}{\sigma_{e}^{3}\sqrt{\pi^{3}}}\exp(-r^{2}/\sigma_{e}^{2}),$
(10)

with width $\sigma_{e}=5$ nm. This choice essentially prevents numerical
artifacts arising when implementing a point-charge particle moving inside the
simulation domain, while also being compatible with the typical spot size in
CL experiments Fiedler  _et al._ (2022). Thereby, we introduce a new spatial
scale, which needs to be taken into account when it is comparable to the
parameters related to the discretization of the computational domain, such as
the mesh element size, as well as the characteristic lengths of the physical
system, e.g. the radius of the NP and the impact parameter. For very large
distances between the electron and the NP (i.e., $b-R\gg\sigma_{e}$), the
source resembles a point charge and we, therefore, expect an excellent
agreement with analytic results. However, in the opposite scenario, the finite
width of the electron beam becomes important, and the corresponding fields do
not accurately match those of a point charge.

##  III Results and discussion

###  III.1 EEL and CL spectroscopy of perfectly spherical NPs

Figure 2:  (a) Scattering cross section, normalized to the geometrical cross
section $\pi R^{2}$, of the metallic NP of Fig. 1 under plane-wave excitation
(gray shaded area). The spectrum is decomposed into the contributions from the
first three multipoles (blue, red and green curves, for $\ell=1,2,3$
respectively). (b) CL probability (yellow shaded area) decomposed into the
contributions from the first four multipoles (blue, red, green, and violet
curves). (c) EEL probability (yellow shaded area) decomposed into the surface,
bulk, and Begrenzung contributions (blue, green and pink curves), with the
inset focusing into the energy window $1-4$ eV. The bulk term is given by Eq.
(7a) for $q_{\mathrm{c}}=0.71$ nm-1. In (b) and (c) we consider an electron
traveling with velocity $v=0.33c$ (kinetic energy $\approx 30$ keV) passing
through the NP of Fig. 1 at impact parameter $b=35$ nm.

In what follows, we analyze the response of a metallic NP excited by a fast
electron beam, as predicted by both the analytic and the numerical approach.
We consider a smooth sphere of radius $R=75$ nm with plasma energy
$\hbar\omega_{\mathrm{p}}=5$ eV and damping rate $\tau^{-1}$ corresponding to
an energy of $\hbar\tau^{-1}=50$ meV. The parameters mimic typical plasmonic
metals and are chosen for the purpose of illustration, while the particular
values have no consequences for our general conclusions. Fig. 2 shows the
scattering spectrum (panel a) under plane-wave illumination, and the CL and
EELS probability (panels b and c, respectively) calculated for a low electron
velocity $v=0.33c$ (kinetic energy $\approx 30$ keV) intersecting the NP at
$b=35$ nm from its center. As suggested by Eq. (4), the total photon emission
probability can be decomposed into contributions of pairwise orthogonal
electric and magnetic multipoles of order $\ell$. Fig. 2b reveals that the CL
spectrum is composed of the contributions of the first four ($\ell=1,2,3,4$)
electric-type modes, appearing at approximately $2$, $2.8$, $3.1$, and $3.2$
eV, associated with the excitation of LSP resonances, while higher-order
($\ell>4$) multipoles contribute negligibly to the spectrum. In comparison
with the scattering cross section of the NP shown in Fig. 2a (calculated using
Mie theory Bohren and Huffman (1983)), the CL spectrum provides very similar
information. This is somewhat expected, since both calculations are based on
the collection of far-field radiation. We observe, nonetheless, two notable
distinct features in the CL spectrum of Fig. 2b. Firstly, the electron source
excites more efficiently the $\ell=2$ and $3$ modes, whereas in the scattering
spectrum of Fig. 2a the dipolar mode peak features the highest intensity. The
relative peak intensities in the CL and EEL spectra depend strongly on the
impact parameter, which determines the arrangement of the polarization charges
in the material Batson  _et al._ (2011). Secondly, we observe a small redshift
of the dipolar ($\ell=1$) mode in CL due to retardation, stemming from the
fact that the speed of the electron is only a fraction of the speed of light.

Figure 3: (a) CL and (b) EEL probability in the interaction between the
metallic NP of Fig. 1 and an electron traveling with velocity $v=0.33c$
(kinetic energy $\approx 30$ keV) considering 5 different impact parameters,
as denoted in the labels and illustrated in the insets of each panel, and
momentum cutoff $q_{\mathrm{c}}=0.71\,\mathrm{nm}^{-1}$. The shaded areas
correspond to results of the analytic Mie-type calculation, whereas the solid
lines show the corresponding DGTD simulations. The vertical gray dotted and
dashed lines in the three leftmost panels of (b) trace the energy of the SPP
($\hbar\omega_{\mathrm{spp}}\approx 3.5$ eV) and the BP mode
($\hbar\omega_{\mathrm{bp}}=5$ eV), respectively.

In Fig. 2c we present the EEL spectrum, decomposed as described in Eq. (6). In
the zoom-in area of the figure we observe the excitation of numerous higher-
order multipoles, appearing as sharp peaks at energies up to $3.5$ eV. With
increasing multipole order, the wavelength of the corresponding mode reduces,
so higher-order modes experience the curved surface of the NP as increasingly
more flat. As a result, they accumulate at the energy corresponding to that of
a surface plasmon polariton (SPP) at a planar interface, at
$\hbar\omega_{\mathrm{spp}}=\hbar\omega_{\mathrm{p}}/\sqrt{2}\approx 3.5$ eV.
Above the SPP energy, the spectrum exhibits a pronounced peak at the BP energy
$\hbar\omega_{\mathrm{bp}}=\hbar\omega_{\mathrm{p}}=5$ eV, pertaining to the
excitation of BPs in the volume of the NP. Since BPs are longitudinal modes,
they do not couple to far-field radiation and, therefore, they can be detected
only in EELS. At the same energy, we observe the expected negative peak
related to the Begrenzung term, reducing the BP peak in the total EEL
probability Lucas and Šunjić (1972).

Having a clear picture of the origin of all spectral features, in Figs. 3a and
b we compare the CL and EEL probability, respectively, of the same metal NP,
as calculated employing the analytic (Mie) and the DGTD method. We test the
agreement between the two calculations probing various impact parameters, that
range from $b=125$ nm to $10$ nm, corresponding to aloof electron trajectories
(violet and blue spectra), grazing (green), and penetrating (yellow and red).
Overall, the DGTD method reproduces the positions of the LSP modes and the
corresponding CL probabilities of the Mie calculations. In Fig. 3a we show an
excellent agreement in the CL spectra, with a relative error at around $1\%$
for aloof and penetrating electron trajectories (see Table 1). The highest
error is acquired when the electron grazes the surface of the NP, passing
exactly at $b=R$ (middle panel in Fig. 3a). This point reflects an important
limit in the capabilities of the DGTD method; in the grazing trajectory, due
to the finite width of the electron beam, half of the Gaussian charge density
distribution lies inside the NP, while the other half lies outside, leading to
numerical inconsistencies. One may avoid this point, which is inevitably
difficult to resolve, by slightly adjusting the impact parameter by half the
Gaussian width.

Regarding compatibility of the EEL spectra, the two rightmost panels in Fig.
3b reveal an excellent agreement for aloof electron trajectories. We
consistently find a higher error for grazing trajectories (middle panel in
Fig. 3b); here, the EELS spectrum calculated with DGTD exhibits a numerical
artifact at the BP energy (gray dashed line), resulting once again from the
fact that only a fraction of the electron charge density distribution
penetrates the NP, exciting only partially the bulk mode. A substantial
deviation between the two methods is found in EEL spectra for both grazing and
penetrating electron trajectories (three leftmost panels in Fig. 3b) between
the SPP and the BP energy, denoted by the gray dotted and dashed lines,
respectively. The disagreement is, naturally, reflected in the large relative
errors presented in Table 1 accordingly. At these impact parameters, the
condition $b-R\gg\sigma_{e}$ is not fulfilled, hence the field of the electron
deviates considerably from that of a point charge. Moreover, the beam width
$\sigma_{e}$ becomes important compared to the mesh element size, since the
surface of the NP is more finely discretized than the surrounding medium (see
Fig. S4 in the SI). Finally, there exists an additional source of error in the
evaluation of the BP contribution, that stems from the rather arbitrary choice
of the transverse momentum cutoff $q_{\mathrm{c}}$ in the analytic approach.
In contrast, in the DGTD implementation there is a respective internal limit,
associated with $\sigma_{e}$. As a result, we consistently find higher
relative errors in EELS in comparison with CL, and for grazing and penetrating
electron trajectories as compared to aloof, probing both low-energy electron
beams as in Fig 3, as well as higher energies that are more realistic for EEL
measurements (see Fig. S.8 in the SI for CL and EEL spectra at energy $200$
keV).

Table 1:  Relative error [see Eq. (S.38) in the SI] between the analytic and the DGTD calculations of $\Gamma_{\mathrm{CL}}$ and $\Gamma_{\mathrm{EELS}}$ for varying $b$. $b$ (nm) |  $\Gamma_{\mathrm{CL}}$ (%) |  $\Gamma_{\mathrm{EELS}}$ (%)
---|---|---
10 | 1.26 | 9.75
35 | 1.06 | 6.17
75 | 8.23 | 28.87
100 | 1.06 | 1.02
125 | 1.07 | 1.93

Admittedly, the accumulation point of high-order multipoles at
$\hbar\omega_{\mathrm{spp}}$ is hard to be resolved in both methods. On the
one hand, the analytic Mie calculation assumes a moving point charge, which
can, in principle, excite an infinite number of multipoles, resulting in a
sharp high-intensity peak at energy $\hbar\omega_{\mathrm{spp}}$. On the other
hand, the finite mesh size and beam width implemented in DGTD imposes a
limitation to the number of multipoles that can be resolved for a given
discretization, since high order multipoles associated with field variation
shorter than the mesh element size at the surface cannot be captured without
the use of very high order polynomials. In EELS and CL experiments, there
exists an analogous limitation, associated with the finite width of the
electron beam employed, as well as the geometric imperfections of the NP. The
versatility of DGTD allows us not only to adjust the electron beam width
according to the experimental setup, but also to mimic NPs with surface
roughness, as we discuss in section III.2. In the analytic approach, the
smearing of higher-order modes and the overall quenching of the sharp peak at
the accumulation point can too be reproduced, once we consider the nonlocal
response of the material; this can be done particularly easily within Mie
theory Christensen  _et al._ (2014); Zouros  _et al._ (2020); Gonçalves and
García de Abajo (2023). Nonlocal effects manifest as increased damping and
uneven energy shifts of high-order modes, and, therefore, lead to the
suppression of the individual modes, as well as the reduction of their overlap
at $\hbar\omega_{\mathrm{spp}}$ Moeferdt (2017).

Figure 4:  (a, c) Area under the curve of the EEL spectra of Fig. 3b, (a) for
the aloof electron trajectory at $b=125$ nm, (c) for the grazing trajectory at
$b=75$ nm, showcasing the convergence of the EEL probability for increasing
values of the multipole cutoff $\ell_{\mathrm{max}}$, plotted versus
$1/\sqrt{\ell_{\mathrm{max}}}$. In both panels, the magenta line is fitted to
the data points marked as black bullets, while the red crosses represent data
points excluded from the fitting. The vertical gray lines serve as guides to
the eye for the position of $\ell_{\mathrm{max}}=1,2\ldots 5$. (b, d) EEL
probability calculated for impact parameters (b) $b=125$ nm, and (d) $b=75$
nm, and for selected values of $\ell_{\mathrm{max}}$, as denoted in the
labels.

The difficulty in resolving the high-order multipoles even in the analytic
calculation is clear in the convergence study presented in Figs. 4 and 5. Due
to multipole orthogonality, increasing $\ell_{\mathrm{max}}$ only ever adds
signal. This suggests monitoring the total EEL spectrum integrated over energy
(area under the curve) as a proxy for convergence. For aloof trajectories, the
asymptotics of the Hankel functions suggest exponential convergence, which is
in agreement with Figs. 4a and b; the EEL probability converges at
$\ell_{\mathrm{max}}=10$. In contrast, Fig. 4c shows that in the case of the
grazing trajectory the convergence order breaks down to a square root law. By
extrapolation (magenta line) we can conclude that even for
$\ell_{\mathrm{max}}=63$ the analytic result is still converged only up to
about $15\%$. Fig. 4d corroborates that the area missing from the converged
value corresponds to the higher-order modes piling up at the SPP energy.
Finally, we note that Fig. 4a exhibits the same square root convergence before
the curve flattens off. The exact value of $\ell_{\mathrm{lmax}}$ where this
transition happens increases as the impact parameter approaches $R$.

Figure 5: (a, c) Area under the curve of the EEL spectrum of the leftmost
panel in Fig. 3b ($b=10$ nm), showcasing the divergence of the EEL probability
for increasing values of (a) the multipole cutoff $\ell_{\mathrm{max}}$, and
(c) the transverse momentum cutoff $q_{\mathrm{c}}$ plotted versus
$\mathrm{ln}(q_{\mathrm{c}})$. The data follow a linear divergence trend with
respect to their corresponding horizontal axes. In both panels, the magenta
line is fitted to the data points marked as black bullets, while the red
crosses represent data points excluded from the fitting. (b, d) EEL
probability calculated at $b=10$ nm for selected values of (b)
$\ell_{\mathrm{max}}$, and (d) $q_{\mathrm{c}}$. In panels (a, b) we scan over
$\ell_{\mathrm{max}}$, while keeping a fixed value $q_{\mathrm{c}}=0.71$ nm-1,
whereas in (c, d) we scan over $q_{\mathrm{c}}$ for a fixed value
$\ell_{\mathrm{max}}=63$, as denoted in the labels. In panel (c) the vertical
lines correspond to the selected $q_{\mathrm{c}}$ values presented in (d),
following the same color coding.

For the penetrating trajectories shown in Fig. 5 we find signs of the same
slow convergence around the SPP accumulation point. This is masked in our
convergence plots by the fact that the BP peak diverges. The source of the
divergence lies in the decomposition of the EEL probability into two competing
contributions stemming from the bulk and the Begrenzung terms. As illustrated
in Figs. 5a and c, the two terms produce divergences of opposite sign; the
negative Begrenzung term diverges linearly for increasing multipole order
$\ell_{\mathrm{max}}$, whereas the positive bulk term diverges logarithmically
for increasing momentum cutoff $q_{\mathrm{c}}$ [see Eqs. (7a) and (7c)]. As
our computational resources do not allow driving $\ell$ and $q_{\mathrm{c}}$
to infinity, once one of the two parameters is truncated to a certain cutoff
value, the other has to be adjusted accordingly. Figs. 5b and d show that the
different values of $\ell$ and $q_{\mathrm{c}}$, respectively, affect the
spectra in the energy window between the SPP and BP modes.

###  III.2 EELS of NPs with surface roughness

The synthesis of spherical metallic NPs, as the ones studied in the present
work, is routinely done with colloidal chemistry, for a large variety of
materials and NP shapes Link and El-Sayed (1999); Grzelczak  _et al._ (2010).
However, despite being able to accurately control the NP size, assuring a
smooth surface is rather challenging. Typically the structures exhibit
protuberances on the surface, which can be responsible for symmetry breaking
Nappa  _et al._ (2005); Salomon  _et al._ (2013), hot spots in dimers Zheng
_et al._ (2015); Yoon  _et al._ (2019) and picocavities Benz  _et al._ (2016);
Urbieta  _et al._ (2018), or energy shifts of the LSPs Rodríguez-Fernández
_et al._ (2009). Within the DGTD method, surface texture can be easily
implemented on top of the perfect spherical mesh and incorporated in the
numerical calculations. Here, we follow the prescription presented in Ref.
Loth  _et al._ (2023) to implement the desired roughness as a radius variation
derived from Gaussian white noise, with a correlation length of $10$ nm and
two different values for the root-mean square (rms) amplitude (see details in
section S.II.F of the SI).

In Fig. 6 we explore the effect of surface roughness on a spherical NP,
described by the same Drude permittivity as in the previous section, now
excited by an electron beam traveling with velocity $v=0.7c$ (kinetic energy
$\approx 200$ keV) at distance $b=100$ nm from its center. We probe meshes of
two degrees of surface roughness on top of the NP of nominal radius $R=75$ nm,
namely $\mathrm{rms}=2$ nm and $4$ nm. Since the breaking of the spherical
symmetry introduces a dependence on the electron propagation direction, and on
the mesh morphology, in Fig. 6 we plot the _average_ EEL probability,
corresponding to the average values obtained for $6$ different meshes. The
resulting spectra for $\mathrm{rms}=2$ nm and $4$ nm (dark red and blue
curves, respectively) deviate notably from that of a smooth sphere (gray
shaded area). Firstly, we observe an increasing redshift of the spectra with
increasing degree of roughness, in agreement with experimental observations of
corrugated plasmonic NPs Rodríguez-Fernández  _et al._ (2009). The energy
shift is most evident for the broad dipolar mode at around $\sim 2$ eV, and is
foremost the result of the area increase of the rough surface as compared to
the smooth NP. Evaluation of this area from the mesh parameters yields an
effective radius of $R_{\mathrm{eff}}=76.2$ nm for $\mathrm{rms}=2$ nm, and
$R_{\mathrm{eff}}=79.6$ nm for $\mathrm{rms}=4$ nm. Indeed the EEL spectra of
smooth spheres of said effective radii reproduce accurately the position of
the dipolar mode (see Fig. S.9 in the SI).

In addition to the redshift, the spectra of the corrugated NPs feature a large
number of low-intensity peaks. These new spectral features are the result of
two factors. Due to the breaking of the spherical symmetry, the prior
degenerate modes associated with the same angular momentum $\ell$ but
different $m$ number, now exhibit a small energy difference. As a result of
the lift of the degeneracy, the sharp peaks observed in the spectrum of the
smooth NP are suppressed and, depending on the size of this energy difference
with respect to the linewidth of the degenerate mode, they are either split or
broadened. Moreover, additional spectral features may arise from hot-spots,
namely protuberances of large curvature that strongly enhance and confine the
incident field Urbieta  _et al._ (2018). It is important to note here that, as
the features become increasingly smaller, a rigorous description of the system
requires the implementation of nonlocal effects in the method Stamatopoulou
and Tserkezis (2022); Wegner  _et al._ (2023), which effectively introduces a
cutoff in the scattering from large wavevector components Wiener  _et al._
(2012); Mortensen (2021).

Figure 6: Average EEL probability in the interaction between a spherical NP
featuring surface roughness and an electron beam passing with velocity
$v=0.7c$ (kinetic energy $\approx 200$ keV) at distance $b=100$ nm. The solid
lines correspond to NPs, whose shapes deviate from that of a perfectly smooth
sphere of radius $R=75$ nm (gray shaded area) by root-mean square roughness
values $\mathrm{rms}=2$ nm (dark red curve) and $\mathrm{rms}=4$ nm (dark blue
curve). The average EEL probability corresponds to the average values obtained
with DGTD for $6$ different rough meshes characterized by the same
$\mathrm{rms}$.

##  IV Conclusion

We have presented an analytic and a numerical method for the study of
spherical structures excited by fast electron beams. Based on Mie theory, we
have derived formulas for the calculation of the EEL and CL probability that
are valid for both aloof and penetrating electron beams, as is typically the
practice in EEL and CL measurements. Focusing on the plasmon oscillations of a
metallic NP as a testbed, we compared the analytic theory with numerical
simulations performed using the DGTD method, and found excellent agreement. We
discussed the applicability and limitations of each method, particularly for
grazing trajectories and at energies near the surface- and bulk-plasmon
resonances, and showcased the flexibility of the DGTD method by studying a NP
with different degrees of surface corrugation, which can lead to resonance
shifts and splittings due to the lifting of mode degeneracy. We thus believe
that both methods are essential and complementary for exploring collective
optical excitations in matter and for interpreting experimental observations.

## Acknowledgments

We thank F. Intravaia, C. Maciel-Escudero, and B. Beverungen for stimulating
discussions. K. B. acknowledges funding by the German Research Foundation
(DFG) in the framework of the Collaborative Research Center 1375 “Nonlinear
Optics down to Atomic Scales (NOA)” (Project No. 398816777). N. A. M. is a
VILLUM Investigator supported by VILLUM Fonden (grant No. 16498). The Center
for Polariton-driven Light-Matter Interactions (POLIMA) is funded by the
Danish National Research Foundation (Project No. DNRF165).

## Author contribution

P. E. S and W. Z. contributed equally to this work. P. E. S performed the
analytic study and W. Z. the DGTD simulations. All authors participated in the
discussion of the results and the writing of the manuscript.

## Appendix

The field expansion coefficients entering Eqs. (4) and (7) are given by

$\displaystyle b_{\ell m}^{\textrm{I}}$
$\displaystyle=T_{M_{\ell}}^{11}b_{\ell
m}^{0,\textrm{I}}+T_{M_{\ell}}^{21}b_{\ell m}^{0,\textrm{II}}-b_{\ell
m}^{0,\textrm{II}},$ (11a) $\displaystyle b_{\ell m}^{\textrm{II}}$
$\displaystyle=T_{M_{\ell}}^{12}b_{\ell
m}^{0,\textrm{I}}+T_{M_{\ell}}^{22}b_{\ell m}^{0,\textrm{II}}-b_{\ell
m}^{0,\textrm{I}}|_{\mathrm{air}},$ (11b) $\displaystyle a_{\ell
m}^{\textrm{I}}$ $\displaystyle=T_{E_{\ell}}^{11}a_{\ell
m}^{0,\textrm{I}}+T_{E_{\ell}}^{21}a_{\ell m}^{0,\textrm{II}}-a_{\ell
m}^{0,\textrm{II}},$ (11c) $\displaystyle a_{\ell m}^{\textrm{II}}$
$\displaystyle=T_{E_{\ell}}^{12}a_{\ell
m}^{0,\textrm{I}}+T_{E_{\ell}}^{22}a_{\ell m}^{0,\textrm{II}}-a_{\ell
m}^{0,\textrm{I}}|_{\mathrm{air}},$ (11d)

where $a/b_{\ell m}^{0,\mathrm{I}}$ and $a/b_{\ell m}^{0,\mathrm{II}}$
correspond to the fields generated by the part of the electron trajectory
lying in region I (inside the NP) and region II (outside the NP),
respectively, as illustrated in Fig. 1. The coefficients are found by the
following expressions

$b_{\ell
m}^{0,\mathrm{II}}=-\frac{ik_{0}^{2}e}{\varepsilon_{0}\omega}\frac{m}{\sqrt{\ell(\ell+1)}}\left[\mathcal{M}_{\ell
m}K_{m}\left(\frac{\omega
b}{v\gamma_{0}}\right)-ik_{0}\int_{-z_{\mathrm{e}}}^{z_{\mathrm{e}}}dz\,\exp{(i\omega
z/v)}h_{\ell}^{+}(k_{0}r){Y^{m}_{\ell}}\left(\theta,0\right)\right],$ (12a)
$a_{\ell
m}^{0,\mathrm{II}}=\frac{ik_{0}^{2}e}{\varepsilon_{0}\omega}\frac{1}{\sqrt{\ell(\ell+1)}}\bigg{[}\frac{\mathcal{N}_{\ell
m}}{\beta\gamma_{0}}K_{m}\left(\frac{\omega
b}{v\gamma_{0}}\right)-\frac{i}{b}\int_{-z_{\mathrm{e}}}^{z_{\mathrm{e}}}dz\,\exp{(i\omega
z/v)}\Big{\\{}\mathcal{H}_{\ell m}^{+}(k_{0}z)+\mathcal{H}_{\ell
m}^{-}(k_{0}z)\Big{\\}}\bigg{]},$ (12b) $b_{\ell
m}^{0,\mathrm{I}}=-\frac{ik_{0}^{2}e}{\varepsilon_{0}\omega}\frac{m}{\sqrt{\ell(\ell+1)}}ik\int_{-z_{\mathrm{e}}}^{z_{\mathrm{e}}}dz\,\exp{(i\omega
z/v)}j_{\ell}(kr){Y_{\ell}^{m}}\left(\theta,0\right),$ (12c) and $a_{\ell
m}^{0,\mathrm{I}}=\frac{ik_{0}^{2}e}{\varepsilon_{0}\omega}\frac{1}{\sqrt{\ell(\ell+1)}}\frac{i}{b}\int_{-z_{\mathrm{e}}}^{z_{\mathrm{e}}}dz\,\exp{(i\omega
z/v)}\Big{\\{}\mathcal{J}_{\ell m}^{-}(kz)+\mathcal{J}_{\ell
m}^{+}(kz)\Big{\\}}.$ (12d)

In Eqs. (11) the notation $a/b_{\ell m}^{0,\textrm{I}}|_{\mathrm{air}}$
indicates evaluation of the terms in Eqs. (12c) and (12d) in air
($\varepsilon=1$, $k=k_{0}$). Moreover, in Eqs. (7) and (12) we have set

$\mathcal{M}_{\ell
m}=i^{\ell+m}\sqrt{\frac{2\ell+1}{\pi}\frac{(\ell-m)!}{(\ell+m)!}}\frac{(2m-1)!!}{(\beta\gamma_{0})^{m}}G_{\ell-m}^{m+1/2}\left(\frac{1}{\beta}\right),$
(13)

and

$\mathcal{N}_{\ell
m}=c_{\ell}^{m}\mathcal{M}_{\ell\,m+1}-c_{\ell}^{-m}\mathcal{M}_{\ell\,m-1},$
(14)

with

$c_{\ell}^{m}=\frac{1}{2}\sqrt{(\ell-m)(\ell+m+1)},$ (15)

where $G_{\ell-m}^{m+1/2}(x)$ is the Gegenbauer polynomial. Eq. (13) holds for
$m\geq 0$, while $\mathcal{M}_{\ell-m}=(-1)^{m}\mathcal{M}_{\ell m}$.
Additionally, in Eqs. (7) and (12), we have set

$\displaystyle\mathcal{F}_{\ell m}^{\pm}(k_{n}z)=\mp c_{\ell}^{\pm
m}\bigg{\\{}\frac{k_{n}b^{2}}{r}f^{\prime}_{\ell}(k_{n}r)Y_{\ell}^{m\pm
1}(\theta,0)$
$\displaystyle\pm\frac{zb}{r^{2}}f_{\ell}(k_{n}r)\left[c_{\ell}^{\pm
m+1}Y_{\ell}^{m\pm 2}(\theta,0)-c_{\ell}^{\pm m}Y_{\ell}^{m}(\theta,0)\right]$
$\displaystyle+(1\pm m)f_{\ell}(k_{n}r)Y_{\ell}^{m\pm 1}(\theta,0)\bigg{\\}}.$
(16)

Eq. (Appendix) holds for any type of spherical Bessel function
$f_{\ell}(k_{n}r)$ evaluated in any medium $n$ (the prime here and on any
other Bessel function denotes the derivative of the function with respect to
the argument). In particular, in Eq. (12d) we use expression (Appendix) for
$\mathcal{F}_{\ell m}^{\pm}(k_{n}z)=\mathcal{J}_{\ell m}^{\pm}(kz)$ and
$f_{\ell}=j_{\ell}$, whereas in Eq. (12b) we use $\mathcal{F}_{\ell
m}^{\pm}(k_{n}z)=\mathcal{H}_{\ell m}^{\pm}(k_{0}z)$ and
$f_{\ell}=h^{+}_{\ell}$.

Finally in Eqs. (11) we have introduced the Mie coefficients

$\displaystyle T_{E_{\ell}}^{22}$ $\displaystyle=\frac{\varepsilon
j_{\ell}(kR)\Psi_{\ell}^{\prime}(k_{0}R)-\Psi_{\ell}^{\prime}(kR)j_{\ell}(k_{0}R)}{h^{+}_{\ell}(k_{0}R)\Psi_{\ell}^{\prime}(kR)-\varepsilon\xi_{\ell}^{\prime}(k_{0}R)j_{\ell}(kR)},$
(17a) $\displaystyle T_{M_{\ell}}^{22}$
$\displaystyle=\frac{j_{\ell}(kR)\Psi_{\ell}^{\prime}(k_{0}R)-\Psi_{\ell}^{\prime}(kR)j_{\ell}(k_{0}R)}{h^{+}_{\ell}(k_{0}R)\Psi_{\ell}^{\prime}(kR)-\xi_{\ell}^{\prime}(k_{0}R)j_{\ell}(kR)},$
(17b) $\displaystyle T_{E_{\ell}}^{21}$
$\displaystyle=-\frac{i\sqrt{\varepsilon}/(k_{0}R)}{h_{\ell}^{+}(k_{0}R)\Psi^{\prime}_{\ell}(kR)-\varepsilon\xi^{\prime}_{\ell}(k_{0}R)j_{\ell}(kR)},$
(17c) $\displaystyle T_{M_{\ell}}^{21}$
$\displaystyle=-\frac{i/(k_{0}R)}{h_{\ell}^{+}(k_{0}R)\Psi^{\prime}_{\ell}(kR)-\xi^{\prime}_{\ell}(k_{0}R)j_{\ell}(kR)},$
(17d) $\displaystyle T_{E_{\ell}}^{11}$
$\displaystyle=\frac{\varepsilon\xi^{\prime}_{\ell}(k_{0}R)h_{\ell}^{+}(kR)-h_{\ell}^{+}(k_{0}R)\xi^{\prime}_{\ell}(kR)}{h_{\ell}^{+}(k_{0}R)\Psi^{\prime}_{\ell}(kR)-\varepsilon\xi^{\prime}_{\ell}(k_{0}R)j_{\ell}(kR)},$
(17e) $\displaystyle T_{M_{\ell}}^{11}$
$\displaystyle=\frac{\xi^{\prime}_{\ell}(k_{0}R)h_{\ell}^{+}(kR)-h_{\ell}^{+}(k_{0}R)\xi^{\prime}_{\ell}(kR)}{h_{\ell}^{+}(k_{0}R)\Psi^{\prime}_{\ell}(kR)-\xi^{\prime}_{\ell}(k_{0}R)j_{\ell}(kR)},$
(17f) $\displaystyle T_{E_{\ell}}^{12}$
$\displaystyle=-\frac{i/(k_{0}R)}{h_{\ell}^{+}(k_{0}R)\Psi^{\prime}_{\ell}(kR)-\varepsilon\xi^{\prime}_{\ell}(k_{0}R)j_{\ell}(kR)},$
(17g) $\displaystyle T_{M_{\ell}}^{12}$
$\displaystyle=-\frac{i/(\sqrt{\varepsilon}k_{0}R)}{h_{\ell}^{+}(k_{0}R)\Psi^{\prime}_{\ell}(kR)-\xi^{\prime}_{\ell}(k_{0}R)j_{\ell}(kR)},$
(17h)

where we have adopted the notation of the Riccati–Bessel functions
$\Psi_{\ell}(x)=xj_{\ell}(x)$ and $\xi_{\ell}(x)=xh^{+}_{\ell}(x)$.

## References

  * García de Abajo (2010) F. J. García de Abajo, Rev. Mod. Phys. 82, 209 (2010).
  * Watanabe (1956) H. Watanabe, J. Phys. Soc. Japan 11, 112 (1956).
  * Ritchie (1957) R. H. Ritchie, Phys. Rev. 106, 874 (1957).
  * Yamamoto  _et al._ (2001) N. Yamamoto, K. Araya, and F. J. García de Abajo, Phys. Rev. B 64, 205419 (2001).
  * Vesseur  _et al._ (2007) E. J. R. Vesseur, R. de Waele, M. Kuttge, and A. Polman, Nano Lett. 7, 2843 (2007).
  * Coenen  _et al._ (2013) T. Coenen, J. van de Groep, and A. Polman, ACS Nano 7, 1689 (2013).
  * Matsukata  _et al._ (2019) T. Matsukata, N. Matthaiakakis, T.-a. Yano, M. Hada, T. Tanaka, N. Yamamoto, and T. Sannomiya, ACS Photonics 6, 2320 (2019).
  * Fiedler  _et al._ (2022) S. Fiedler, P. E. Stamatopoulou, A. Assadillayev, C. Wolff, H. Sugimoto, M. Fujii, N. A. Mortensen, S. Raza, and C. Tserkezis, Nano Lett. 22, 2320 (2022).
  * Lagos  _et al._ (2017) M. J. Lagos, A. Trügler, U. Hohenester, and P. E. Batson, Nature 543, 529 (2017).
  * Hage  _et al._ (2019) F. S. Hage, D. M. Kepaptsoglou, Q. M. Ramasse, and L. J. Allen, Phys. Rev. Lett. 122, 016103 (2019).
  * Maciel-Escudero  _et al._ (2020) C. Maciel-Escudero, A. Konečná, R. Hillenbrand, and J. Aizpurua, Phys. Rev. B 102, 115431 (2020).
  * Polman  _et al._ (2019) A. Polman, M. Kociak, and F. J. García de Abajo, Nat. Mater. 18, 1158 (2019).
  * García de Abajo and Di Giulio (2021) F. J. García de Abajo and V. Di Giulio, ACS Photonics 8, 945 (2021).
  * Egerton (2011) R. F. Egerton,  _Electron energy-loss spectroscopy in the electron microscope_ (Springer Science & Business Media, New York, 2011).
  * Wu  _et al._ (2018) Y. Wu, G. Li, and J. P. Camden, Chem. Rev. 118, 2994 (2018).
  * Batson and Silcox (1983) P. Batson and J. Silcox, Phys. Rev. B 27, 5224 (1983).
  * Schmidt  _et al._ (2012) F.-P. Schmidt, H. Ditlbacher, U. Hohenester, A. Hohenau, F. Hofer, and J. R. Krenn, Nano Lett. 12, 5780 (2012).
  * Koh  _et al._ (2009) A. L. Koh, K. Bao, I. Khan, W. E. Smith, G. Kothleitner, P. Nordlander, S. A. Maier, and D. W. McComb, ACS Nano 3, 3015 (2009).
  * Barwick  _et al._ (2008) B. Barwick, H. S. Park, O.-H. Kwon, J. S. Baskin, and A. H. Zewail, Science 322, 1227 (2008).
  * Barwick  _et al._ (2009) B. Barwick, D. J. Flannigan, and A. H. Zewail, Nature 462, 902 (2009).
  * Ferrell and Echenique (1985) T. L. Ferrell and P. M. Echenique, Phys. Rev. Lett. 55, 1526 (1985).
  * Garcia-Molina  _et al._ (1985) R. Garcia-Molina, A. Gras-Marti, and R. H. Ritchie, Phys. Rev. B 31, 121 (1985).
  * Rivacoba  _et al._ (1992) A. Rivacoba, N. Zabala, and P. M. Echenique, Phys. Rev. Lett. 69, 3362 (1992).
  * Zabala  _et al._ (1989) N. Zabala, A. Rivacoba, and P. M. Echenique, Surf. Sci. 209, 465 (1989).
  * Echenique  _et al._ (1987) P. M. Echenique, J. Bausells, and A. Rivacoba, Phys. Rev. B 35, 1521 (1987).
  * Ritchie and Howie (1988) R. H. Ritchie and A. Howie, Phil. Mag. A 58, 753 (1988).
  * Kröger (1970) E. Kröger, Z. Phys. 235, 403 (1970).
  * García de Abajo (1999) F. J. García de Abajo, Phys. Rev. B 59, 3095 (1999).
  * Hage  _et al._ (2020) F. S. Hage, G. Radtke, D. M. Kepaptsoglou, M. Lazzeri, and Q. M. Ramasse, Science 367, 1124 (2020).
  * Yankovich  _et al._ (2019) A. B. Yankovich, B. Munkhbat, D. G. Baranov, J. Cuadra, E. Olsén, H. Lourenço-Martins, L. H. G. Tizei, M. Kociak, E. Olsson, and T. Shegai, Nano Lett. 19, 8171 (2019).
  * Pogorzelski and Yeh (1973) R. Pogorzelski and C. Yeh, Phys. Rev. A 8, 137 (1973).
  * Ginzburg (1996) V. L. Ginzburg, Phys. Usp. 39, 973 (1996).
  * Yamamoto  _et al._ (1996) N. Yamamoto, H. Sugiyama, and A. Toda, Proc. R. Soc. London. Ser. A 452, 2279 (1996).
  * Gallinet  _et al._ (2015) B. Gallinet, J. Butet, and O. J. F. Martin, Laser Photonics Rev. 9, 577 (2015).
  * García de Abajo and Howie (2002) F. J. García de Abajo and A. Howie, Phys. Rev. B 65, 115418 (2002).
  * Hohenester (2014) U. Hohenester, Comput. Phys. Commun. 185, 1177 (2014).
  * Pomplun  _et al._ (2007) J. Pomplun, S. Burger, L. Zschiedrich, and F. Schmidt, Phys. Status Solidi B 244, 3419 (2007).
  * Burger  _et al._ (2012) S. Burger, J. Pomplun, and F. Schmidt, “Finite element methods for computational nano-optics,” in  _Encyclopedia of Nanotechnology_, edited by B. Bhushan (Springer Netherlands, Dordrecht, 2012) pp. 837–843.
  * Das  _et al._ (2012) P. Das, T. K. Chini, and J. Pond, J. Phys. Chem. C 116, 15610 (2012).
  * Cao  _et al._ (2015) Y. Cao, A. Manjavacas, N. Large, and P. Nordlander, ACS Photonics 2, 369 (2015).
  * Hesthaven and Warburton (2007) J. S. Hesthaven and T. Warburton,  _Nodal discontinuous Galerkin methods: algorithms, analysis, and applications_ (Springer Science & Business Media, New York, 2007).
  * Matyssek  _et al._ (2011) C. Matyssek, J. Niegemann, W. Hergert, and K. Busch, Photon. Nanostruct. Fundam. Appl. 9, 367 (2011).
  * Busch  _et al._ (2011) K. Busch, M. König, and J. Niegemann, Laser Photonics Rev. 5, 773 (2011).
  * Husnik  _et al._ (2013) M. Husnik, F. von Cube, S. Irsen, S. Linden, J. Niegemann, K. Busch, and M. Wegener, Nanophotonics 2, 241 (2013).
  * Prämassing  _et al._ (2021) M. Prämassing, T. Kiel, S. Irsen, K. Busch, and S. Linden, Phys. Rev. B 103, 115403 (2021).
  * Maradudin and Mills (1975) A. A. Maradudin and D. L. Mills, Phys. Rev. B 11, 1392 (1975).
  * Maradudin and Méndez (2007) A. A. Maradudin and E. R. Méndez, Sci. Prog. 90, 161 (2007).
  * Mie (1908) G. Mie, Ann. Phys. 330, 377 (1908).
  * Jackson (1998) J. D. Jackson,  _Classical Electrodynamics_, 3rd ed. (John Wiley & Sons, New York, 1998).
  * Lucas and Šunjić (1972) A. A. Lucas and M. Šunjić, Prog. Surf. Sci. 2, 75 (1972).
  * Bohren and Huffman (1983) C. F. Bohren and D. R. Huffman,  _Absorption and Scattering of Light by Small Particles_ (John Wiley & Sons, New York, 1983).
  * Batson  _et al._ (2011) P. E. Batson, A. Reyes-Coronado, R. G. Barrera, A. Rivacoba, P. M. Echenique, and J. Aizpurua, Nano Lett. 11, 3388 (2011).
  * Christensen  _et al._ (2014) T. Christensen, W. Yan, S. Raza, A.-P. Jauho, N. A. Mortensen, and M. Wubs, ACS Nano 8, 1745 (2014).
  * Zouros  _et al._ (2020) G. P. Zouros, G. D. Kolezas, N. A. Mortensen, and C. Tserkezis, Phys. Rev. B 101, 085416 (2020).
  * Gonçalves and García de Abajo (2023) P. A. D. Gonçalves and F. J. García de Abajo, Nano Lett. 23, 4242 (2023).
  * Moeferdt (2017) M. Moeferdt,  _Nonlocal and nonlinear properties of plasmonic nanostructures within the hydrodynamic Drude model_ , Ph.D. thesis, Humboldt Universität zu Berlin (Germany) (2017).
  * Link and El-Sayed (1999) S. Link and M. A. El-Sayed, J. Phys. Chem. B 103, 4212 (1999).
  * Grzelczak  _et al._ (2010) M. Grzelczak, J. Vermant, E. M. Furst, and L. M. Liz-Marzán, ACS Nano 4, 3591 (2010).
  * Nappa  _et al._ (2005) J. Nappa, G. Revillod, I. Russier-Antoine, E. Benichou, C. Jonin, and P. Brevet, Phys. Rev. B 71, 165407 (2005).
  * Salomon  _et al._ (2013) A. Salomon, M. Zielinski, R. Kolkowski, J. Zyss, and Y. Prior, J. Phys. Chem. C 117, 22377 (2013).
  * Zheng  _et al._ (2015) P. Zheng, S. K. Cushing, S. Suri, and N. Wu, Phys. Chem. Chem. Phys. 17, 21211 (2015).
  * Yoon  _et al._ (2019) J. H. Yoon, F. Selbach, L. Schumacher, J. Jose, and S. Schlücker, ACS Photonics 6, 642 (2019).
  * Benz  _et al._ (2016) F. Benz, M. K. Schmidt, A. Dreismann, R. Chikkaraddy, Y. Zhang, A. Demetriadou, C. Carnegie, H. Ohadi, B. de Nijs, R. Esteban, J. Aizpurua, and J. J. Baumberg, Science 354, 726 (2016).
  * Urbieta  _et al._ (2018) M. Urbieta, M. Barbry, Y. Zhang, P. Koval, D. Sánchez-Portal, N. Zabala, and J. Aizpurua, ACS Nano 12, 585 (2018).
  * Rodríguez-Fernández  _et al._ (2009) J. Rodríguez-Fernández, A. M. Funston, J. Pérez-Juste, R. A. Álvarez-Puebla, L. M. Liz-Marzán, and P. Mulvaney, Phys. Chem. Chem. Phys. 11, 5909 (2009).
  * Loth  _et al._ (2023) F. Loth, T. Kiel, K. Busch, and P. T. Kristensen, J. Opt. Soc. Am. B 40, B1 (2023).
  * Stamatopoulou and Tserkezis (2022) P. E. Stamatopoulou and C. Tserkezis, Opt. Mater. Express 12, 1869 (2022).
  * Wegner  _et al._ (2023) G. Wegner, D.-N. Huynh, N. A. Mortensen, F. Intravaia, and K. Busch, Phys. Rev. B 107, 115425 (2023).
  * Wiener  _et al._ (2012) A. Wiener, A. F. Fernández-Domínguez, A. P. Horsfield, J. B. Pendry, and S. A. Maier, Nano Lett. 12, 3308 (2012).
  * Mortensen (2021) N. A. Mortensen, Nanophotonics 10, 2563 (2021).


# Cosmological chirality and magnetic fields from parity violating particle
decays

Tanmay Vachaspati∗, Alexander Vilenkin† ∗Physics Department, Arizona State
University, Tempe, AZ 85287, USA.
†Institute of Cosmology, Department of Physics and Astronomy, Tufts
University, Medford, MA 02155, USA.

###### Abstract

We estimate the chirality of the cosmological medium due to parity violating
decays of standard model particles, focusing on the example of tau leptons.
The non-trivial chirality is however too small to make a significant
contribution to the cosmological magnetic field via the chiral-magnetic
effect.

##  I Introduction

The last few decades have seen growing interest in cosmic magnetic fields on
several fronts [1]. Several ideas have been proposed that can generate
magnetic fields in cosmology, some of which are directly tied to known
particle physics [2, 3, 4] and its possible extensions [5, 6, 7, 8, 1]. The
magneto-hydrodynamical (MHD) evolution of cosmological magnetic fields is now
understood quite well on the basis of analytical arguments [9, 10] and direct
simulations [11]. There are claims for an indirect lower bound on the
cosmological magnetic field strength [12, 13, 14, 15, 16], though not without
debate [17, 18], and more direct evidence [19]. Concurrently there are claims
of a parity violating signature that can be used to measure the magnetic field
helicity spectrum [20, 21] though with no significant detections as yet [22,
23].

In parallel to these developments, motivated by heavy-ion collision
experiments [24], there has been renewed interest in chiral effects in
plasmas, namely the chiral-magnetic [25] and chiral-vortical [26] effects (CME
and CVE respectively). The CME and CVE have also been applied to the evolution
of cosmological and astrophysical magnetic fields [5, 27, 28, 29, 30, 31, 32,
33]. In this paper we discuss how CME and CVE can effectively arise in
standard cosmology with standard particle interactions due to the parity-
violating decays of heavy leptons and quarks. The basic idea is that the
standard model has a number of unstable particles that decay at various
cosmological epochs, primarily due to the weak interactions. Since the weak
interactions violate parity, the decay products are chiral and this provides a
net particle helicity to the cosmological medium. The net particle helicity in
principle leads to electric currents via the CME that can generate magnetic
helicity. However, accounting only for decays of standard model particles, the
net particle helicity is too small to significantly affect cosmological
magnetic fields and their helicity.

We start by describing the physical effect in some detail in the context of
the tau lepton in Sec. II, where we also estimate the induced electric
currents. We find an upper bound to the magnetic helicity that can be
generated due to chiral effects in Sec. III. Our conclusions are summarized
and discussed in Sec. IV.

##  II Chirality production in tau decays

To illustrate the physics of the effect, in this section we will discuss the
decay of tau leptons in the background of a magnetic field and fluid
vorticity. Except for small differences, the physics carries over to the case
of decays of other particles.

###  II.1 Particle decay

Tau leptons decay into electrons (or muons) and neutrinos

$\tau^{-}\to e^{-}+\nu_{\tau}+{\bar{\nu}}_{e}$ (1)

and anti-tau into positrons and neutrinos

$\tau^{+}\to e^{+}+{\bar{\nu}}_{\tau}+\nu_{e}$ (2)

These decays violate parity since they proceed primarily by the weak
interactions. Therefore the tau predominantly decays into a relativistic left-
handed electron, while an anti-tau decays into a relativistic right-handed
positron. Due to the lepton asymmetry of the universe there are more taus than
anti-taus, and the cosmological medium gains net left-handed chirality as
tau’s decay.

The decay product electrons are chiral since they are produced by the weak
interactions, but chirality is not preserved for massive particles. Instead,
as emphasized in Ref. [34] in the context of supernovae and neutron stars,
chirality is nearly equal to helicity for ultrarelativistic particles, so it
is better to think of the final electrons as being in a definite helicity
state. Helicity can only change due to particle interactions. We shall adopt
this view in what follows.

The $\tau$ mass is $m_{\tau}=1777~{}{\rm MeV}$ and the $\tau$ lifetime in its
rest frame is $\tau_{\tau}=2.9\times 10^{-13}~{}{\rm s}$. However, the
decaying taus are constantly reproduced by reactions inverse to (1),
(2),111Tau-particles are also produced and destroyed in scattering reactions
like $\tau^{-}+{\nu}_{e}\to e^{-}+\nu_{\tau}$. We disregard them in what
follows, since they do not change the order of magnitude of the effect. so the
number density of taus, $n_{\tau}$, remains comparable to that of photons
until the time

$t_{\tau}\sim 10^{-7}~{}{\rm s},$ (3)

when the cosmic temperature drops to $T\sim m_{\tau}$. At later times
$n_{\tau}$ decreases exponentially.

The particle helicity density, $n_{\chi}$, is produced in tau decays and is
dissipated by helicity flipping scatterings and due to the chiral anomaly. The
latter is proportional to $\alpha^{3}B^{2}$ [35], where $\alpha\approx 1/137$
is the fine structure constant and $B$ the magnetic field strength, and is
much slower than helicity flipping scatterings for vanishing or weak magnetic
fields. We will ignore the anomalous flipping for now but will discuss it in
Sec. LABEL:Bgen when we consider the effect of particle chirality on the
generation of magnetic fields. The evolution of particle helicity density can
be described by the kinetic equation in the relaxation time approximation,

$\frac{d}{dt}(a^{3}n_{\chi})=\frac{a^{3}}{\tau_{d}}(\delta n_{\tau}-\delta
n_{\tau}^{\rm eq})-\frac{a^{3}n_{\chi}}{\tau_{\chi}},$ (4)

where

$\delta n_{\tau}=n_{\tau}^{+}-n_{\tau}^{-},$ (5)

$n_{\tau}^{-}$ and $n_{\tau}^{+}$ are the densities of tau and anti-tau
particles, respectively, $\delta n_{\tau}^{\rm eq}$ is the equilibrium value
of $\delta n_{\tau}$, $\tau_{d}\sim(T/m_{\tau})\tau_{\tau}$ is the decay time
of taus (assuming that $T>m_{\tau}$ and with time dilation taken into account)
and $\tau_{\chi}^{-1}$ is the electron helicity flipping rate. At $T\gg
m_{e}$, the helicity flipping rate is suppressed by a factor $m_{e}^{2}/T^{2}$
compared to the scattering rate $\alpha T$ [36] (earlier estimates of the
scattering rate were suppressed by another factor of $\alpha$ [34]),

$\tau_{\chi}\sim\frac{1}{\alpha T}\frac{T^{2}}{m_{e}^{2}}.$ (6)

The excess of anti-tau’s over tau’s, $\delta n_{\tau}$, decreases due to tau
decay and is described by the equation,

$\frac{d}{dt}(a^{3}\delta n_{\tau})=\frac{a^{3}}{\tau_{d}}(\delta
n_{\tau}^{\rm eq}-\delta n_{\tau}).$ (7)

At temperatures below the elecroweak phase transition, $T\lesssim T_{\rm
EW}\sim 100$ GeV, we have $\tau_{d}\ll t$, where $t$ is the cosmic time222This
is easily verified using the relation $t\sim m_{\rm P}/\sqrt{N}T^{2}$, where
$m_{\rm P}$ is the Planck mass and $N$ is the number of particle species in
equilibrium.. This means that the equilibrium density of taus establishes very
quickly (compared to the Hubble time), and the approximate solution of (7) is
$\delta n_{\tau}\approx\delta n_{\tau}^{\rm eq}$. Inserting (7) in (4) and
then using $\delta n_{\tau}\approx\delta n_{\tau}^{\rm eq}$ we have

$\frac{d}{dt}(a^{3}n_{\chi})=-\frac{d}{dt}\left(a^{3}\delta n_{\tau}^{\rm
eq}\right)-\frac{a^{3}n_{\chi}}{\tau_{\chi}}.$ (8)

With a given $\delta n_{\tau}^{\rm eq}$, this equation can be solved in
quadratures, but we shall instead find an approximate solution. Since we are
in the regime where $\tau_{\chi}\ll t$, the term on the left-hand side can be
neglected and we obtain

$n_{\chi}\approx-\tau_{\chi}T^{3}\frac{d}{dt}\left(\frac{\delta n_{\tau}^{\rm
eq}}{T^{3}}\right),$ (9)

where we have used $aT\approx{\rm const}$.

Once we determine the equilibrium excess of anti-taus over taus, denoted by
$\delta n_{\tau}^{\rm eq}$, we can estimate the chirality density of the
universe due to tau decays using (9).

###  II.2 Equilibrium density

The equilibrium density $\delta n_{\tau}^{\rm eq}$ is given by

$\delta n_{\tau}^{\rm
eq}=\frac{1}{2\pi^{2}}\int_{0}^{\infty}dpp^{2}\left[f\left(\frac{E-\mu_{\tau}}{T}\right)-f\left(\frac{E+\mu_{\tau}}{T}\right)\right],$
(10)

where $f(x)=(e^{x}+1)^{-1}$ is the Fermi distribution,
$E=\sqrt{p^{2}+m_{\tau}^{2}}$, and $\mu_{\tau}$ is the chemical potential of
$\tau$ particles. At $T\gg m_{\tau},\mu_{\tau}$ we can expand the integrand in
Eq. (10) in powers of $m_{\tau}^{2}/p^{2}$ and $\mu_{\tau}/T$. The
integrations are then easily performed and we find

$\delta n_{\tau}^{\rm
eq}\approx\frac{\mu_{\tau}T^{2}}{6}\left(1-\frac{3m_{\tau}^{2}}{2\pi^{2}T^{2}}\right).$
(11)

We assume that the baryon and/or lepton asymmetry of the universe was
generated at $T\gg T_{EW}$ by some interactions beyond the Standard Model, for
example by $(B-L)$-violating leptoquark decays. This asymmetry was then
redistributed between the Standard Model leptons and quarks by sphaleron
processes, so at $T\ll T_{EW}$ we expect the chemical potentials of light
baryons and leptons to be of the order $\mu/T\sim\eta_{B}$ [37, 38], where
$\eta_{B}\sim 10^{-9}$ is the observed baryon to photon ratio. In the high-
temperature regime, when $T$ is large compared to all relevant particle
masses, we have $\mu_{\tau}/T\approx{\rm const}$, with a mass correction
${\cal O}(m^{2}/T^{2})$ [39]. Then Eq. (11) becomes

$\frac{\delta n_{\tau}^{\rm eq}}{T^{3}}\approx
C\eta_{B}-K\eta_{B}\frac{m_{\tau}^{2}}{T^{2}},$ (12)

where $C$ and $K$ are ${\cal O}(1)$ numerical constants333This estimate
assumes that taus are the heaviest particles present in equilibrium at
temperature $T$. If a heavier particle is present in equilibrium, it too will
contribute to the mass correction and may change the estimate.. The mass
correction term in (12) can be qualitatively understood as follows. As the
temperature decreases, it becomes energetically favorable to transfer the
conserved $\tau$-lepton number from $\tau$-particles to $\tau$-neutrinos. The
excess $\tau$-lepton number is also decreased as a result [39].

Substituting Eq. (12) in (9) we obtain

$n_{\chi}\approx-3K\eta_{B}\tau_{\chi}m_{\tau}^{2}{\dot{T}}.$ (13)

With ${\dot{T}}=-T/2t$, $t\sim m_{\rm P}/T^{2}$ and $\tau_{\chi}$ from Eq.
(6), this gives (omitting numerical factors)

$n_{\chi}\sim\frac{\eta_{B}m_{\tau}^{2}}{\alpha m_{e}^{2}}\frac{T}{m_{\rm
P}}n_{\gamma},$ (14)

where $n_{\gamma}\sim T^{3}$ is the photon number density.

This estimate was derived assuming $T\gg m_{\tau}$, but it still applies at
$T\sim m_{\tau}$. Reactions (1), (2) remain in equilibrium when $T$ drops well
below $m_{\tau}$. In this regime, $\delta n_{\tau}$ and $n_{\chi}$ decrease
exponentially.

Similar formulae can be written down for the decay of other unstable
particles. The largest helicity is injected by the decay of the heaviest
particle into the lightest particle and at the highest temperature.

##  III Magnetic helicity

As noted in Ref. [32], the maximum magnetic helicity that can be obtained from
particle helicity can be derived from the chiral anomaly equation, which can
be written as a conservation law,

$n_{\chi}+\frac{4\pi}{\alpha}h={\rm constant}.$ (15)

where $h=\langle{\bf A}\cdot{\bf B}\rangle$ is the magnetic helicity. Assuming
that the initial magnetic helicity and the final particle helicity vanish, we
get

$h_{\rm max}=\frac{\alpha}{4\pi}n_{\chi}\sim\frac{\eta_{B}m_{\tau}^{2}}{4\pi
m_{e}^{2}}\frac{T}{m_{\rm P}}n_{\gamma}$ (16)

where we have used (14). We compare $h_{\rm max}$ to magnetic helicity that
could be induced due to baryogenesis [3, 4],

$h_{B}\sim\frac{\eta_{B}}{\alpha}n_{\gamma}\sim 10^{-5}\,{\rm cm}^{-3}\sim
10^{-45}\,{\rm G}^{2}\,{\rm Mpc}$ (17)

where we have used the known cosmic baryon number density and are using
natural units. Then

$h_{\rm max}\sim\frac{\alpha m_{\tau}^{2}}{4\pi m_{e}^{2}}\frac{T}{m_{\rm
P}}h_{B}\sim 10^{-10}h_{B}$ (18)

where we have used $T\sim 100\,{\rm GeV}$ in the numerical estimate. Even if
the decay of top quarks with mass $\sim 175\,{\rm GeV}$ to down quarks with
mass $\sim 1\,{\rm MeV}$ is considered, $h_{\rm max}\sim 10^{-6}h_{B}$.
Comparing to observations, even with the most conservative lower bound of
$10^{-19}\,{\rm G}$ on Mpc scales, we get an estimate for the observed
helicity $\sim 10^{-38}\,{\rm G}^{2}\,{\rm Mpc}$.

##  IV Conclusions

We have shown that the decays of certain standard model particles can lead to
a chiral cosmological medium around the epoch of the electroweak phase
transition. The final result for the chiral asymmetry due to tau-lepton decays
is given in (14). However, the asymmetry is suppressed by the baryon to
entropy ratio ($\eta_{B}\sim 10^{-9}$) and the effect on magnetic field
helicity generation is very weak as we have shown in Sec. III. Nonetheless it
is of interest that the cosmological medium was chiral at the earliest moments
even within the standard model of particle physics.

##  V Acknowledgements

We thank the participants of the Nordita workshop on “Quantum Anomalies and
Chiral Magnetic Phenomena”, especially Axel Brandenburg and Kohei Kamada for
feedback. We also thank Matt Baumgart, Cecilia Lunardini, Igor Shovkovy, and
Tracy Slatyer for discussions. TV’s work is supported by the U.S. Department
of Energy, Office of High Energy Physics, under Award No. DE-SC0019470 at
Arizona State University. AV is supported by the National Science Foundation
Award No. 1820872.

## References

  * [1] Tanmay Vachaspati.  Progress on Cosmological Magnetic Fields.  10 2020.
  * [2] T. Vachaspati.  Magnetic fields from cosmological phase transitions.  Phys.Lett., B265:258–261, 1991.
  * [3] John M. Cornwall.  Speculations on primordial magnetic helicity.  Phys. Rev., D56:6146–6154, 1997.
  * [4] Tanmay Vachaspati.  Estimate of the primordial magnetic field helicity.  Phys. Rev. Lett., 87:251302, 2001.
  * [5] M. Joyce and Mikhail E. Shaposhnikov.  Primordial magnetic fields, right-handed electrons, and the Abelian anomaly.  Phys. Rev. Lett., 79:1193–1196, 1997.
  * [6] Michael McNeil Forbes and Ariel R. Zhitnitsky.  Primordial galactic magnetic fields from domain walls at the QCD phase transition.  Phys. Rev. Lett., 85:5268–5271, 2000.
  * [7] Trevor Stevens, Mikkel B. Johnson, Leonard S. Kisslinger, and Ernest M. Henley.  Non-Abelian Higgs model of magnetic field generation during a cosmological first-order electroweak phase transition.  Phys. Rev., D85:063003, 2012.
  * [8] F. Miniati, G. Gregori, B. Reville, and S. Sarkar.  Axion-driven cosmic magnetogenesis during the qcd crossover.  Physical Review Letters, 121(2), Jul 2018.
  * [9] Robi Banerjee and Karsten Jedamzik.  The Evolution of cosmic magnetic fields: From the very early universe, to recombination, to the present.  Phys. Rev., D70:123003, 2004.
  * [10] Karsten Jedamzik and Gunter Sigl.  The Evolution of the Large-Scale Tail of Primordial Magnetic Fields.  Phys.Rev., D83:103005, 2011.
  * [11] Axel Brandenburg, Tina Kahniashvili, Sayan Mandal, Alberto Roper Pol, Alexander G. Tevzadze, and Tanmay Vachaspati.  Evolution of hydromagnetic turbulence from the electroweak phase transition.  Phys. Rev., D96(12):123528, 2017.
  * [12] A. Neronov and I. Vovk.  Evidence for Strong Extragalactic Magnetic Fields from Fermi Observations of TeV Blazars.  Science, 328:73–, April 2010.
  * [13] Warren Essey, Shin’ichiro Ando, and Alexander Kusenko.  Determination of intergalactic magnetic fields from gamma ray data.  Astropart. Phys., 35:135–139, 2011.
  * [14] K. Dolag, M. Kachelriess, S. Ostapchenko, and R. Tomas.  Lower limit on the strength and filling factor of extragalactic magnetic fields.  Astrophys. J., 727:L4, 2011.
  * [15] Matthew Wood, Jonathan Biteau, Regina Caputo, Mattia Di Mauro, and Manuel Meyer.  Preliminary Results of the $Fermi$ High-Latitude Extended Source Catalog.  PoS, ICRC2017:647, 2018.
  * [16] M. Ackermann et al.  The Search for Spatial Extension in High-latitude Sources Detected by the $Fermi$ Large Area Telescope.  Astrophys. J. Suppl., 237(2):32, 2018.
  * [17] Avery E. Broderick, Philip Chang, and Christoph Pfrommer.  The Cosmological Impact of Luminous TeV Blazars I: Implications of Plasma Instabilities for the Intergalactic Magnetic Field and Extragalactic Gamma-Ray Background.  Astrophys.J., 752:22, 2012.
  * [18] Timothy C. Arlen, Vladimir V. Vassiliev, Thomas Weisgarber, Scott P. Wakely, and S. Yusef Shafi.  Intergalactic Magnetic Fields and Gamma Ray Observations of Extreme TeV Blazars.  Astrophys. J., 796(1):18, 2014.
  * [19] Wenlei Chen, James H. Buckley, and Francesc Ferrer.  Search for GeV ? -Ray Pair Halos Around Low Redshift Blazars.  Phys. Rev. Lett., 115:211103, 2015.
  * [20] Hiroyuki Tashiro, Wenlei Chen, Francesc Ferrer, and Tanmay Vachaspati.  Search for CP Violating Signature of Intergalactic Magnetic Helicity in the Gamma Ray Sky.  Mon. Not. Roy. Astron. Soc., 445(1):L41–L45, 2014.
  * [21] Wenlei Chen, Borun D. Chowdhury, Francesc Ferrer, Hiroyuki Tashiro, and Tanmay Vachaspati.  Intergalactic magnetic field spectra from diffuse gamma rays.  Mon. Not. Roy. Astron. Soc., 450(4):3371–3380, 2015.
  * [22] Julia Asplund, Guölaugur Jóhannesson, and Axel Brandenburg.  On the measurement of handedness in fermi large area telescope data, 2020\.
  * [23] M. Kachelriess and B. C. Martinez.  Searching for primordial helical magnetic fields, 2020.
  * [24] Dmitri E. Kharzeev.  The Chiral MagnetoHydroDynamics of QCD fluid at RHIC and LHC.  J. Phys., G38:124061, 2011.
  * [25] A. Vilenkin.  Equilibrium Parity Violating Current in a Magnetic Field.  Phys. Rev., D22:3080–3084, 1980.
  * [26] A. Vilenkin.  Macroscopic Parity Violating Effects: Neutrino Fluxes from Rotating Black Holes and in Rotating Thermal Radiation.  phys. rev., D20:1807–1812, 1979.
  * [27] Alexey Boyarsky, Jurg Frohlich, and Oleg Ruchayskiy.  Self-consistent evolution of magnetic fields and chiral asymmetry in the early Universe.  Phys. Rev. Lett., 108:031301, 2012.
  * [28] Hiroyuki Tashiro, Tanmay Vachaspati, and Alexander Vilenkin.  Chiral Effects and Cosmic Magnetic Fields.  Phys. Rev., D86:105033, 2012.
  * [29] Maxim Dvornikov and Victor B. Semikoz.  Lepton asymmetry growth in the symmetric phase of an electroweak plasma with hypermagnetic fields versus its washing out by sphalerons.  Phys. Rev. D, 87(2):025023, 2013.
  * [30] Maxim Dvornikov and Victor B. Semikoz.  Instability of magnetic fields in electroweak plasma driven by neutrino asymmetries.  JCAP, 05:002, 2014.
  * [31] Maxim Dvornikov and Victor B. Semikoz.  Influence of the turbulent motion on the chiral magnetic effect in the early Universe.  Phys. Rev. D, 95(4):043538, 2017.
  * [32] Axel Brandenburg, Jennifer Schober, Igor Rogachevskii, Tina Kahniashvili, Alexey Boyarsky, Jurg Frohlich, Oleg Ruchayskiy, and Nathan Kleeorin.  The turbulent chiral-magnetic cascade in the early universe.  Astrophys. J., 845(2):L21, 2017.
  * [33] Youhei Masada, Kei Kotake, Tomoya Takiwaki, and Naoki Yamamoto.  Chiral magnetohydrodynamic turbulence in core-collapse supernovae.  Phys. Rev., D98(8):083018, 2018.
  * [34] Dorota Grabowska, David B. Kaplan, and Sanjay Reddy.  Role of the electron mass in damping chiral plasma instability in Supernovae and neutron stars.  Phys. Rev., D91(8):085035, 2015.
  * [35] Daniel G. Figueroa and Mikhail Shaposhnikov.  Anomalous non-conservation of fermion/chiral number in Abelian gauge theories at finite temperature.  JHEP, 04:026, 2018.
  * [36] A. Boyarsky, V. Cheianov, O. Ruchayskiy, and O. Sobol.  Evolution of the Primordial Axial Charge across Cosmic Times.  Phys. Rev. Lett., 126:021801, 2021.
  * [37] V. A. Kuzmin, V. A. Rubakov, and M. E. Shaposhnikov.  Anomalous Electroweak Baryon Number Nonconservation and GUT Mechanism for Baryogenesis.  Phys. Lett., B191:171–173, 1987.
  * [38] Jeffrey A. Harvey and Michael S. Turner.  Cosmological baryon and lepton number in the presence of electroweak fermion number violation.  Phys. Rev., D42:3344–3349, 1990.
  * [39] A. I. Bochkarev, S. Yu. Khlebnikov, and M. E. Shaposhnikov.  Sphalerons and Baryogenesis: Electroweak CP Violation at High Temperatures.  Nucl. Phys., B329:493–518, 1990.