# Supplemental material
Observation of interatomic Coulombic decay induced
by double excitation of helium in nanodroplets

B. Bastian  Wilhelm Ostwald Institute for Physical and Theoretical Chemistry,
University of Leipzig, Linnéstraße 2, 04103 Leipzig, Germany  J. D. Asmussen
L. Ben Ltaief  H. B. Pedersen  Department of Physics and Astronomy, Aarhus
University, 8000 Aarhus C, Denmark  K. Sishodia  S. De  Indian Institute of
Technology Madras, Chennai 600036, India  S. R. Krishnan  Department of
Physics and QuCenDiEM-group, Indian Institute of Technology Madras, Chennai
600036, India  C. Medina  Institute of Physics, University of Freiburg, 79104
Freiburg, Germany  N. Pal  R. Richter  Elettra-Sincrotrone Trieste, 34149
Basovizza, Trieste, Italy  N. Sisourat  Sorbonne Université, CNRS, Laboratoire
de Chimie Physique Matière et Rayonnement, UMR 7614, F-75005 Paris, France  M.
Mudrich<EMAIL_ADDRESS>Department of Physics and Astronomy, Aarhus
University, 8000 Aarhus C, Denmark

(September 1, 2024)

Figure 1: Interactions at a double excitation (Fano) resonance in a helium
nanodroplet. An incident photon $h\nu$ can (1) doubly excite a helium atom or
(2) directly ionize. A doubly excited helium atom will quickly autoionize (1a)
or undergo single deexcitation (1b). The autoionization (1a) and direct
ionization (2) pathways are indistinguishable and give rise to the Fano
resonance profile when varying the photon energy. The resulting photoelectron
may escape the droplet with $35.8\text{\,}\mathrm{e}\mathrm{V}$ kinetic energy
or excite a neighbor atom by inelastic scattering (dashed) and escape as a
slow $14.6\text{\,}\mathrm{e}\mathrm{V}$ electron. In the single deexcitation
pathway (1b), instead the doubly excited atom decays to a singly excited one
while a neigboring atom is ionized by energy transfer and a
$14.6\text{\,}\mathrm{e}\mathrm{V}$ ICD electron is ejected. The aim of this
work is to disentangle the inelastic scattering and ICD contributions to the
formation of slow electrons.

##  I Disentanglement of ICD and inelastically scattered electrons

To spectroscopically disentangle the ICD and inelastic scattering
contributions to electron-energy loss (EEL) electrons around
$E_{e}=$15\text{\,}\mathrm{e}\mathrm{V}$$ electron kinetic energy, the ratio
of EEL to total electron counts is divided by its off-resonant
($E_{h\nu}\leq$59.9\text{\,}\mathrm{e}\mathrm{V}$$) average which gives the
relative count ratio

$r(E_{e})=\frac{c_{\text{\tiny
EEL}}(E_{e})}{c_{\text{total}}(E_{e})}\left<\frac{c_{\text{\tiny
EEL}}}{c_{\text{total}}}\right>_{\text{\\!off}}^{-1}$ (1)

where $c$ denotes absolute counts and $\langle\leavevmode\nobreak\
\rangle_{\text{\\!off}}$ the average over the off-resonant photon energies. If
only inelastic scattering contributed to EEL, $c_{\text{\tiny EEL}}$ would be
proportional to $c_{\text{total}}$ because the inelastic scattering cross
sections are nearly constant across the resonance. Consequently, the ratio
would be constant and $r(E_{e})=1$ across the tuning range of the photon
energy studied here. The ICD channel adds additional counts $c_{\text{\tiny
EEL}}(E_{e})$ which results in values $r(E_{e})>1$ on the resonance. The ratio
of ICD to total EEL electrons is then $\left(r-1\right)/r$ which gives the
number of ICD electrons as a function of EEL or total electron counts

$c_{\text{\tiny ICD}}(E_{e})=\frac{r(E_{e})-1}{r(E_{e})}\,c_{\text{\tiny
EEL}}(E_{e})=\left(r(E_{e})-1\right)\,\left<\frac{c_{\text{\tiny
EEL}}}{c_{\text{total}}}\right>_{\text{\\!off}}\,c_{\text{total}}(E_{e})$ (2)

and the probability of an electron to originate from the ICD process

$p_{\text{\tiny ICD}}(E_{e})=\frac{c_{\text{\tiny
ICD}}(E_{e})}{c_{\text{total}}(E_{e})}=\frac{r(E_{e})-1}{r(E_{e})}\,\frac{c_{\text{\tiny
EEL}}(E_{e})}{c_{\text{total}}(E_{e})}=\left(r(E_{e})-1\right)\,\left<\frac{c_{\text{\tiny
EEL}}}{c_{\text{total}}}\right>_{\text{\\!off}}\;.$ (3)

Uncertainties are determined by Gaussian error propagation assuming $\sqrt{c}$
errors for the EEL and total electron counts and using the standard deviation
of the off-resonant $c_{\text{\tiny EEL}}/c_{\text{total}}$ values for the
error of their average. To evaluate the hemispherical electron analyzer (HEA)
data for which the total electrons were not measured, the ICD probability is
factorized into

$p_{\text{\tiny ICD}}(E_{e})=\frac{c_{\text{\tiny
ICD}}(E_{e})}{c_{\text{total}}(E_{e})}=\frac{c_{\text{\tiny
ICD}}(E_{e})}{c_{\text{\tiny EEL}}(E_{e})}\frac{c_{\text{\tiny
EEL}}(E_{e})}{c_{\text{total}}(E_{e})}=\frac{c_{\text{\tiny
ICD}}(E_{e})}{c_{\text{\tiny EEL}}(E_{e})}\,p_{\text{\tiny EEL}}(E_{e})$ (4)

and the second factor ($p_{\text{\tiny EEL}}$) is taken from VMI data with a
similar average droplet size.

##  II Baseline correction for count ratios from the HEA data

As photoelectrons were not detected in the HEA measurement, measured total
electron scans from a previous beamtime with the same droplet source but a
different detector were used [1]. The resulting skewed baseline was corrected
with a constant slope which is justified by comparing to the present VMI data
where the ratios stem from a single measurement, see Fig. 2 and Fig. 2 (d) in
the main text. The derived normalized EEL to total electron ratios and total
ICD counts are presented Fig. 3 (b) and Fig. 4 (b) in the main text,
respectively.

Figure 2: Linear correction of the count ratio of low kinetic energy
($18\text{\,}\mathrm{e}\mathrm{V}23\text{\,}\mathrm{e}\mathrm{V}$) to total
electrons. VMI and HEA data are shown for the He nanodroplet (HND) expansion
conditions specified in the legend. The HEA data are corrected by a linear fit
to the off-resonant photon energy range around the double excitation
resonance. The agreement of the corrected data with the result from velocity
map imaging (VMI) data justifies this approach.

Figure 3: (a) Average HEA high-resolution EEL spectrum with fit residuals
(gray crosses, right axis) and individual Gaussians (see labels in the legend)
of the multi-peak fit that are not shown in Fig. 4 (a) in the main text. Black
circles are averaged over 39 spectra recorded in the range $h\nu=58.4$ – 62.2
eV. Vertical dashed lines show atomic excitations with values from NIST [2].
Shaded areas correspond to He∗He+ states with boundaries computed for ICD
decay of the He∗∗ droplet excitation into the He${}_{2}^{+*}$ states shown in
(c). (b) Photon energy dependence of the normalized peak area to total
electron ratios for different states with colors corresponding to (a). In two
cases (magenta and gray), the sum of two peak area is considered. Those
components assigned to molecular excitations display pronounced maxima in the
range of the Fano resonance. The ratios for the molecular 4–5
${}^{2}\Sigma/^{1}\Pi$ excitations are given by the second $y$-axis on the
right. (c) Estimate of the total ICD counts in the electron impact excitation
(EEL) and electron impact ionization (EII) range of the electron kinetic
energy. Integration ranges of the electron energy loss are specified in the
legend. The 1s3p/1s4p range is excluded because it overlaps with the EII
region. (d) Photon energy dependence of ICD counts at different levels of
molecular excitation in the 21 – 23 eV range. Higher excitation is
preferentially observed at higher photon energies. (e) and (f) Ratios of EEL
(18 – 23 eV e- kinetic energy loss) to photoelectron counts ($-$2 – 2 eV loss)
as a function of photon energy at different HND sizes specified in the legend.
Total photoelectrons at the $n=2$ and $n=3$ Fano resonances (assigned
according to Ref. [1] Fig. 1) are shown as solid red lines (right axes). The
lower photon energies (e) feature the discussed ICD resonance around 60.5 eV
photon energy. Assuming a similar ICD efficiency for the higher photon
energies (f) would give an increase of the order of 2 % which is in the noise
level.

##  III Count ratios from Fano and Lorentzian fits

To evaluate the ICD probability in (3) or (4) and related quantities, we fit
the measured electron counts as a function of photon energy with the Fano
profile [3]

$c(h\nu)=s_{d}+s_{r}\,\frac{\left(q+\frac{h\nu-
E_{r}}{\Gamma}\right)^{2}}{1+\left(\frac{h\nu-E_{r}}{\Gamma}\right)^{2}}$ (5)

with photon energy $h\nu$, asymmetry (Fano) parameter $q$, resonance energy
and width $E_{r}$ and $\Gamma$, the resonant amplitude $s_{r}$ for double
excitation and the direct ionization amplitude $s_{d}$. In this case, $q$ is
negative and the maximum, resonance and minimum points are

$\displaystyle E_{r}+\Gamma/q,$ $\displaystyle E_{r}$ and $\displaystyle
E_{r}-\Gamma\,q$ (6) with the amplitudes $\displaystyle
s_{d}+\left(q^{2}+1\right)s_{r},$ $\displaystyle s_{d}+q^{2}s_{r}$ and
$\displaystyle s_{d}\,.$ (7)

In case of ICD electrons only the resonant pathway is involved and the shape
is found to reduce to a symmetric Lorentzian profile with $q=0$,
$s_{d}+s_{r}\approx 0$ and the maximum amplitude $s_{d}$ is found at the
resonance energy. To compute count ratios (relative probabilities), we
evaluate all amplitudes at the resonance energies obtained from the fit which
gives $s_{d}+q^{2}s_{r}$. The resonance energy of all fits to the total, EEL
and ICD electrons is close to 60.4 eV. All fit results are summarized in Table
III.

Table 1: Fit results from Fano and Lorentzian fits. The photon energy range
from 59.8 eV to 60.2 eV was excluded for the fits because of the sharply
peaked photoline from an atomic He fraction in the droplet beam. The He+
coincidences are most strongly affected so the fits for the total coincidence
electrons are not reliable and impossible at 16 K. Counts ($s_{d}$, $s_{r}$
and $s_{d}+q^{2}s_{r}$) are tabulated in units of $10^{6}$ for all electrons
(all $e^{-}$) and $10^{4}$ for He${}_{n}^{+}$–$\,e^{-}$ coincidences and HEA
data.

Det. | T/K | p/bar | Data | Electrons |  $q$ |  $E_{r}$/eV |  $\Gamma$/eV |  $s_{d}$ |  $s_{r}$ | $s_{d}+q^{2}s_{r}$
---|---|---|---|---|---|---|---|---|---|---
VMI | 10 | 30 |  all $e^{-}$ | total |  $-8.21\pm 1.21$ |  $60.349\pm 0.006$ |  $0.217\pm 0.005$ |  $0.800\pm 0.015$ |  $0.016\pm 0.004$ | $1.9\pm 0.4$
VMI | 10 | 30 |  all $e^{-}$ | EEL |  $-8.63\pm 1.32$ |  $60.354\pm 0.006$ |  $0.220\pm 0.005$ |  $0.201\pm 0.004$ |  $0.004\pm 0.001$ | $0.50\pm 0.12$
VMI | 10 | 30 |  all $e^{-}$ | ICD |  $0$ |  $60.415\pm 0.003$ |  $0.220\pm 0.009$ |  $0.011\pm 0.000$ |  $-0.022\pm 0.001$ | $0.0110\pm 0.0003$
VMI | 12 | 30 |  all $e^{-}$ | total |  $-2.36\pm 0.07$ |  $60.402\pm 0.005$ |  $0.233\pm 0.005$ |  $1.361\pm 0.011$ |  $0.272\pm 0.015$ | $2.88\pm 0.12$
VMI | 12 | 30 |  all $e^{-}$ | EEL |  $-3.10\pm 0.14$ |  $60.409\pm 0.006$ |  $0.220\pm 0.006$ |  $0.161\pm 0.002$ |  $0.023\pm 0.002$ | $0.38\pm 0.03$
VMI | 12 | 30 |  all $e^{-}$ | ICD |  $0$ |  $60.428\pm 0.006$ |  $0.183\pm 0.012$ |  $0.056\pm 0.002$ |  $-0.059\pm 0.002$ | $0.056\pm 0.002$
VMI | 14 | 30 |  all $e^{-}$ | total |  $-3.03\pm 0.18$ |  $60.365\pm 0.009$ |  $0.245\pm 0.010$ |  $4.441\pm 0.060$ |  $0.567\pm 0.064$ | $9.6\pm 0.9$
VMI | 14 | 30 |  all $e^{-}$ | EEL |  $-3.91\pm 0.27$ |  $60.401\pm 0.007$ |  $0.219\pm 0.008$ |  $0.413\pm 0.006$ |  $0.040\pm 0.005$ | $1.02\pm 0.11$
VMI | 14 | 30 |  all $e^{-}$ | ICD |  $0$ |  $60.439\pm 0.010$ |  $0.173\pm 0.019$ |  $0.196\pm 0.013$ |  $-0.210\pm 0.013$ | $0.196\pm 0.013$
VMI | 16 | 30 |  all $e^{-}$ | total |  $-2.34\pm 0.13$ |  $60.391\pm 0.014$ |  $0.276\pm 0.013$ |  $0.522\pm 0.008$ |  $0.121\pm 0.011$ | $1.18\pm 0.09$
VMI | 16 | 30 |  all $e^{-}$ | EEL |  $-2.41\pm 0.12$ |  $60.491\pm 0.011$ |  $0.243\pm 0.012$ |  $0.040\pm 0.001$ |  $0.011\pm 0.001$ | $0.104\pm 0.009$
VMI | 16 | 30 |  all $e^{-}$ | ICD |  $0$ |  $60.483\pm 0.014$ |  $0.147\pm 0.025$ |  $0.028\pm 0.003$ |  $-0.030\pm 0.003$ | $0.028\pm 0.003$

VMI & 10 30 He+ coinc. total $-3.39\pm 1.96$ $60.300\pm 0.043$ $0.344\pm
0.033$ $0.118\pm 0.005$ $0.007\pm 0.007$ $0.19\pm 0.12$ VMI 10 30 He+ coinc.
EEL $-24.56\pm 59.89$ $60.335\pm 0.032$ $0.263\pm 0.020$ $0.032\pm 0.004$
$0.000\pm 0.000$ $0.032\pm 0.004$ VMI 10 30 He+ coinc. ICD $0$ $60.418\pm
0.002$ $0.235\pm 0.008$ $0.002\pm 0.000$ $-0.004\pm 0.000$ $0.00182\pm
0.00004$ VMI 12 30 He+ coinc. total $-2.04\pm 0.12$ $60.435\pm 0.013$
$0.277\pm 0.015$ $1.130\pm 0.013$ $0.174\pm 0.018$ $1.85\pm 0.11$ VMI 12 30
He+ coinc. EEL $-4.07\pm 0.32$ $60.398\pm 0.009$ $0.255\pm 0.012$ $0.277\pm
0.006$ $0.029\pm 0.004$ $0.76\pm 0.10$ VMI 12 30 He+ coinc. ICD $0$ $60.420\pm
0.004$ $0.205\pm 0.009$ $0.108\pm 0.002$ $-0.113\pm 0.003$ $0.108\pm 0.002$
VMI 14 30 He+ coinc. total $-2.12\pm 0.33$ $60.457\pm 0.035$ $0.300\pm 0.038$
$2.001\pm 0.019$ $0.097\pm 0.025$ $2.44\pm 0.18$ VMI 14 30 He+ coinc. EEL
$-6.04\pm 1.09$ $60.376\pm 0.014$ $0.236\pm 0.018$ $0.482\pm 0.013$ $0.018\pm
0.006$ $1.1\pm 0.3$ VMI 14 30 He+ coinc. ICD $0$ $60.429\pm 0.006$ $0.204\pm
0.013$ $0.207\pm 0.007$ $-0.219\pm 0.007$ $0.207\pm 0.007$ VMI 16 30 He+
coinc. total $-3.41\pm 0.47$ $60.154\pm 0.010$ $0.040\pm 0.006$ $1.021\pm
0.097$ $0.298\pm 0.104$ $4.5\pm 1.5$ VMI 16 30 He+ coinc. EEL $-3.51\pm 0.30$
$60.453\pm 0.019$ $0.289\pm 0.021$ $0.123\pm 0.003$ $0.017\pm 0.003$ $0.33\pm
0.05$ VMI 16 30 He+ coinc. ICD $0$ $60.459\pm 0.013$ $0.200\pm 0.026$
$0.084\pm 0.005$ $-0.089\pm 0.005$ $0.084\pm 0.005$ VMI 10 30 He${}_{2}^{+}$
coinc. total $-10.11\pm 1.70$ $60.299\pm 0.005$ $0.242\pm 0.004$ $0.380\pm
0.005$ $0.004\pm 0.001$ $0.79\pm 0.17$ VMI 10 30 He${}_{2}^{+}$ coinc. EEL
$-10.76\pm 6.73$ $60.312\pm 0.020$ $0.258\pm 0.015$ $0.054\pm 0.003$ $0.001\pm
0.001$ $0.17\pm 0.19$ VMI 10 30 He${}_{2}^{+}$ coinc. ICD $0$ $60.417\pm
0.003$ $0.247\pm 0.012$ $0.003\pm 0.000$ $-0.006\pm 0.000$ $0.0025\pm 0.0001$
VMI 12 30 He${}_{2}^{+}$ coinc. total $-2.39\pm 0.10$ $60.384\pm 0.007$
$0.232\pm 0.007$ $4.210\pm 0.042$ $0.757\pm 0.056$ $8.5\pm 0.5$ VMI 12 30
He${}_{2}^{+}$ coinc. EEL $-3.19\pm 0.21$ $60.400\pm 0.009$ $0.240\pm 0.009$
$0.486\pm 0.007$ $0.060\pm 0.007$ $1.10\pm 0.11$ VMI 12 30 He${}_{2}^{+}$
coinc. ICD $0$ $60.432\pm 0.007$ $0.191\pm 0.013$ $0.158\pm 0.006$ $-0.166\pm
0.006$ $0.158\pm 0.006$ VMI 14 30 He${}_{2}^{+}$ coinc. total $-3.01\pm 0.26$
$60.359\pm 0.011$ $0.214\pm 0.012$ $6.410\pm 0.081$ $0.586\pm 0.095$ $11.7\pm
1.3$ VMI 14 30 He${}_{2}^{+}$ coinc. EEL $-4.54\pm 0.84$ $60.379\pm 0.018$
$0.236\pm 0.020$ $0.644\pm 0.015$ $0.030\pm 0.011$ $1.3\pm 0.3$ VMI 14 30
He${}_{2}^{+}$ coinc. ICD $0$ $60.455\pm 0.011$ $0.187\pm 0.022$ $0.237\pm
0.016$ $-0.256\pm 0.016$ $0.237\pm 0.016$ VMI 16 30 He${}_{2}^{+}$ coinc.
total $-2.06\pm 0.08$ $60.456\pm 0.010$ $0.242\pm 0.010$ $2.175\pm 0.027$
$0.528\pm 0.037$ $4.4\pm 0.2$ VMI 16 30 He${}_{2}^{+}$ coinc. EEL $-2.86\pm
0.23$ $60.478\pm 0.015$ $0.214\pm 0.016$ $0.202\pm 0.005$ $0.034\pm 0.005$
$0.48\pm 0.06$ VMI 16 30 He${}_{2}^{+}$ coinc. ICD $0$ $60.480\pm 0.014$
$0.152\pm 0.025$ $0.127\pm 0.013$ $-0.135\pm 0.012$ $0.127\pm 0.013$ VMI 10 30
He${}_{3}^{+}$ coinc. total $-4.23\pm 0.59$ $60.361\pm 0.011$ $0.261\pm 0.008$
$0.042\pm 0.001$ $0.002\pm 0.001$ $0.08\pm 0.02$ VMI 10 30 He${}_{3}^{+}$
coinc. EEL $-10.82\pm 20.34$ $60.340\pm 0.073$ $0.341\pm 0.060$ $0.003\pm
0.001$ $0.000\pm 0.000$ $0.003\pm 0.001$ VMI 10 30 He${}_{3}^{+}$ coinc. ICD
$0$ $60.416\pm 0.003$ $0.225\pm 0.011$ $0.000\pm 0.000$ $-0.000\pm 0.000$
$0.00023\pm 0.00001$ VMI 12 30 He${}_{3}^{+}$ coinc. total $-2.34\pm 0.11$
$60.412\pm 0.007$ $0.228\pm 0.008$ $0.569\pm 0.007$ $0.105\pm 0.009$ $1.14\pm
0.07$ VMI 12 30 He${}_{3}^{+}$ coinc. EEL $-2.29\pm 0.29$ $60.468\pm 0.021$
$0.239\pm 0.022$ $0.053\pm 0.002$ $0.010\pm 0.002$ $0.105\pm 0.017$ VMI 12 30
He${}_{3}^{+}$ coinc. ICD $0$ $60.440\pm 0.007$ $0.191\pm 0.013$ $0.017\pm
0.001$ $-0.018\pm 0.001$ $0.017\pm 0.001$ VMI 14 30 He${}_{3}^{+}$ coinc.
total $-2.93\pm 0.28$ $60.397\pm 0.013$ $0.219\pm 0.014$ $0.775\pm 0.011$
$0.073\pm 0.013$ $1.40\pm 0.16$ VMI 14 30 He${}_{3}^{+}$ coinc. EEL $-6.36\pm
3.73$ $60.412\pm 0.046$ $0.266\pm 0.050$ $0.066\pm 0.004$ $0.002\pm 0.002$
$0.15\pm 0.12$ VMI 14 30 He${}_{3}^{+}$ coinc. ICD $0$ $60.466\pm 0.012$
$0.191\pm 0.024$ $0.026\pm 0.002$ $-0.028\pm 0.002$ $0.026\pm 0.002$ VMI 16 30
He${}_{3}^{+}$ coinc. total $-1.95\pm 0.09$ $60.488\pm 0.012$ $0.249\pm 0.012$
$0.261\pm 0.004$ $0.073\pm 0.006$ $0.54\pm 0.03$ VMI 16 30 He${}_{3}^{+}$
coinc. EEL $-1.61\pm 0.17$ $60.588\pm 0.025$ $0.205\pm 0.026$ $0.018\pm 0.001$
$0.008\pm 0.001$ $0.039\pm 0.005$ VMI 16 30 He${}_{3}^{+}$ coinc. ICD $0$
$60.492\pm 0.014$ $0.151\pm 0.025$ $0.013\pm 0.001$ $-0.013\pm 0.001$
$0.013\pm 0.001$ HEA 14 50 all $e^{-}$ EEL $-2.94\pm 0.15$ $60.405\pm 0.010$
$0.223\pm 0.011$ $0.356\pm 0.006$ $0.063\pm 0.006$ $0.90\pm 0.08$ HEA 14 50
all $e^{-}$ ICD $0.01\pm 0.02$ $60.381\pm 0.015$ $0.24\pm 0.03$ $0.133\pm
0.008$ $-0.141\pm 0.007$ $0.133\pm 0.008$

##  IV Electron trapping probability

Trapping of slow ICD electrons in helium droplets with increasing size affects
the probability of their observation [4, 5]. It is modelled as a product of
the size independent ICD probability $p_{0}$ and the likelihood for electron
ejection before trapping (see fit in Fig. 3 in the main text),

$\displaystyle p^{\rm obs}_{\rm ICD}$
$\displaystyle=p_{0}\exp\left(-\frac{\bar{d}}{l_{\rm
slow}}\right)=p_{0}\exp\left(-\frac{3R}{4l_{\rm slow}}\right)$ (8)

where $\bar{d}=\frac{3}{4}R$ is the average travelling distance to the surface
in a droplet with radius $R$ and $l_{\rm slow}$ is the mean free path of the
electron. We have thus neglected the angular deflection by elastic and
inelastic collisions that slow down the electron and roughly estimate the
ejection probability alone from the average distance $\bar{d}$ from the
surface. The droplet radius is computed from the number of atoms $N_{\rm He}$
in the droplet by the expression

$\displaystyle R/\text{nm}$ $\displaystyle=0.22\sqrt[3]{N_{\rm He}}$ (9)

by Gomez et al. [6]. The expression for the mean distance $\bar{d}$ is derived
as follows. From Fig. 4 we find the distance $d$ of an electron from the
surface of a spherical droplet with radius $R$,

$\displaystyle d$
$\displaystyle=\sqrt{R^{2}-a^{2}\sin^{2}(\theta)}-a\cos(\theta)$ (10)

where $a$ is the distance from the electron to the center of the droplet. The
isotropic likelihood that the electron originates from a distance $a$ from the
center and moves toward a surface element $[\theta,\theta+{\rm
d}\theta]\times[\phi,\phi+{\rm d}\phi]$ is given by

$\displaystyle\frac{3a^{2}}{R^{3}}\,{\rm d}a\leavevmode\nobreak\
\times\leavevmode\nobreak\ \frac{1}{4\pi}\sin(\theta)\,{\rm d}\theta\,{\rm
d}\phi$ (11)

and we obtain the distance averaged over all positions and directions by
integration,

$\displaystyle\bar{d}$ $\displaystyle=\frac{3}{4\pi
R^{3}}\int_{0}^{R}\int_{0}^{\pi}\int_{0}^{2\pi}a^{2}\left(\sqrt{R^{2}-a^{2}\sin^{2}(\theta)}-a\cos(\theta)\right)\sin(\theta)\,{\rm
d}\phi\,{\rm d}\theta\,{\rm d}a$ (12)
$\displaystyle=\frac{3}{2R^{3}}\int_{0}^{R}\int_{0}^{\pi}\left(\sqrt{R^{2}-a^{2}\sin^{2}(\theta)}-a\cos(\theta)\right)a^{2}\sin(\theta)\,{\rm
d}\theta\,{\rm d}a\quad\text{(the $\cos(\theta)\sin(\theta)$ term vanishes by
symmetry)}$
$\displaystyle=\frac{3R}{2}\int_{0}^{1}\int_{0}^{\pi}\sqrt{1-k^{2}\sin^{2}(\theta)}\;k^{2}\sin(\theta)\,{\rm
d}\theta\,{\rm
d}k=\frac{3R}{2}\int_{0}^{1}\int_{-1}^{1}\sqrt{1-k^{2}\left(1-\mu^{2}\right)}\;k^{2}\,{\rm
d}\mu\,{\rm d}k$
$\displaystyle=\frac{3R}{2}\int_{0}^{1}k\int_{-k}^{k}\sqrt{x^{2}+1-k^{2}}\;\,{\rm
d}x\,{\rm
d}k=\frac{3R}{2}\int_{0}^{1}\frac{k}{2}\left[x\sqrt{x^{2}+1-k^{2}}+(1-k^{2})\ln\left(x+\sqrt{x^{2}+1-k^{2}}\right)\right]_{-k}^{k}\,{\rm
d}k$
$\displaystyle=\frac{3R}{2}\int_{0}^{1}\frac{k}{2}\left(k-(-k)+(1-k^{2})\ln\left(\frac{1+k}{1-k}\right)\right)\,{\rm
d}k=\frac{3R}{4}\int_{0}^{1}\left(2k^{2}+\left(k-k^{3}\right)\ln\left(\frac{1+k}{1-k}\right)\right)\,{\rm
d}k=\frac{3R}{4}\;.$

$R$$d$$a$$a\sin\theta$$a\cos\theta$$\theta$ Figure 4: Geometrical construction
of the distance $d$ of a point in a sphere from a point on its surface. The
distance of the inner point from the center of the sphere with radius $R$ is
denoted by $a$. One finds
$\left(d+a\cos\left(\theta\right)\right)^{2}=R^{2}-a^{2}\sin^{2}\left(\theta\right)$.

##  V Electron impact excitation probability

The loss of fast electrons is predominated by electron impact excitation (EIE)
and ionization (EII) that we describe in a simple model by the average
travelling distance $\bar{d}$ to the surface that has been derived in the
previous section and a mean free path $l_{\rm fast}$ in analogy to Eq. (8).
The EIE probability is then given by

$\displaystyle P_{\rm EIE}/\text{\small\%}$ $\displaystyle=p_{\rm
EIE}\left(1-\exp\left(-\frac{\bar{d}}{l_{\rm fast}}\right)\right)=p_{\rm
EIE}\left(1-\exp\left(-\frac{3R}{4l_{\rm fast}}\right)\right)$ (13)

where $p_{\rm EIE}$ is the fraction of EIE from the total inelastic scattering
cross section (EIE = EEL $-$ ICD). The fit in Fig. 5 (a) gives the mean free
path for inelastic scattering $l_{\rm fast}=(17.6\pm 1.4)$ nm and a $(40\pm
2)$ % fraction of EIE of the 1s2s/1s2p states relative to the total inelastic
cross section in close agreement with $17\text{\,}\mathrm{n}\mathrm{m}$ and
$38\text{\,}\mathrm{\char 37\relax}$ in Sec. VI derived from the cross
sections in Ref. [7]. Surprisingly, the loss of slow electrons thus seems to
be negligible in the investigated droplet size range, in opposition to the
size dependence of the ICD probability that indicates suppressed ejection of
slow electrons. This discrepancy is only partially explained by the fact that
fast electrons typically travel a substantial distance to the droplet surface
before losing kinetic energy by inelastic scattering. A consistent explanation
of the observed size dependences would thus require an improved model
considering the different pathways for electrons propagating in a droplet
instead of using a single average distance from the surface.

Figure 5: Droplet-size dependence for total electrons and coincidences with
He${}_{n}^{+}$ cations, respectively. The points for $n=2,3$ are somewhat
shifted to lower and higher droplet sizes for better visibility. (a) EIE
probability relative to total electrons (all e-) or coincidences
(He${}_{n}^{+}$, $n=2,3$) with the EIE model fitted to the total electron
data. (b) Relative ratio of ICD and EIE (electron impact excitation of
1s2s/1s2p states) counts. The reen solid line connects the total electron
points to guide the eye.

The ratios of EIE coincidences to total coincidences in Fig. 5 (a) do increase
with growing EEL probability at larger droplet sizes, but the growth is
significantly weaker than for all electrons because the cations are more
easily trapped in the droplets than electrons [4]. At all sizes, the ICD/EIE
ratio in Fig. 5 (b) is within uncertainties identical for total electrons and
e-–He${}_{n}^{+}$, $n=1,2,3$, coincidences. There is consequently no evidence
that the ICD process at short distances in the He∗He+ potential energy curves
affects the ion formation and ejection dynamics.

##  VI Electron impact cross section

To estimate the mean free path $l_{\rm fast}$ for inelastic scattering, we sum
up the electron impact ionization and excitation cross sections from the
initial 1 1S state at $E_{e}=$35.8\text{\,}\mathrm{e}\mathrm{V}$$ electron
energy (corresponding to $E_{h\nu}=$60.4\text{\,}\mathrm{e}\mathrm{V}$$ photon
energy) according to graphs 4 – 6 in Ref. [7] which gives the cross section
$\sigma_{\rm total}\approx$0.27\text{\,}\AA^{2}$$. With the helium density
$n_{\rm He}=N_{\rm He}/(4/3\pi
R^{3})\approx$22.42\text{\,}\mathrm{n}\mathrm{m}^{-3}$$ from Eq. (9) we obtain
the corresponding mean free path

$\displaystyle l_{\rm fast}=\frac{1}{\sigma_{\rm exc}\,n_{\rm
He}}\approx$17\text{\,}\mathrm{n}\mathrm{m}$.$ (14)

In the main text, EEL denotes all electrons with 18 to 23 eV energy loss,
excluding excitation of 1s3p and 1s4p states that overlap with electron impact
ionization. Subtracting the ICD contribution gives the EIE yield that
consequently only includes excitation of 1s2s and 1s2p states with the cross
section $\sigma_{\rm 1s2s/1s2p}\approx$0.10\text{\,}\AA^{2}$$ in Ref. [7]. Its
ratio to the total inelastic cross section is

$\displaystyle p_{\rm EIE}=\frac{\sigma_{\rm 1s2s/1s2p}}{\sigma_{\rm
total}}\approx 0.38\,.$ (15)

##  VII Theoretical estimate of the ICD probability

Total and partial decay widths of the He(2s2p)–He(1s2) ${}^{1}\Sigma_{u}$ and
${}^{1}\Pi_{u}$ states have been computed as described in the main text at
five different distances. The resulting autoionization (AI) probabilities are
summarized in Table 2. The differences of 100 % and the AI probabilities are
the total probabilities for the decay into different ICD channels as a
function of distance which is plotted in Fig. 6. We are not aware of available
data on the relative excitation of the ${}^{1}\Sigma_{u}$ and ${}^{1}\Pi_{u}$
states and the respective ICD probabilities are not differing much.
Particularly at larger distances, they are close to zero such that the
computational method can only give an estimate with limited accuracy.
Therefore, we average the probabilities at each distance with equal weight and
the underlying assumption of similar cross sections has only a small effect on
the resulting average.

Table 2: Likelihood of AI for ${}^{1}\Sigma_{u}$ and ${}^{1}\Pi_{u}$ states at different pair distances $r$.  r (a. u.) | 6 | 5.5 | 5 | 4.5 | 4
---|---|---|---|---|---
r (Å) | 3.18 | 2.91 | 2.65 | 2.38 | 2.12
He(2s2p) – He(1s2) ${}^{1}\Sigma_{u}$ $\rightarrow$ AI (%) | 98 | 97 | 94 | 92 | 85
He(2s2p) – He(1s2) ${}^{1}\Pi_{u}$ $\rightarrow$ AI (%) | 98 | 96 | 99 | 95 | 91

Figure 6: Computed ICD probabilities for the He(2s2p) – He(1s2)
${}^{1}\Sigma_{u}$ and ${}^{1}\Pi_{u}$ states with an exponential fit to their
average (magenta squares). Green circles represent the pair-distance
distribution function $g(r)$ [8, 9].

To derive the overall ICD probability, the distance dependent average
probabilities are fitted with an exponential curve. This curve is weighted
with the pair-distance distribution function $g(r)$ [8, 9] to obtain the total
ICD probability relative to the resonant double excitation cross section (AI
plus ICD) which is $(3.0\pm 2.8)\,\%$. The large uncertainty is due to the
standard deviations of the fit parameters $a=670\pm 380$, $b=1.9\pm 0.3$,
$p(r)=a\exp\left(-br\right)$. The ICD probability relative to the total
ionization cross section is obtained by multiplying this result with the ratio
$f_{\text{res}}$ of the resonant cross section to the total electron counts
$c(h\nu)$ in Eq. (5) evaluated at the double excitation resonance
($h\nu=E_{r}$),

$f_{\text{res}}=\frac{s_{r}q^{2}}{s_{d}+s_{r}q^{2}}\left(1+\frac{2}{\pi
q^{2}}\right)^{-1}$ (16)

which is a product of two ratios. According to Ref. [10], we identify a
constant background contribution $s_{d}$ and the unperturbed continuum cross
section $s_{r}$ in the Fano profile in Eq. (5). Then the total cross section
$c(h\nu)$ is a sum of the background $s_{d}$ and a transition probability —
which we denote as $c_{r}(h\nu)$ — that stems from the Fano interference of
the continuum with the discrete level (double excitation leading to AI or
ICD). At the resonance, the latter simplifies to $s_{r}q^{2}$ and its ratio to
the total cross section, $c_{r}(E_{r})/c(E_{r})$, is the first factor in Eq.
(16). For the second factor, we use the identification of $\frac{1}{2}\pi
q^{2}$ in Ref. [11] as the ratio of the transition probabilities to the
modified discrete level ($A$) and to a band width of unperturbed continuum
states ($B$). Assuming that the sum $A+B$ approximates the transition
probability $c_{r}(h\nu)$, its fraction of the resonant pathway is

$\frac{A}{c_{r}(h\nu)}\approx\frac{A}{A+B}=\frac{A/B}{A/B+1}=\frac{\frac{1}{2}\pi
q^{2}}{\frac{1}{2}\pi q^{2}+1}=\frac{1}{1+\frac{2}{\pi q^{2}}}$ (17)

which is the second factor needed to compute $f_{\text{res}}$. With the fit
results for total electrons in Table III we obtain the ratios in Table 3 for
four different nozzle temperatures. The values do not significantly differ
from each other and we use the average $f_{\text{res}}=0.51\pm 0.03$ as final
result.

Table 3: Fractions of resonant to total ionization according to Eq. 16 for VMI data at different nozzle temperatures.  Nozzle temperature | 10 K | 12 K | 14 K | 16 K
---|---|---|---|---
$f_{\text{res}}$ |  $0.57\pm 0.09$ |  $0.47\pm 0.02$ |  $0.50\pm 0.04$ | $0.49\pm 0.04$

## References

  * LaForge  _et al._ [2016] A. C. LaForge, D. Regina, G. Jabbari, K. Gokhberg, N. V. Kryzhevoi, S. R. Krishnan, M. Hess, P. O’Keeffe, A. Ciavardini, K. C. Prince, R. Richter, F. Stienkemeier, L. S. Cederbaum, T. Pfeifer, R. Moshammer, and M. Mudrich, Phys. Rev. A 93, 050502 (2016).
  * Kramida  _et al._ [2013] A. Kramida, Y. Ralchenko, J. Reader, and N. A. Team, NIST Atomic Spectra Database (version 5.1) (National Institute of Standards and Technology, Gaithersburg, MD) (2013).
  * Hertel and Schulz [2008] I. V. Hertel and C.-P. Schulz,  _Atoms, molecules and optical physics 1. Atomic physics and foundations of spectroscopy_ (Springer-Verlag Berlin Heidelberg, 2008).
  * Asmussen  _et al._ [2023a] J. D. Asmussen, L. Ben Ltaief, K. Sishodia, A. R. Abid, B. Bastian, S. Krishnan, H. B. Pedersen, and M. Mudrich, J. Chem. Phys. 159, 034301 (2023a).
  * Asmussen  _et al._ [2023b] J. D. Asmussen, K. Sishodia, B. Bastian, A. R. Abid, L. Ben Ltaief, H. B. Pedersen, S. De, C. Medina, N. Pal, R. Richter, T. Fennel, S. Krishnan, and M. Mudrich, Nanoscale 15, 14025 (2023b).
  * Gomez  _et al._ [2011] L. F. Gomez, E. Loginov, R. Sliter, and A. F. Vilesov, J. Chem. Phys. 135, 154201 (2011).
  * Ralchenko  _et al._ [2008] Y. Ralchenko, R. Janev, T. Kato, D. Fursa, I. Bray, and F. De Heer, At. Data Nucl. Data Tables 94, 603 (2008).
  * Chin and Krotscheck [1995] S. A. Chin and E. Krotscheck, Phys. Rev. B 52, 10405 (1995).
  * Peterka  _et al._ [2007] D. S. Peterka, J. H. Kim, C. C. Wang, L. Poisson, and D. M. Neumark, J. Phys. Chem. A 111, 7449 (2007).
  * Mezei  _et al._ [2015] J. Mezei, I. F. Schneider, and C. Jungen, EPJ Web of Conferences 84, 04005 (2015).
  * Fano [1961] U. Fano, Phys. Rev. 124, 1866 (1961).


[]Tamás G. Kovács

# High temperature $U(1)_{A}$ breaking in the chiral limit

###### Abstract

We solve the long-standing problem concerning the fate of the chiral
$U(1)_{A}$ symmetry in QCD-like theories at high temperature in the chiral
limit. We introduce a simple instanton based random matrix model that
precisely reproduces the properties of the lowest part of the lattice overlap
Dirac spectrum. We show that in the chiral limit the instanton gas splits into
a free gas component with a density proportional to $m^{N_{f}}$ and a gas of
instanton-antiinstanton molecules. While the latter do not influence the
chiral properties, for any nonzero quark mass the free gas component produces
a singular spectral peak at zero that dominates Banks-Casher type spectral
sums. By calculating these we show that the difference of the pion and delta
susceptibility vanishes only for three or more massless flavors, however, the
chiral condensate is zero already for two massless flavors.

##  1 Introduction

The approximate $SU(2)_{A}\times U(1)_{A}$ chiral symmetry of light quarks is
an essential feature of quantum-chromodynamics (QCD). While the $SU(2)_{A}$
part, spontaneously broken at low temperatures, is expected to be restored
above the crossover temperature $T_{c}$, the flavor singlet $U(1)_{A}$ part is
anomalous due to the presence of instantons. The fate of the flavor singlet
part is the subject of a long debate that was started by the seminal paper of
Pisarski and Wilczek, who discussed possible scenarios for the finite
temperature chiral transition in the chiral limit [1]. Their results were
based on the $\epsilon$-expansion, and there has been an ongoing effort to
check it using lattice QCD.

Figure 1: The distribution of the lowest Dirac eigenvalue on configurations
containing exactly one instanton-antiinstanton pair on the lattice and in the
random matrix model in a spatial volume of $32^{3}$ in lattice units (left
panel). The same comparison for the second lowest eigenvalue on a larger
volume of $40^{3}$ (right panel). For comparison we also show the distribution
used for the fit, the one in the previous figure.

In the chiral limit the order parameter of the spontaneously broken chiral
symmetry, $\langle\bar{\psi}\psi\rangle$ can be written in terms of the
spectral density $\rho(\lambda)$ of the quark Dirac operator as

$\lim_{m\rightarrow 0}\langle\bar{\psi}\psi\rangle\propto\rho(0).$

This is the Banks-Casher relation [2], showing that the low lying spectrum of
the Dirac operator is intimately related to how chiral symmetry is realized.
Indeed, for some time the standard lore was that at the transition temperature
to the quark-gluon plasma, the spontaneously broken chiral symmetry is
restored, and $\rho(0)$ vanishes, up to small explicit breaking effects due to
the finite, but small light quark masses.

However, this view was challenged when instead of a small $\rho(0)$, a sharp
rise was found in the spectral density of the overlap Dirac operator at zero
[3]. The reason this spike in the spectral density went unnoticed before, was
that unlike the overlap, the staggered and Wilson lattice Dirac operators,
used exclusively previously, did not respect chiral symmetry, and could not
properly resolve the smallest Dirac eigenvalues, the ones that make up the
spike of the density. For some time this finding remained largely ignored,
mostly because it was considered to be a quenched artifact, the result of
ignoring the quark determinant in the path integral. Indeed, the determinant
disfavors eigenvalues of small magnitude, and is expected to suppress the
spike in the spectral density at zero. Later on, the spike was found to
persist even in the presence of light dynamical quarks [4, 6, 5, 7]. However,
some doubts could still remain, as the lattice fermions used in these works
were not chirally symmetric, and their poor resolution of the small Dirac
eigenvalues might have lead to an improper suppression of the spike. Indeed,
results by the JLQCD collaboration using chiral lattice fermions suggest that
toward the chiral limit the spectral spike completely disappears already at a
nonzero quark mass [9, 8]. Another possibility, the one we advocate in the
present paper, is that indeed, chiral fermions are needed for a proper
suppression of the spike, but the spike is still present at any finite quark
mass, however, it becomes undetectably small with the currently used lattice
volumes and statistics.

It seems to us that using the currently available lattice technology it is not
possible to explore the fate of the spike in the spectral density and also the
fate of the $U(1)_{A}$ symmetry for light dynamical quarks. In the present
work we suggest a different approach, based on the finding that in the
quenched theory the statistics of the eigenvalues in the spectral spike are
consistent with mixing zero modes of a free instanton gas, a proposal already
put forward in Ref. [3] and recently confirmed in more detail [10].

##  2 Random matrix model

We propose a random matrix model for the description of the mixing zero modes,
the zero mode zone (ZMZ) of the free instanton gas in quenched QCD. In a free
instanton gas the number distribution of instantons and antiinstantons are
independent and identical Poisson distributions with mean $\chi_{0}V/2$, where
$\chi_{0}$ is the topological susceptibility. A random matrix is constructed
by first drawing the number of instantons $n_{\mbox{\scriptsize i}}$ and
antiinstantons $n_{\mbox{\scriptsize a}}$ from the Poisson distributions. The
size of the matrix is $(n_{\mbox{\scriptsize i}}+n_{\mbox{\scriptsize
a}})\times(n_{\mbox{\scriptsize i}}+n_{\mbox{\scriptsize a}})$. Since zero
modes of like charge objects are protected from mixing by the index theorem,
the matrix has two diagonal blocks of zeros of size $n_{\mbox{\scriptsize
i}}\times n_{\mbox{\scriptsize i}}$ and $n_{\mbox{\scriptsize a}}\times
n_{\mbox{\scriptsize a}}$. At high temperatures the instanton zero mode wave
functions decay exponentially with exponent $\pi T$ [11], and we assume that
the mixing matrix elements of instanton and antiinstanton zero modes also
follow this exponential decay with the distance of the given instanton and
antiinstanton. In a noninteracting gas, the location of the topological
objects are chosen randomly within a three-dimensional111At high temperature
the instantons typically occupy the whole available space in the temporal
direction, so we ignore the temporal location of the instantons. box of size
$L$. In this way the remaining off-diagonal blocks of the matrix are filled
with elements of the form $w_{kl}=A\cdot\exp(-\pi Tr_{kl})$, where $r_{kl}$ is
the distance of the randomly placed instanton $k$ and antiinstanton $l$. This
completes the construction of one random matrix, and we can easily generate
ensembles of such random matrices.

Figure 2: The distribution of the lowest Dirac eigenvalue on the lattice
ensemble described in the text, compared to that in the matrix model with the
best fit parameter $A$ (left panel). In both cases we used only the
configurations with exactly one instanton-antiinstanton pair , and to better
resolve the spike in the spectral density, we plotted the distribution of the
natural log of the eigenvalues. The right panel shows a similar comparison,
but with the second lowest eigenvalue, and on a larger lattice with $L=40$.

The model has two parameters, the topological susceptibility $\chi_{0}$ and
the prefactor $A$. To determine these parameters we consider an ensemble of
20k quenched lattice configurations generated with the Wilson action at
$T=1.11T_{c}$ on lattices of size $32^{3}\times 8$. Computing the lowest
eigenvalues of the overlap Dirac operator on this ensemble and counting the
exact zero eigenvalues yields the topological susceptibility. For fitting the
single remaining parameter $A$, we consider the distribution of the lowest
nonzero Dirac eigenvalue on those configurations that have exactly two complex
conjugate eigenvalues in the ZMZ, corresponding to an instanton-antiinstanton
pair. In Fig. 2 (left panel) we compare the distribution of the lattice
eigenvalues with that of the random matrix model with the best fit parameter
$A$. It is already remarkable that the whole distribution can be fitted with
just this one parameter, but now that the model is fixed we can compare
different properties of the Dirac spectrum on the lattice and in the matrix
model. The comparison can be made also for the distribution of the second
lowest eigenvalue or the full spectral density, and also on different volumes.
As an illustration, in Fig. 2 (right panel) we show the comparison for the
second lowest eigenvalue on a larger volume. These tests demonstrate that the
random matrix model properly describes the details of the lattice overlap
spectrum. Simulations of the random matrix model on larger volumes, not
accessible to direct lattice simulations, indicate the in the thermodynamic
limit the spectral density is singular at zero [12].

##  3 Including dynamical quarks

On the lattice, including dynamical quarks means that in addition to the
quenched Boltzmann weights, each configuration gets another weight factor
proportional to the quark Dirac determinant $\det(D+m)^{N_{f}}$. For
simplicity we assume $N_{f}$ degenerate quark flavors. The determinant of the
lattice Dirac operator can be split into the contribution of the zero mode
zone and that of the bulk as

$\det(D+m)=\prod_{i\in ZMZ}(\lambda_{i}+m)\cdot\prod_{i\in\mbox{\tiny
bulk}}(\lambda_{i}+m).$ (1)

As can be seen in Fig. 1, at high temperatures the ZMZ and the bulk are well
separated, therefore eigenvalues in the bulk are not expected to be correlated
with the ones in the ZMZ. Our main concern here is to study how the
determinant suppresses the eigenvalues in the ZMZ. It is thus a good
approximation to ignore the bulk contribution to the determinant, as that will
only give a trivial factor in the path integral for quantities depending on
the ZMZ. Especially for small quark mass, the suppression of the eigenvalues
in the ZMZ will be driven by the contribution of those small eigenvalues to
the determinant. This is exactly the part of the Dirac spectrum that is
faithfully represented by the random matrix model, so including the
determinant of the random matrices in the statistical weight will properly
describe the suppression of the spectral spike by dynamical quarks. Thus the
random matrix model we propose for the ZMZ of QCD with dynamical quarks is the
one described in the previous section, with the additional weight
$\det(D+m)^{N_{f}}$ for each instanton configuration, where $D$ is the random
matrix corresponding to the given instanton configuration.

The numerical simulation of this model reveals that the topological
susceptibility is suppressed by light dynamical quarks as

$\chi(m)=\chi_{0}m^{N_{f}},$ (2)

where $\chi_{0}$ is the quenched susceptibility, the one obtained without
including the quark determinant. To understand this behavior we note that
throughout the simulations we found that the eigenvalues in the spectral spike
always satisfied $|\lambda|\ll m$, even for the smallest quark mass of
$m=0.01$ (in lattice units). The reason for this is that as the quarks become
lighter, the determinant suppresses the number of instantons, the typical
instanton-antiinstanton distance grows, and the matrix elements become
exponentially small, resulting in ever smaller eigenvalues. As a result, in
the chiral limit the magnitude of the eigenvalues in the ZMZ decreases much
faster than the quark mass, and they will always remain smaller than the quark
mass. If the eigenvalues are much smaller than the quark mass then to a very
good approximation the determinant of a matrix with $n_{\mbox{\scriptsize i}}$
instantons and $n_{\mbox{\scriptsize a}}$ antiinstantons can be written as

$\det(D+m)=\prod_{i}(\lambda_{i}+m)\approx m^{n_{\mbox{\scriptsize
i}}+n_{\mbox{\scriptsize a}}}.$ (3)

Figure 3: The spectral density of the random matrix model with two degenerate
dynamical quark flavors of mass $m=0.05$ (in lattice units) compared to the
spectral density of the matrix model of the free instanton gas with the same
topological susceptibility (left panel). A comparison of the density of
nearest opposite charged topological objects between the quenched and the
dynamical random matrix ensemble (right panel).

With $N_{f}$ quark flavors, each (anti)instanton contributes a suppression
factor $m^{N_{f}}$ to the determinant, and the distribution of (anti)instanton
numbers is still Poissonian, but with a susceptibility suppressed as
$\chi_{0}\rightarrow\chi_{0}m^{N_{f}}$. This explains the quark mass
dependence of the susceptibility in Eq. (2). The fact that the Poisson
distribution of the number of topological objects is preserved also implies
that even in the presence of light dynamical quarks, the lowest part of the
Dirac spectrum can still be understood as the zero mode zone of a free
instanton gas. To demonstrate this, in the left panel of Fig. 3 we compare the
spectral density of the random matrix model with two degenerate dynamical
quark flavors of mass $m=0.05$ (in lattice units) with that of the matrix
model of the free instanton gas (without the determinant) with the same
topological susceptibility. The two curves exactly agree, except for the
largest eigenvalues, where the model with dynamical quarks shows an excess of
eigenvalues. Large eigenvalues in the matrix model indicate that there might
be large matrix elements which in turn would imply that in the dynamical case
there are nearby instanton-antiinstanton pairs. To check that, in the right
panel of Fig. 3 we compare the density of the nearest opposite charged objects
as a function of their distance for the dynamical and quenched matrix
ensembles of the left panel of the same figure. Indeed, in the dynamical case
we find an excess of tightly bound instanton-antiinstanton pairs at a distance
smaller than the instanton size $1/T$. These instanton-antiinstanton
“molecules” are held together by the attractive force due to light dynamical
quarks.

##  4 Chiral condensate and $U(1)_{A}$ breaking susceptibility

We have seen that with light dynamical quarks the instanton gas has two
components. There is a dilute gas of free instantons, responsible for the
small Dirac eigenvalues, i.e. the spectral spike. Besides that, there is a
component of tightly bound instanton-antiinstanton molecules. Our simulations
reveal that in the chiral limit the eigenvalues corresponding to these two
components behave differently. While the free instanton eigenvalues in the
spike become smaller as the gas becomes more dilute, the eigenvalues
corresponding to the molecules remain at a fixed scale in the spectrum. This
is because the size of the molecules does not change with the quark mass. An
important consequence of this is that in the chiral limit any nonvanishing
contribution to Banks-Casher type sums can only come from the free instanton
generated part of the spectrum.

The simplest of these sums is the one providing the chiral condensate. In the
chiral limit the chiral condensate can be written in terms of the Dirac
spectrum as

$\langle\bar{\psi}\psi\rangle\approx\langle\sum_{i}\frac{m}{m^{2}+\lambda_{i}^{2}}\rangle\approx\underbrace{\left(\mbox{\parbox{43.05542pt}{\tiny
avg.\ number of instantons in free
gas}}\right)}_{\chi_{0}m^{N_{\mbox{\scriptsize
f}}}V}\cdot\frac{1}{m}\;=\;m^{N_{\mbox{\scriptsize
f}}-1}\chi_{\mbox{\scriptsize 0}}V.$ (4)

Here we used the fact that the eigenvalues corresponding to the molecular
component of the instanton gas remain at a fixed scale, and in the chiral
limit they do not contribute to the sum. In contrast, the magnitude of
eigenvalues generated by the free gas component becomes smaller in the chiral
limit so rapidly that $\lim_{m\rightarrow 0}\frac{\lambda}{m}=0$, and all
these eigenvalues will contribute a term $1/m$ to the sum. So the chiral limit
of the condensate is given by the contribution of the free instanton gas. This
result shows that for two and more flavors of vanishing mass, the chiral
condensate goes to zero, and the spontaneously broken chiral symmetry is
restored, as expected. We also note that Eq. (4) is consistent with the
expansion of the QCD free energy density in terms of the quark mass around the
chiral point [13], and also consistent with the resulting quasi-instanton
picture of [14]. In fact, our instanton-based matrix model provides a physical
explanation for the quasi-instanton picture, and shows that it is valid not
only in the small quark mass limit.

Another quantity of interest in the chiral limit is the susceptibility
$\chi_{\pi}-\chi_{\delta}$, a nonzero value of which signals the breaking of
the $U(1)_{A}$ symmetry [6, 5, 9]. Similarly to the chiral condensate, in the
chiral limit this quantity can also be written as a sum over the Dirac
spectrum as

$\chi_{\pi}-\chi_{\delta}\approx\langle\sum_{i}\frac{m^{2}}{(m^{2}+\lambda_{i}^{2})^{2}}\rangle\approx
m^{N_{\mbox{\scriptsize f}}}\chi_{\mbox{\scriptsize
0}}V\cdot\frac{1}{m^{2}}\;=\;m^{N_{\mbox{\scriptsize
f}}-2}\chi_{\mbox{\scriptsize 0}}V,$ (5)

showing that even for two flavors of vanishing mass, the symmetry remains
effectively broken.

## Acknowledgments

This work was supported by Hungarian National Research, Development and
Innovation Office NKFIH Grant No. KKP126769 and NKFIH excellence grant
TKP2021-NKTA-64.

## References

  * [1] R. D. Pisarski and F. Wilczek, Phys. Rev. D 29, 338-341 (1984) doi:10.1103/PhysRevD.29.338
  * [2] T. Banks and A. Casher, Nucl. Phys. B 169, 103-125 (1980) doi:10.1016/0550-3213(80)90255-2
  * [3] R. G. Edwards, U. M. Heller, J. E. Kiskis and R. Narayanan, Phys. Rev. D 61, 074504 (2000) doi:10.1103/PhysRevD.61.074504.
  * [4] A. Alexandru and I. Horváth, Phys. Rev. D 92, no.4, 045038 (2015) doi:10.1103/PhysRevD.92.045038.
  * [5] H. T. Ding, S. T. Li, S. Mukherjee, A. Tomiya, X. D. Wang and Y. Zhang, Phys. Rev. Lett. 126, no.8, 082001 (2021) doi:10.1103/PhysRevLett.126.082001.
  * [6] O. Kaczmarek, L. Mazur and S. Sharma, Phys. Rev. D 104, no.9, 094518 (2021) doi:10.1103/PhysRevD.104.094518.
  * [7] X. L. Meng et al. [$\chi$QCD and CLQCD], [arXiv:2305.09459 [hep-lat]]. [8]
  * [8] S. Aoki et al. [JLQCD], PTEP 2022, no.2, 023B05 (2022) doi:10.1093/ptep/ptac001
  * [9] S. Aoki et al. [JLQCD], Phys. Rev. D 103, no.7, 074506 (2021) doi:10.1103/PhysRevD.103.074506.
  * [10] R. A. Vig and T. G. Kovacs, Phys. Rev. D 103, no.11, 114510 (2021) doi:10.1103/PhysRevD.103.114510.
  * [11] D. J. Gross, R. D. Pisarski and L. G. Yaffe, Rev. Mod. Phys. 53, 43 (1981) doi:10.1103/RevModPhys.53.43.
  * [12] T. G. Kovacs, [arXiv:2311.04208 [hep-lat]].
  * [13] T. Kanazawa and N. Yamamoto, JHEP 01, 141 (2016) doi:10.1007/JHEP01(2016)141.
  * [14] T. Kanazawa and N. Yamamoto, Phys. Rev. D 91, 105015 (2015) doi:10.1103/PhysRevD.91.105015.