#  Separation of Infrared and Bulk in Thermal QCD

Xiao-Lan Meng1,2, Peng Sun3, Andrei Alexandru4, Ivan Horváth5,6, Keh-Fei Liu6,
Gen Wang7, and Yi-Bo Yang1,2,8,9 ($\chi$QCD & CLQCD Collaboration) 1CAS Key
Laboratory of Theoretical Physics, Institute of Theoretical Physics, Chinese
Academy of Sciences, Beijing 100190, China
2University of Chinese Academy of Sciences, School of Physical Sciences,
Beijing 100049, China
3Institute of Modern Physics, Chinese Academy of Sciences, Lanzhou, 730000,
China
4The George Washington University, Washington, DC 20052, USA
5Nuclear Physics Institute CAS, 25068 Rez (Prague), Czech Republic
6University of Kentucky, Lexington, KY 40506, USA
7Aix-Marseille Université, Université de Toulon, CNRS, CPT, Marseille, France
8School of Fundamental Physics and Mathematical Sciences, Hangzhou Institute
for Advanced Study, UCAS, Hangzhou 310024, China
9International Centre for Theoretical Physics Asia-Pacific, Beijing/Hangzhou,
China

###### Abstract

New thermal phase of QCD, featuring scale invariance in the infrared (IR), was
proposed to exist both in the pure-glue (Nf=0) and the “real-world” (Nf=2+1)
settings. Among key aspects of the proposal is that the system in this IR
phase separates into two independent components: the scale-invariant IR part
and the non-invariant bulk. Such scenario requires non-analyticities in the
theory and, in case of pure-glue QCD, they were found to arise via Anderson-
like mobility edges in Dirac spectra ($\lambda_{\scriptscriptstyle{\text{\rm
IR}}}\\!=\\!0$, $\pm\lambda_{\text{A}}\\!\neq\\!0$) manifested in dimension
function $d_{\scriptscriptstyle{\text{\rm IR}}}(\lambda)$. Here we present
first evidence that this mechanism is also at work in “real-world QCD” (Nf=2+1
theory at physical quark masses and $a\\!=\\!0.105\,$fm), supporting the
existence of the proposed IR regime. Dimensional jump between zero-modes and
lowest near-zero modes very close to unity ($d_{\scriptscriptstyle{\text{\rm
IR}}}\\!=\\!3$ vs $d_{\scriptscriptstyle{\text{\rm IR}}}\\!\simeq\\!2$) was
found.

1\. Introduction: Starting with the pre-QCD times of Hagedorn Hagedorn (1965,
1971) and early lattice QCD calculations in the pure-glue setting McLerran and
Svetitsky (1981); Kuti  _et al._ (1981); Engels  _et al._ (1981), the question
of thermal transition in strongly-interacting matter has become one of the
highly researched topics in nuclear and particle physics. Apart from well-
motivated need to understand strong interactions, the interest in the issue
was fueled, to a large extent, by the potential significance of its resolution
to the physics of early universe.

Hagedorn in fact set the basic scenario, wherein the thermal transformation
process in strong interactions boils down to a single “instability
temperature” which, nowadays in QCD, is commonly referred as the critical
temperature ($T_{c}$). Due to the non-perturbative nature of the problem,
lattice QCD became the workhorse for investigations in this area. Advances in
lattice QCD techniques and computational resources led to a major conclusion,
namely that true phase transition does not occur in “real-world” QCD. Rather,
an analytic crossover takes place in the temperature range 150-200 MeV, with
$T_{c}\approx 155$ MeV for the case of chiral crossover Aoki  _et al._ (2006).

Transitionless outlook meant a setback to QCD’s role in cosmology, but an
important new twist appeared around the same time. Experiments at RHIC and LHC
concluded that the state of strongly interacting matter with properties akin
to near-perfect fluid exists in certain range of temperatures Arsene  _et al._
(2005); Back  _et al._ (2005); Adams  _et al._ (2005); Adcox  _et al._ (2005);
Muller  _et al._ (2012). Among other things, this invites questions about how
can such an exotic state of matter arise without a true phase transition.

To this end, some of us presented evidence of an unusual change in QCD Dirac
spectra at temperatures well above the chiral crossover Alexandru and Horváth
(2019): the anomalous accumulation of infrared (IR) Dirac modes, first seen in
high-$T$ phase of pure-glue QCD Edwards  _et al._ (2000) and shown to persist
into the continuum and infinite-volume limits Alexandru and Horváth (2015),
dramatically increases and starts to exhibit signs of scale-invariant
behavior. This sharp change was found in both the pure-glue and real-world
($N_{f}$ = 2+1 at physical masses) QCD Alexandru and Horváth (2019). It was
proposed that, at the associated temperature $T_{\scriptscriptstyle{\text{\rm
IR}}}$, thermal state of strong interactions reconfigures by forming two
independent components separated by new scale
$\Lambda_{{\scriptscriptstyle{\text{\rm IR}}}}(T)\lesssim T$: the bulk
governing distances $\ell\\!<\\!1/\Lambda_{\scriptscriptstyle{\text{\rm IR}}}$
and the IR part describing $\ell\\!>\\!1/\Lambda_{\scriptscriptstyle{\text{\rm
IR}}}$ via scale-invariant glue Alexandru and Horváth (2019). In pure-glue
case, $T_{\scriptscriptstyle{\text{\rm IR}}}$ coincides with $T_{c}$ of
Polyakov-line phase transition. In real-world QCD, it was also proposed to be
a true phase transformation occurring at
$200\\!<\\!T_{\scriptscriptstyle{\text{\rm IR}}}\\!<\\!250$ MeV Alexandru and
Horváth (2019). Its presence may clarify the physics of near-perfect fluid and
enhance the role of “QCD transition” in cosmology.

The 2-component scenario was first evoked by a clean IR-bulk separation in
Dirac spectra (bimodality of mode density $\rho(\lambda)$), which is very
suggestive of decoupling Alexandru and Horváth (2019). But more detail was
needed. Indeed, how would the scale invariant and non-invariant physics
coexist and would it imply a non-analytic running of the coupling constant at
$\Lambda_{\scriptscriptstyle{\text{\rm IR}}}$? Concrete proposal was presented
in Refs. Alexandru and Horváth (2021, 2022), ascribing the origin of non-
analyticity to two Anderson-like mobility edges (critical points) in Dirac
spectra. The first one at $\lambda_{\scriptscriptstyle{\text{\rm A}}}>0$ was
found previously Garcia-Garcia and Osborn (2006, 2007); Kovacs and Pittler
(2010); Giordano  _et al._ (2014); Ujfalusi  _et al._ (2015), and its purpose
here is to shield the IR component from the influence of the bulk. Indeed,
bulk fluctuations (including UV) will not affect the IR component owing to the
intervening non-analyticity. The second mobility edge was found recently
Alexandru and Horváth (2022), and is strictly IR
($\lambda_{\scriptscriptstyle{\text{\rm IR}}}\\!=\\!0$). Its role is to
facilitate the long-range physics of the IR component.

A suitable tool to express this scenario is the function
$d_{\scriptscriptstyle{\text{\rm IR}}}(\lambda)$, namely the spatial IR
dimension of Dirac modes at eigenvalue $\lambda$ Alexandru and Horváth (2021).
Indeed, $d_{{\scriptscriptstyle{\text{\rm IR}}}}$ is a proper dimensional
construct to probe the infrared Horváth  _et al._ (2023). The key conceptual
step granting its uses in quantum theory is the assignment of a meaningful
measure-based dimension to a region effectively selected by probabilities.
This has recently become possible via effective number theory Horváth and
Mendris (2020); Horváth (2021); Horváth and Mendris (2019): replacing ordinary
counting in definition of box-counting dimension for fixed sets by effective
counting for probability distributions leads to such measure-based dimension
Horváth  _et al._ (2023). For Dirac modes in thermal QCD the prescription is
as follows. In lattice-regularized Euclidean setup the number of sites
$N(L)\equiv(L/a)^{3}/(aT)$ (UV cutoff $1/a$, IR cutoff $1/L$, temperature $T$)
grows as $L^{3}$ at fixed $a$, conveying that IR dimension of space is
$D_{\scriptscriptstyle{\text{\rm IR}}}\\!=\\!3$. But Dirac eigenmode
$D\psi(x)\\!=\\!\lambda\psi(x)$ entails probabilities
$P\\!=\\!(p_{1},p_{2},\ldots,p_{N})$,
$p_{i}\\!\equiv\\!\psi^{+}(x_{i})\psi(x_{i})$, and sites have to be counted
effectively in order to quantify the volume $\psi$ actually extends into,
namely Horváth and Mendris (2020)

$N\,\longrightarrow\,\mathscr{N}_{\star}[\psi]\,=\,\mathscr{N}_{\star}[P]\,=\,\sum_{i=1}^{N}\min\,\\{Np_{i},1\\}.\;\,$
(1)

The IR scaling of QCD-averaged effective volume at given $\lambda$ then
determines $d_{\scriptscriptstyle{\text{\rm IR}}}(\lambda)$ at UV cutoff $a$,
namely Alexandru and Horváth (2021); Horváth  _et al._ (2023)

$\langle\,\mathscr{N}_{\star}\,\rangle_{L,\lambda,a}\,\propto\,L^{d_{\scriptscriptstyle{\text{\rm
IR}}}(\lambda,a)}\quad\text{for}\quad L\to\infty.\quad$ (2)

Using the overlap Dirac operator Neuberger (1998) due to its superior chiral
and topological properties, an unusual $d_{\scriptscriptstyle{\text{\rm
IR}}}(\lambda)$ was found in IR phase of pure-glue QCD Alexandru and Horváth
(2019). Indeed, the function is non-analytic at both
$\lambda_{\scriptscriptstyle{\text{\rm IR}}}$ and $\lambda_{A}$, with spectral
region of low-$d$ ($d_{\scriptscriptstyle{\text{\rm IR}}}\\!\leq\\!1$) modes
between them Alexandru and Horváth (2021). Moreover, in contrast to exact zero
modes, which are $d_{\scriptscriptstyle{\text{\rm IR}}}\\!=\\!3$, the lowest
near-zero modes $(\lambda\\!\to\\!0^{+})$ are close to other topological value
$d_{\scriptscriptstyle{\text{\rm IR}}}\\!=\\!2$. Such jump at
$\lambda_{\scriptscriptstyle{\text{\rm IR}}}\\!=\\!0$ is surprising since the
proposed origin of anomalous IR mode accumulation is the conventional mixing
of topological lumps Edwards  _et al._ (2000); Vig and Kovacs (2021) which, in
absence of additional (unknown) effects, leads to
$d_{\scriptscriptstyle{\text{\rm IR}}}\\!=\\!3$ in both cases. The jump could
thus offer valuable clues on IR phase dynamics, and could be used to detect
the transition into IR phase.

In this work, we make a key step toward this proposal becoming a reality: we
present evidence supporting the existence of the above unusual pattern also in
“real-world QCD”. In particular, in $N_{f}\\!=\\!2+1$ theory at physical quark
masses (see setup below) we obtain

$d_{\scriptscriptstyle{\text{\rm IR}}}(0)=2.98(09),\quad\quad
d_{\scriptscriptstyle{\text{\rm IR}}}(0^{+})=2.03(16)$ (3)

This finding lends support to topological origin of exotic IR-phase dynamics
Alexandru and Horváth (2021, 2022) (see also Edwards  _et al._ (2000); Vig and
Kovacs (2021); Cardinali  _et al._ (2021)), and significantly strenghtens its
connection to non-analyticities of Anderson-like origin Alexandru and Horváth
(2022). Regarding the latter, the observed topology aspects may have close
ties to critical Anderson dynamics Horváth and Markoš (2022a, 2023, b) which
is entirely unexpected.

2\. Simulation Setup: We lattice-regularize $N_{f}\\!=\\!2+1$ QCD using
tadpole-improved clover fermion action (1-step stout link smearing with
parameter 0.125) and tadpole-improved Symanzik gauge action at
$a=0.105~{}\mathrm{fm}$ and $m_{\pi}\simeq 135$ MeV. At temperature
$T\\!=\\!234\,$MeV, numerous spatial volumes (up to $L\\!=\\!10.1\,$fm) were
simulated by CLQCD collaboration (see Table 1), allowing for reliable
$d_{\scriptscriptstyle{\text{\rm IR}}}(\lambda)$ calculations. More detailed
ensemble description is given in ref . We note in passing that ensembles with
similar quark and gauge actions were already used in previous zero-temperature
calculations Zhang  _et al._ (2021); Liu  _et al._ (2022); Xing  _et al._
(2022); Liu  _et al._ (2023).

Table 1: UV cutoff $a$, pion mass $m_{\pi}$, lattice volumes $n_{L}^{3}\times n_{T}$ and temperature $T$ of lattice QCD ensembles studied. $a$(fm) |  $m_{\pi}$(MeV) | $n_{L}$ | $n_{T}$ |  $T$(MeV)
---|---|---|---|---
0.105 | 135 | 24/28/32/40/48/64/96 | 8 | 234

Glue fields $U$ of this theory will be studied via their effect on the overlap
Dirac operator $D_{\text{ov}}[U]$. We construct $D_{\text{ov}}$ using the
original square-root prescription Neuberger (1998) at $\rho\\!=\\!1.5$ with
1-step HYP smearing of $U$. To determine the low-lying eigensystem, we select
the chiral sector containing zero mode(s), calculate the eigenvectors of
$D_{\rm ov}^{\dagger}D_{\rm ov}$ in it using Arnoldi method, and then
construct non-zero modes Sorensen (1992); Lehoucq and Sorensen (1996); Giusti
_et al._ (2003); Alexandru  _et al._ (2011); ref . Transformation $D\equiv
D_{\textrm{ov}}/(1-\frac{a}{2\rho}D_{\textrm{ov}})$ Chiu and Zenkin (1999)
yields purely imaginary eigenvalues
($D\psi_{\lambda}(x)\\!=\\!i\lambda\psi_{\lambda}(x)$) and the associated
spectral density is
$\rho(\lambda)=T\sum_{i}\delta(\lambda-\lambda_{i})/L^{3}$. Further technical
details can be found in the supplemenatary material ref .

Eigenmodes with $\lambda$ up to $\sim 500\,$MeV were computed for all $L$ in
Table 1. Densities $\rho(\lambda)$ were renormalized in
$\overline{\textrm{MS}}$ at 2 GeV, using
$Z_{m}\\!=\\!Z^{-1}_{S}\\!=\\!0.907(26)$ obtained by interpolating the results
at 11 UV cutoffs He  _et al._ (2022). We wish to focus on a temperature in the
range $200\\!<\\!T\\!<\\!250$ MeV, where the system was originally predicted
to reach the IR phase at certain $T_{\scriptscriptstyle{\text{\rm IR}}}$. In
Fig. 1 we show $\rho(\lambda)$ at $T\\!=\\!234~{}$MeV. The striking bimodal
structure exhibits features previously associated with IR phase Alexandru and
Horváth (2019), including a fully-formed region of severe depletion. We also
include the $T\\!\simeq\\!0$ result from identical simulation setup on
$48^{3}\times 96$ lattice and $\rho(\lambda)$ obtained using stochastic
estimator Cossu  _et al._ (2016). The difference is indeed quite remarkable.
Note that, for large enough $\lambda$, the two densities come together which
is expected for all $T\\!\ll\\!1/a$.

Figure 1: Spectral density $\rho(\lambda)$ at $T\\!=\\!234\,$MeV (circles) and
$T\\!\simeq\\!0$ (triangles), both at $L=5.0$ fm.

3\. The Results: We now examine in detail whether the unusual
$d_{\scriptscriptstyle{\text{\rm IR}}}(\lambda)$ in IR phase of pure-glue QCD
Alexandru and Horváth (2021) is also present in real-world QCD at
$T\\!=\\!234~{}$MeV. To that end, we utilize and extend the techniques of
early studies. A useful concept is the “finite-volume”
$d_{\scriptscriptstyle{\text{\rm IR}}}$, namely Alexandru and Horváth (2021)

$d_{\scriptscriptstyle{\text{\rm
IR}}}(L,s)=\frac{1}{\ln(s)}\ln\frac{\mathscr{N}_{\star}(L)}{\mathscr{N}_{\star}(L/s)}\quad,\quad
s>0$ (4)

since then $d_{\scriptscriptstyle{\text{\rm
IR}}}\\!=\\!\lim_{L\to\infty}d_{\scriptscriptstyle{\text{\rm IR}}}(L,s)$
independently of $s$. Estimating the limit from linear extrapolations in $1/L$
work well in Anderson models, at least for extended states and at criticality
Horváth and Markoš (2022a). Here we utilize this, and point out that the
procedure is equivalent to direct fitting of $\mathscr{N}_{\star}(L)$ to the
form $b\,L^{d_{\scriptscriptstyle{\text{\rm IR}}}}e^{-c/L}$ ref , which is
technically more convenient. Using the data from our five largest systems in
such fits, we obtained $d_{\scriptscriptstyle{\text{\rm IR}}}(\lambda)$ shown
in Fig. 2. Despite some differences (see below), its behavior is strikingly
similar to pure-glue case (Fig. 1 of Ref. Alexandru and Horváth (2021)).

Important commonality is the discontinuity feature at
$\lambda_{\scriptscriptstyle{\text{\rm IR}}}\\!=\\!0$, suggesting that exact
zero-modes ($d_{\scriptscriptstyle{\text{\rm IR}}}(0)\\!\simeq\\!3$) differ
from lowest near-zero modes ($d_{\scriptscriptstyle{\text{\rm
IR}}}(0^{+})\\!\simeq\\!2$) in a robust qualitative manner. This is made
explicit by the inset of Fig. 2 focusing on very deep IR, and yielding
linearly $\lambda\\!\rightarrow\\!0^{+}$ extrapolated
$d_{\scriptscriptstyle{\text{\rm IR}}}(0^{+})\\!=\\!2.03(16)$, which is more
than $5\sigma$ smaller than 3. Explaining this difference in terms of the
underlying IR glue may prove important for deciphering the nature of IR phase.

Like in pure-glue case, we find a clear low-$d$
($d_{\scriptscriptstyle{\text{\rm IR}}}\leq 1$) plateau, here in the range of
about $10-220~{}$MeV. This roughly coincides with the region of strongly
suppressed $\rho(\lambda)$ (Fig. 1 vs Fig. 2). Dimensional structure of
plateaus will be further clarified using the multidimensionality technique of
Ref. Horváth and Markoš (2022b). Such studies are forthcoming.

The onset of rise toward dimension 3 past $\sim 230\,$MeV confirms the
viability of our scenario with mobility edges
$\lambda_{\scriptscriptstyle{\text{\rm IR}}}$ and $\pm\lambda_{A}$. However,
the discontinuity of $d_{\scriptscriptstyle{\text{\rm IR}}}(\lambda)$ at
presumed $\lambda_{A}$ is not apparent in our data, contrary to both the pure-
glue case Alexandru and Horváth (2021) and the situation in Anderson models
Horváth and Markoš (2022a). Resolution of this difference may provide an
additional new insight into the IR-phase dynamics.

Figure 2: Function $d_{\rm IR}(\lambda)$ at $T\\!=\\!234\,$MeV. The inset
zooms in on deep infrared with $\lambda\to 0^{+}$ extrapolation shown.

To illustrate the quality of scaling in various $\lambda$-regimes shown in
Fig. 2, we plot in Fig. 3 the fraction of volume occupied by the effective
support of the state, namely
$f_{\star}\\!\equiv\\!\mathscr{N}_{\star}/N\\!=\\!\mathscr{N}_{\star}/(n_{\scriptscriptstyle\text{L}}^{3}n_{\scriptscriptstyle\text{T}})$.
Since $\mathscr{N}_{\star}(L)\\!\propto\\!L^{d_{\scriptscriptstyle{\text{\rm
IR}}}}e^{-c/L}$ is used to extract $d_{\scriptscriptstyle{\text{\rm IR}}}$, we
have $f_{\star}(L)\\!\propto\\!L^{d_{\scriptscriptstyle{\text{\rm
IR}}}-3}e^{-c/L}$ and these fits are shown in Fig. 3. The displayed
$\chi^{2}$/dof for modes in different regimes do indeed confirm very good
scaling behavior. Note how functions $f_{\star}(L)$ in Fig. 3 visually
separate the bulk modes and near-bulk modes from IR modes. Indeed, although
zero-modes are $d_{\scriptscriptstyle{\text{\rm IR}}}\\!=\\!3$, and hence
occupy a finite fraction of volume in thermodynamic limit
($\lim_{L\to\infty}f_{\star}(L)\\!>\\!0$), its magnitude is much smaller than
that of typical bulk modes. At the same time, for
$d_{\scriptscriptstyle{\text{\rm IR}}}\\!<\\!3$ modes of IR component the
fraction vanishes ($\lim_{L\to\infty}f_{\star}(L)\\!=\\!0$).

Fig. 3 reveals that the lowest near-zero modes we studied
($\lambda\\!=\\!0.22\,$MeV, $d_{\scriptscriptstyle{\text{\rm IR}}}\\!<\\!3$)
have larger $f_{\star}$ at studied volumes than those of zero modes
($\lambda\\!=\\!0$, $d_{\scriptscriptstyle{\text{\rm IR}}}\\!=\\!3)$. But
given their $d_{\scriptscriptstyle{\text{\rm IR}}}$, this order has to reverse
at sufficiently large volume. We can read off the graphs in Fig. 3 that this
happens at $L\approx 20\,$fm. Such deep IR thresholds simply do not appear in
other QCD regimes. Only at larger $L$ will modes at $\lambda\\!=\\!0.22\,$MeV
become “sparser” than zero modes. Note that the qualitative difference between
zero and near-zero modes is expressed here by the opposite convexity
properties of their $f_{\star}(L;\lambda)$.

Figure 3: Function $f_{\star}(L)$ for various $\lambda$ at $T\\!=\\!234\,$MeV.
Figure 4: Typical color-coded $\log_{10}(\min\,\\{Np(x_{i}),1\\})$ in a 2d
plane containing $x_{i}$ with maximal probability. Modes in different
$\lambda$-regimes are shown. $T\\!=\\!234$ MeV and $L\\!=\\!10.1\,$fm.

Finally, we wish to gain some visual insight into the spatial geometry of
modes. In definition (1) of $\mathscr{N}_{\star}$, uniform probability
$p_{u}\\!=\\!1/N$ enters as a reference value: points $x_{i}$ with
$p(x_{i})\\!=\\!\psi_{\lambda}^{\dagger}(x_{i})\psi_{\lambda}(x_{i})\\!\geq\\!p_{u}$
are guaranteed to be in effective support, and we refer to them as “core”. We
wish to set up a sea-level representation that visualizes it sharply. Plotting
$\min\\{Np(x_{i}),1\\}$, namely the contribution of $x_{i}$ to effective
count, accomplishes that. In Fig. 4 we color-code this input (on a logarithmic
scale) and show its typical behavior on a plane containing the global
probability maximum. The black regions mark the core. The panels represent
different $\lambda$-regimes on the same glue background. The bulk mode at
$\lambda\\!=\\!330\,$MeV (right) resembles to some extent modes at low
temperatures. Indeed, its core spreads out contiguously over large distances
and its granularity (composition from distinct lumps) is not very obvious. To
the contrary, the plateau mode ($\lambda\\!=\\!100\,$MeV) is usually dominated
by a well-formed lump as shown. The near-zeromodes ($\lambda\\!=\\!0.22\,$MeV)
maintain the granularity, but involve multiple lumpy features forming a larger
spatial structure. The zero-modes (left) at this volume are in fact quite
similar but, due to $d_{\scriptscriptstyle{\text{\rm IR}}}\\!=\\!3$ vs
$d_{\scriptscriptstyle{\text{\rm IR}}}\\!=\\!2$ difference, will become
infinitely more “space-filling” in thermodynamic limit.

Additional results are described in our supplement ref .

4\. Summary and Discussion: Remarkable property of QCD IR phase Alexandru and
Horváth (2019) is that it requires the presence of non-analyticities not only
at the transition point $T_{\scriptscriptstyle{\text{\rm IR}}}$, but at any
temperature within the phase. It was proposed and verified Alexandru and
Horváth (2021, 2022) that in pure-glue QCD the system arranges for this by
reconfiguring itself into two independent components (IR & bulk), sharply
separated in Dirac spectrum. The needed non-analyticity enters via Anderson-
like mobility edges $\lambda_{\scriptscriptstyle{\text{\rm IR}}}\\!=\\!0$ and
$\pm\lambda_{\text{A}}\\!\neq\\!0$, encoded by dimension function
$d_{\scriptscriptstyle{\text{\rm IR}}}(\lambda)$. Our present results suggest
that key elements of this scenario also materialize in “real-world” QCD. Thus,
in certain regards, thermal state in IR phase of strong interactions resembles
the Tisza-Landau two-fluid model of liquid helium. The proposed 2-component
nature of thermal state may in fact be the most essential attribute of IR
phase.

We wish to point out certain aspects of our results. (i) The computed
dimension $d_{\scriptscriptstyle{\text{\rm IR}}}$ of near-zero modes is in
close vicinity of “topological” value 2, thus inviting a systematic inquiry
into its possible origin in certain topological feature of underlying glue
fields. At the same time, recent findings of possible topological behavior in
3d Anderson model Horváth and Markoš (2022b) also involve dimension $2$ but no
glue fields. (ii) In the existing QCD data there is no clear evidence yet for
critical value $d_{\scriptscriptstyle{\text{\rm IR}}}\\!\approx\\!8/3$, which
was suggested to be a generic feature of Anderson models Horváth and Markoš
(2022a). (iii) Unlike in the case of $\lambda_{\scriptscriptstyle{\text{\rm
IR}}}$, in the vicinity of $\lambda_{A}$ we did not find an obviuous dimension
jump. This differs from situation in pure-glue QCD Alexandru and Horváth
(2021) and from that at critical points of Anderson models Horváth and Markoš
(2022a). Taken together, points (i-iii) constitute an intriguing complex
puzzle to be solved by future studies. Indeed, the satisfactory understanding
of Anderson-like features in QCD require the resolution of these issues. This
resolution will involve control over the usual lattice QCD systematics, which
is important given that some of these effects are possibly enhanced in
dynamical simulations Zhao  _et al._ (2023).

Recently a number of lattice QCD papers focused on the same temperature range
as the one investigated here and in other recent IR phase studies (see e.g.
Dick  _et al._ (2015); Ding  _et al._ (2020); Aoki  _et al._ (2021); Kaczmarek
_et al._ (2021); Kotov  _et al._ (2021); Chen  _et al._ (2022); Kaczmarek  _et
al._ (2023)). Their physics goals are mostly different and tend to involve the
chiral limit, such as in studies of U${}_{\text{A}}$(1) problem or chiral
phase transition. Other related developments include an approximate color-spin
symmetry Rohrhofer  _et al._ (2019); Glozman  _et al._ (2022), as well as
recent Refs. Liu (2023); Athenodorou  _et al._ (2022); Kehr  _et al._ (2023).
Conversely, the present CLQCD data could be used to study these other
problems.

## Acknowledgments

We thank the CLQCD collaborations for providing their gauge configurations,
and also to Ting-Wai Chiu, Heng-Tong Ding and Jian Liang for valuable
discussions. This work is supported in part by the Strategic Priority Research
Program of Chinese Academy of Sciences, Grant No. XDB34030303 and XDPB15, NSFC
grants No. 12293062, 12293065 and 12047503, and also a NSFC-DFG joint grant
under Grant No. 12061131006 and SCHA 458/22. This work is also supported in
part by the U.S. DOE Grant No. DE-SC0013065, DOE Grant No. DE-FG02-95ER40907,
and DOE Grant No. DE-AC05-06OR23177 which is within the framework of the TMD
Topical Collaboration. This publication received funding from the French
National Research Agency under the contract ANR-20-CE31-0016. The numerical
calculation were carried out on the ORISE Supercomputer through HIP
programming model Bi  _et al._ (2020), and HPC Cluster of ITP-CAS. This
research used resources of the Oak Ridge Leadership Computing Facility at the
Oak Ridge National Laboratory, which is supported by the Office of Science of
the U.S. Department of Energy under Contract No. DE-AC05-00OR22725. This work
used Stampede time under the Extreme Science and Engineering Discovery
Environment (XSEDE), which is supported by National Science Foundation Grant
No. ACI-1053575. We also thank the National Energy Research Scientific
Computing Center (NERSC) for providing HPC resources that have contributed to
the research results reported within this paper. We acknowledge the facilities
of the USQCD Collaboration used for this research in part, which are funded by
the Office of Science of the U.S. Department of Energy.

## References

  * Hagedorn (1965) R. Hagedorn, Nuovo Cim. Suppl. 3, 147 (1965).
  * Hagedorn (1971) R. Hagedorn, CERN report 71-12 (1971), 10.5170/CERN-1971-012.
  * McLerran and Svetitsky (1981) L. D. McLerran and B. Svetitsky, Phys. Lett. B 98, 195 (1981).
  * Kuti  _et al._ (1981) J. Kuti, J. Polonyi, and K. Szlachanyi, Phys. Lett. B 98, 199 (1981).
  * Engels  _et al._ (1981) J. Engels, F. Karsch, H. Satz, and I. Montvay, Phys. Lett. B 101, 89 (1981).
  * Aoki  _et al._ (2006) Y. Aoki, G. Endrodi, Z. Fodor, S. Katz, and K. Szabo, Nature 443, 675 (2006), arXiv:hep-lat/0611014 [hep-lat] .
  * Arsene  _et al._ (2005) I. Arsene _et al._ (BRAHMS), Nucl. Phys. A757, 1 (2005), arXiv:nucl-ex/0410020 [nucl-ex] .
  * Back  _et al._ (2005) B. B. Back _et al._ , Nucl. Phys. A757, 28 (2005).
  * Adams  _et al._ (2005) J. Adams _et al._ (STAR), Nucl. Phys. A757, 102 (2005), arXiv:nucl-ex/0501009 [nucl-ex] .
  * Adcox  _et al._ (2005) K. Adcox _et al._ (PHENIX), Nucl. Phys. A757, 184 (2005), arXiv:nucl-ex/0410003 [nucl-ex] .
  * Muller  _et al._ (2012) B. Muller, J. Schukraft, and B. Wyslouch, Ann. Rev. Nucl. Part. Sci. 62, 361 (2012), arXiv:1202.3233 [hep-ex] .
  * Alexandru and Horváth (2019) A. Alexandru and I. Horváth, Phys. Rev. D 100, 094507 (2019), arXiv:1906.08047 [hep-lat] .
  * Edwards  _et al._ (2000) R. G. Edwards, U. M. Heller, J. E. Kiskis, and R. Narayanan, Phys.Rev. D61, 074504 (2000), arXiv:hep-lat/9910041 [hep-lat] .
  * Alexandru and Horváth (2015) A. Alexandru and I. Horváth, Phys. Rev. D92, 045038 (2015), arXiv:1502.07732 [hep-lat] .
  * Alexandru and Horváth (2021) A. Alexandru and I. Horváth, Phys. Rev. Lett. 127, 052303 (2021), arXiv:2103.05607 [hep-lat] .
  * Alexandru and Horváth (2022) A. Alexandru and I. Horváth, Phys. Lett. B 833, 137370 (2022), arXiv:2110.04833 [hep-lat] .
  * Garcia-Garcia and Osborn (2006) A. M. Garcia-Garcia and J. C. Osborn, Nucl. Phys. A 770, 141 (2006), arXiv:hep-lat/0512025 .
  * Garcia-Garcia and Osborn (2007) A. M. Garcia-Garcia and J. C. Osborn, Phys. Rev. D 75, 034503 (2007), arXiv:hep-lat/0611019 .
  * Kovacs and Pittler (2010) T. G. Kovacs and F. Pittler, Phys. Rev. Lett. 105, 192001 (2010), arXiv:1006.1205 [hep-lat] .
  * Giordano  _et al._ (2014) M. Giordano, T. G. Kovacs, and F. Pittler, Phys. Rev. Lett. 112, 102002 (2014), arXiv:1312.1179 [hep-lat] .
  * Ujfalusi  _et al._ (2015) L. Ujfalusi, M. Giordano, F. Pittler, T. G. Kovács, and I. Varga, Phys. Rev. D 92, 094513 (2015), arXiv:1507.02162 [cond-mat.dis-nn] .
  * Horváth  _et al._ (2023) I. Horváth, P. Markoš, and R. Mendris, Entropy 25, 482 (2023), arXiv:2205.11520 [hep-lat] .
  * Horváth and Mendris (2020) I. Horváth and R. Mendris, Entropy 22, 1273 (2020), arXiv:1807.03995 [quant-ph] .
  * Horváth (2021) I. Horváth, Quantum Rep. 3, 534 (2021), arXiv:1809.07249 [quant-ph] .
  * Horváth and Mendris (2019) I. Horváth and R. Mendris, MDPI Proc. 13, 8 (2019), arXiv:1907.01606 [quant-ph] .
  * Neuberger (1998) H. Neuberger, Phys. Lett. B417, 141 (1998), arXiv:hep-lat/9707022 [hep-lat] .
  * Vig and Kovacs (2021) R. A. Vig and T. G. Kovacs, Phys. Rev. D 103, 114510 (2021), arXiv:2101.01498 [hep-lat] .
  * Cardinali  _et al._ (2021) M. Cardinali, M. D’Elia, and A. Pasqui, (2021), arXiv:2107.02745 [hep-lat] .
  * Horváth and Markoš (2022a) I. Horváth and P. Markoš, Phys. Rev. Lett. 129, 106601 (2022a), arXiv:2110.11266 [cond-mat.dis-nn] .
  * Horváth and Markoš (2023) I. Horváth and P. Markoš, Phys. Lett. A 467, 128735 (2023), arXiv:2207.13569 [cond-mat.dis-nn] .
  * Horváth and Markoš (2022b) I. Horváth and P. Markoš, (2022b), arXiv:2212.09806 [cond-mat.dis-nn] .
  * (32) Supplementary materials .
  * Zhang  _et al._ (2021) Q.-A. Zhang, J. Hua, F. Huang, R. Li, Y. Li, C.-D. Lu, P. Sun, W. Sun, W. Wang, and Y.-B. Yang, (2021), arXiv:2103.07064 [hep-lat] .
  * Liu  _et al._ (2022) H. Liu, J. He, L. Liu, P. Sun, W. Wang, Y.-B. Yang, and Q.-A. Zhang, (2022), arXiv:2207.00183 [hep-lat] .
  * Xing  _et al._ (2022) H. Xing, J. Liang, L. Liu, P. Sun, and Y.-B. Yang, (2022), arXiv:2210.08555 [hep-lat] .
  * Liu  _et al._ (2023) H. Liu, L. Liu, P. Sun, W. Sun, J.-X. Tan, W. Wang, Y.-B. Yang, and Q.-A. Zhang, (2023), arXiv:2303.17865 [hep-lat] .
  * Sorensen (1992) D. C. Sorensen, SIAM J. on Matrix Analysis and Applications 13, 357 (1992).
  * Lehoucq and Sorensen (1996) R. B. Lehoucq and D. C. Sorensen, SIAM J. Matrix Anal. Appl. 17, 789 (1996).
  * Giusti  _et al._ (2003) L. Giusti, C. Hoelbling, M. Luscher, and H. Wittig, Comput. Phys. Commun. 153, 31 (2003), arXiv:hep-lat/0212012 [hep-lat] .
  * Alexandru  _et al._ (2011) A. Alexandru, M. Lujan, C. Pelissier, B. Gamari, and F. X. Lee, in  _Proceedings, SAAHPC’11_ (2011) pp. 123–130, arXiv:1106.4964 [hep-lat] .
  * Chiu and Zenkin (1999) T.-W. Chiu and S. V. Zenkin, Phys. Rev. D59, 074501 (1999), arXiv:hep-lat/9806019 [hep-lat] .
  * He  _et al._ (2022) F. He, Y.-J. Bi, T. Draper, K.-F. Liu, Z. Liu, and Y.-B. Yang ($\chi$QCD), Phys. Rev. D 106, 114506 (2022), arXiv:2204.09246 [hep-lat] .
  * Cossu  _et al._ (2016) G. Cossu, H. Fukaya, S. Hashimoto, T. Kaneko, and J.-I. Noaki, PTEP 2016, 093B06 (2016), arXiv:1607.01099 [hep-lat] .
  * Zhao  _et al._ (2023) D.-J. Zhao, G. Wang, F. He, L. Jin, P. Sun, Y.-B. Yang, and K. Zhang ($\chi$QCD), Phys. Rev. D 107, L091501 (2023), arXiv:2207.14132 [hep-lat] .
  * Dick  _et al._ (2015) V. Dick, F. Karsch, E. Laermann, S. Mukherjee, and S. Sharma, Phys. Rev. D91, 094504 (2015).
  * Ding  _et al._ (2020) H. T. Ding, S. T. Li, S. Mukherjee, A. Tomiya, X. D. Wang, and Y. Zhang, (2020), arXiv:2010.14836 [hep-lat] .
  * Aoki  _et al._ (2021) S. Aoki, Y. Aoki, G. Cossu, H. Fukaya, S. Hashimoto, T. Kaneko, C. Rohrhofer, and K. Suzuki (JLQCD), Phys. Rev. D 103, 074506 (2021), arXiv:2011.01499 [hep-lat] .
  * Kaczmarek  _et al._ (2021) O. Kaczmarek, L. Mazur, and S. Sharma, Phys. Rev. D 104, 094518 (2021).
  * Kotov  _et al._ (2021) A. Y. Kotov, M. P. Lombardo, and A. Trunin, Physics Letters B 823, 136749 (2021).
  * Chen  _et al._ (2022) Y.-C. Chen, T.-W. Chiu, and T.-H. Hsieh (TWQCD), Phys. Rev. D 106, 074501 (2022), arXiv:2204.01556 [hep-lat] .
  * Kaczmarek  _et al._ (2023) O. Kaczmarek, R. Shanker, and S. Sharma, (2023), arXiv:2301.11610 [hep-lat] .
  * Rohrhofer  _et al._ (2019) C. Rohrhofer, Y. Aoki, G. Cossu, H. Fukaya, C. Gattringer, L. Ya. Glozman, S. Hashimoto, C. B. Lang, and S. Prelovsek, Phys. Rev. D100, 014502 (2019).
  * Glozman  _et al._ (2022) L. Y. Glozman, O. Philipsen, and R. D. Pisarski, Eur. Phys. J. A 58, 247 (2022), arXiv:2204.05083 [hep-ph] .
  * Liu (2023) K.-F. Liu, (2023), arXiv:2302.11600 [hep-ph] .
  * Athenodorou  _et al._ (2022) A. Athenodorou, C. Bonanno, C. Bonati, G. Clemente, F. D’Angelo, M. D’Elia, L. Maio, G. Martinelli, F. Sanfilippo, and A. Todaro, JHEP 10, 197 (2022), arXiv:2208.08921 [hep-lat] .
  * Kehr  _et al._ (2023) R. Kehr, D. Smith, and L. von Smekal, (2023), arXiv:2304.13617 [hep-lat] .
  * Bi  _et al._ (2020) Y.-J. Bi, Y. Xiao, W.-Y. Guo, M. Gong, P. Sun, S. Xu, and Y.-B. Yang,  _Proceedings, Lattice 2019_ , PoS LATTICE2019, 286 (2020), arXiv:2001.05706 [hep-lat] .

##  I Supplemental Materials

###  I.1 Visualization of Spatial Distributions

In the present subsection we wish to extend our visualization of mode
distributions in various $\lambda$-regimes shown in Fig. 4 of the manuscript.
In particular, Fig. 5 shows examples of modes at the same four values of
$\lambda$, but at all IR cutoffs (sizes $L$ of the system) considered in this
work. Note that $\lambda\\!=\\!\,$0, 0.22 and 100 MeV belong to the IR part of
the spectrum, while $\lambda\\!=\\!\,$330 MeV is in near-bottom part of the
bulk.

Interesting aspect of observing the typical geometry at fixed $\lambda$ for
increasing $L$ is an evolution in degree of granularity. Indeed, note that for
the plateau mode ($\lambda\\!=\\!100\,$MeV), increasing $L$ confirms the
picture of a single solid lump present in the volume. On the other hand, for
zero-modes and near-zero modes, the apparent (visually observable) degree of
granularity increases with increasing $L$, reflecting that their effective
supports keeps spreading out in the volume. In fact, all qualitative aspects
we observe agree with metal-to-critical picture of transition to IR phase, put
forward in Ref. Alexandru and Horváth (2022). The associated details will be
worked out in dedicated follow-up publications.

Figure 5:  Typical color-coded $\log_{10}(\min\,\\{Np(x_{i}),1\\})$ in a 2d
plane containing point $x_{i}$ with maximal probability. Modes in different
$\lambda$-regimes on given equilibrium glue background at $T\\!=\\!234$ MeV
are shown for all spatial sizes $L$ studied.

###  I.2 Accuracy of Overlap Eigenmodes

In this section we focus on the accuracy of the low-lying eigenvectors used in
this study. For efficiency reasons, we compute the low-lying eigenvalues and
eigenvectors of $M\equiv\frac{1}{\rho^{2}}D_{\rm ov}^{\dagger}D_{\rm
ov}=\frac{1}{\rho}(D_{\rm ov}+D_{\rm ov}^{\dagger})$ in the chiral sector
which includes the exact zero mode. The eigenvalues and eigenvectors of $D$
are simply related to the eigenvectors of $M$: non-zero modes of
$D_{\text{ov}}$ come in pairs, $\lambda_{R}\pm i\lambda_{I}$, and they span an
invariant two-dimensional space for $M$ with $\lambda_{M}=\lambda_{R}$ and for
near zero-modes $\lambda_{R}\approx\lambda_{I}^{2}/(2\rho)$. When computing
the spectrum of $M$, zero modes appear as Ritz pairs with eigenvalues of the
order of $\epsilon$, the precision of the sign-function approximation used to
implement $D_{\text{ov}}$. When we have near-zero modes with
$(a\lambda_{I})^{2}<\epsilon$, it is impossible to distinguish them from zero
modes. Using a polynomial approximation for the sign functions, the best
precision we are able to achieve is $\epsilon\approx 10^{-12}$, and
consequently we can only confidently resolve eigenvalues with
$a\lambda_{I}>10^{-6}$, which in physical units correspond to
$\lambda_{I}>2\,\text{eV}$. For the volumes used in this study, the near-zero
eigenmodes satisfy this condition.

Another concern is the mixing between nearly-degenerate eigenvectors. For
eigenvector observables (like $f_{\star}$) that are smooth as a function of
$\lambda$, this is less of a concern. However, at discontinuities mixing could
introduce systematic effects. This could potentially be a problem at
$\lambda=0$ since the zero modes and near-zero modes behave differently. We
argue here that this is not the case.

To see this, consider two eigenvectors of the projected operator $D$ with
$Dv_{1}=i\lambda_{1}v_{1}$ and $Dv_{2}=i\lambda_{2}v_{2}$. A mixed vector
$v=v_{1}\cos\theta+v_{2}\sin\theta$ has a Ritz “eigenvalue”
$i\lambda=v^{\dagger}Dv$. The residue for this vector is $\delta=\|Dv-i\lambda
v\|$ is going to be $|\cos\theta\sin\theta(\lambda_{1}-\lambda_{2})|$. For our
case we take $\lambda_{1}\approx 0$ and the near-zero have
$\lambda_{2}>0.1\,\text{MeV}$ and our residues, even in the worst case, are
$\delta>10^{-7}$. This implies that the mixing angle is at most
$\theta\sim\delta/\lambda_{2}=2\times 10^{-3}$. Given that $f_{\star}$ varies
slowly (the difference between zero modes and near-zero modes is less than
two), this mixing will have negligible effect given our statistical errors.

###  I.3 Fitting of Effective Fractions

Our procedure to extract $d_{\scriptscriptstyle{\text{\rm IR}}}$ assumes an
approximately linear (in $1/L$) approach of “finite-volume” dimension
$d_{\scriptscriptstyle{\text{\rm IR}}}(L,s)$ in Eq. (4) to its
$L\\!\to\\!\infty$ limit. This was suggested by Ref. Horváth and Markoš
(2022a) in the context of Anderson models. One can easily check that this is
equivalent to direct fits of $\mathscr{N}_{\star}(L)$ to the form
$\mathscr{N}_{\star}(L)\\!\propto\\!L^{d_{\scriptscriptstyle{\text{\rm
IR}}}}e^{-c/L}$ and
$f_{\star}(L)\\!\propto\\!L^{d_{\scriptscriptstyle{\text{\rm IR}}}-3}e^{-c/L}$
for effective volume fraction
$f_{\star}(L)\\!\equiv\\!\mathscr{N}_{\star}(L)/N$. Fits shown in the left
panel of Fig. 6 support the validity of our approach.

Figure 6: Function $f_{\star}(L)$ (left) and the associated
$\bar{d}_{\scriptscriptstyle{\text{\rm IR}}}(\bar{L})$ (right) for various
$\lambda$ at $T\\!=\\!234\,$MeV. See text for explanations.

Here we wish to check the nature of finite-$L$ correction more directly. Given
the above scaling form, $L$-dependence of IR dimension can be expressed as

$d_{\scriptscriptstyle{\text{\rm IR}}}(L,s)=d_{\scriptscriptstyle{\text{\rm
IR}}}+\frac{s-1}{\ln(s)}\,\frac{c}{L}\quad\;\;\longrightarrow\quad\;\;\bar{d}_{\scriptscriptstyle{\text{\rm
IR}}}(\bar{L})=d_{\scriptscriptstyle{\text{\rm
IR}}}+\frac{c}{\bar{L}}\quad\,\text{with}\quad\,\bar{d}_{\scriptscriptstyle{\text{\rm
IR}}}(\bar{L}(L,s))\equiv d_{\scriptscriptstyle{\text{\rm
IR}}}(L,s)\;\;,\;\;\bar{L}(L,s)\equiv L\,\frac{\ln(s)}{s-1}\;\;$ (5)

Introduction of variable $\bar{L}$ thus makes it possible to combine
$d_{\scriptscriptstyle{\text{\rm IR}}}$ results from all pairs of distinct
lattices and follow their trends. Indeed, according to the above, value
$\bar{d}_{\scriptscriptstyle{\text{\rm IR}}}(1/\bar{L})$ from each pair should
fall on an indicated straight line, at least near $1/\bar{L}\\!=\\!0$. To
check this, we show in the left plot of Fig. 6 our $f_{\star}$ data for five
largest volumes at selected values of $\lambda$, and in the right plot the
associated functions $\bar{d}_{\scriptscriptstyle{\text{\rm IR}}}(1/\bar{L})$.
Note that there are 10 data points for each $\lambda$ in the latter case since
this is the number of possible lattice pairs. Displayed fits are indeed
consistent with linear nature of $d_{\scriptscriptstyle{\text{\rm
IR}}}(1/\bar{L})$ near $1/\bar{L}\\!=\\!0$. The qualitative difference between
exact zero modes and lowest near-zero modes is expressed in the right plot by
the crossed lines representing $\lambda\\!=\\!0$ and
$\lambda\\!=\\!0.22\,$MeV. Their finite-volume corrections are in fact of
opposite sign.

Figure 7: Function $f_{\star}(L)$ with different simulated size $L$ at $T$=234
MeV, for $\lambda\in[0.1,2.0]$ MeV (left), $\lambda\in[20,200]$ MeV (middle)
and $\lambda\in[280,330]$ MeV (right).

We also provide the $f_{\star}(1/L;\lambda)$ data for more values of $\lambda$
in Fig. 7. As shown in the left panel, the
$f_{\star}(L>4~{}\mathrm{fm};0<\lambda<2~{}\mathrm{MeV})$ becomes higher when
$\lambda$ is smaller, and thus the corresponding $d_{\rm IR}$ is also larger.
The tendency converges at $\lambda\sim$ 0.2 MeV (orange band), which
corresponds to $d_{\rm IR}=2$ and is consistent with the $\lambda=$0.13 MeV
case (dark red band) within the uncertainty. However, this limit is
significantly different from $f_{\star}(;0)$ as illustrated by the red band.
Thus, $f_{\star}$ and also $d_{\rm IR}$ will be discontinued at $\lambda=0$.

In contrast, the middle panel of Fig. 7 shows that $d_{\rm IR}$ changes
smoothly with large $\lambda$ for $\lambda\in[20,200]$ MeV, corresponding to
$d_{\rm IR}\sim 1$ within 2$\sigma$. The changes on $d_{\rm IR}$ is also
smooth in the range of $\lambda\in[280,330]$ MeV where $d_{\rm IR}$ converges
to 3 with increasing $\lambda$, as in the right panel. Therefore the
discontinuation of $d_{\rm IR}$ would only occurs at $\lambda=0$, given our
statistical precision at $a=0.105$ fm and $T$=234 MeV.

Figure 8: Function $f_{\star}(L)$ at $T\\!=\\!234\,$MeV in different spectral
regions, for all simulated sizes $L$. Left: zero modes and near-zero modes.
Right: plateau, just below $\lambda_{A}$ and at the bottom of the bulk. Shaded
regions are excluded from the displayed fits.

Finally, we give a justification for using the five largest lattices in our
fits i.e. systems with $L\\!>\\!3\,$fm. To that end, we show in Fig. 8
functions $f_{\star}(1/L)$ for all simulated volumes, together with previously
shown fits. Shaded areas mark the volumes excluded from these fits. One can
see that in case of zero modes and near-zero modes (left plot), the systems in
shaded region do not follow the fit curves, and were thus excluded from fits
in all spectral regions.


# To Rain or Not to Rain: Correlating GOES Flare Class and Coronal Rain
Statistics

E I Mason Predictive Science Inc.
9990 Mesa Rim Rd, Suite 170
San Diego, CA 92121, USA  K Kniezewski United States Naval Academy
121 Blake Road
Annapolis, MD 21402, USA

(Received March 8, 2022; Accepted September 21, 2022)

###### Abstract

Post-flare arcades are well-known components of solar flare evolution, which
have been observed for several decades. Coronal rain, cascades of
catastrophically-cooled plasma, outline the loops and provide eye-catching
evidence of the recent flare. These events are acknowledged to be common, but
the scientific literature does not include any statistical overview
documenting just how common the phenomenon actually is. This study reviews
Solar Dynamics Observatory Atmospheric Imaging Assembly (SDO AIA) observations
of 241 flares collected from the Space Weather Prediction Center (SWPC)
database between 2011 and 2018. The flares cover the entire strength range of
the C, M, and X GOES classes, and are distributed evenly across the SDO-
observed majority of Solar Cycle 24. We find that post-flare arcade rain
occurs for nearly all X and most M-class flares, but that it tapers off
rapidly within C-class flares. There appears to be a cut-off point around C5,
below which the occurrence of post-flare arcade rain drops significantly.
There is also a general positive correlation between GOES class and the
average duration of post-flare rain events. Post-flare arcade rain events in
X- and M-class flares appear to track with the sunspot number, providing a
potential new tool for estimating, if not predicting, solar cycle strength.
Furthermore, arcades are observed to persist for up to several days after the
originating flare, transitioning from hosting post-flare rain to typical
quiescent active region condensations. These results open up further avenues
for future research, including new methods to estimate energy deposition and
to gain greater insight into steady active region heating.

††journal: ApJ

##  1 Introduction

This paper addresses coronal rain events in solar flares, investigating the
frequency with which such events occur, their duration, and other salient
details. We correlate these statistics to overarching open questions
concerning the Sun: flare evolution, coronal mass ejection (CME) creation,
solar cycle characteristics, and localized coronal heating.

White-light study of solar flares began in 1859 with the independent but
concurrent observations of Carrington (1859) and Hodgson (1859), and
investigation of these explosive events spread to both longer and shorter
wavelengths – sometimes intentionally, sometimes serendipitously – as
technology advanced. The current model for flares, known as the standard model
and succinctly summarized by Holman (2016), holds that magnetic reconnection
across field lines arrayed along a polarity inversion line releases energy,
heating the plasma to temperatures in excess of 30 MK and accelerating
particles along the newly-rearranged magnetic field. Flares are often cited as
sources of solar energetic particles (e.g., Reames (1990), Cliver et al.
(2012), Reames (2017)), radio bursts (for starters, see Wild & Smerd (1972),
Aschwanden et al. (1985), Nelson & Melrose (1985), and Cane et al. (2002)),
and CMEs (Feynman & Hundhausen (1994), Švestka (2001), Andrews (2003), Yashiro
et al. (2005), Cane & Lario (2006), among many others).

According to the flare Standard Model (Holman (2016)), shearing along a
polarity inversion line (PIL) elongates field lines; eventually, some form of
disturbance can cause these lines to distort, triggering reconnection along
the PIL. This builds up a flux rope (which carries twist, translated upward
into the corona from the shearing) above an arcade of loops that are arrayed
roughly perpendicular with respect to the PIL. If the flare is eruptive, the
flux rope is then ejected into the outer corona and the solar wind as a CME
(Lin & Forbes (2000), Mittal & Narain (2010)); otherwise, it is considered a
confined flare. The conditions which determine whether or not a flare is
confined are not well-known (Sun et al. (2015), Zuccarello et al. (2017),
Mason et al. (2021)). These events carry magnetic field as well as plasma into
the heliosphere, and can strongly affect systems both on the Earth and in
space. Early observations tied their existence almost solely to solar flares,
but it was subsequently found that CMEs also frequently form from more
quiescent events like prominence eruptions (Gosling (1993)). However, strong
(and particularly long-duration) flares are well-correlated with CMEs, as
discussed in Sheeley et al. (1983) and Yashiro et al. (2005).

The arcade loops that form below the flux rope are filled with extremely hot,
dense plasma through electron beam heating, in a process termed ”chromospheric
evaporation”. This results in a loop which is severely outside the normal
coronal conditions, resulting in thermal instability (Parker (1953), Priest
(1978), Antiochos (1980), and Klimchuk (2019), among many others). The
temperature and density in the loop are unsupported by the energy being added
by whatever ambient background or footpoint heating is occurring in the arcade
loop, so the loop begins radiating rapidly. Here we reach the end of what is
well-understood about flare arcades: if the temperature profile of the loop
was uniform, we would expect the arcade loop to cool and evacuate in its
entirety, as was found by Reep et al. (2020) for loops heated purely by the
electron beams. Contrary to this, the recent work of Ruan et al. (2021)
successfully reproduced post-flare rain with a 2.5D MHD simulation without
even including beams. That model gets sufficient heating and condensations
simply from the magnetic energy conversion during flare reconnection.
Observations show that individual condensations occur within arcade loops;
this implies instead that a localized temperature inversion is critical to the
post-flare arcade rain formation mechanism, to seed nascent rain blobs.

In this work, we focus on these diminutive condensations. It is important to
point out and clearly define the phenomenon which we study here, because it
delineates an important distinction for observations. Any typical active
region will host significant numbers of coronal condensations on its loops
during its lifetime, and these may increase in number in the time after the
impulsive phase of a flare. However, since our goal in this paper is to
quantify that rain which is an immediate effect of flare reconnection, the
most natural observational selection was to only consider condensations formed
within the flare arcade. One of the main findings of this study, however, was
the discovery that _the arcades themselves can often become simply other
coronal loops that host quiescent rain, given enough time_ : this frequently
occurs with medium-sized flares, where the arcade grows to a point but does
not vanish. Since the rain which can occur continuously in these loops several
days after the flare is almost certainly not causally linked to the flare,
this rain also cannot be considered. Therefore, in our observations we also
needed a marker for the end of the rain duration. We define it to be at the
point that the raining arcade loops ceased to increase in height (best for
off-limb flares), or when the flare ribbons stopped evolving (best for on-disk
flares). Both of these measures – which are two observables of the same
physical process – signal that the flare reconnection has ceased and that
further rain on the same loops is likely due to another heating source. Using
this narrow definition, we consider only the rain directly tied to flare
reconnection, producing results that are the most relevant to flare dynamics.
This is understood by our use of the common terms _post-flare rain_ or _post-
flare arcade rain_ throughout this paper.

One of the greatest benefits of the Solar Dynamics Observatory Atmospheric
Imaging Assembly (SDO AIA; Pesnell et al., 2011; Lemen et al., 2011) has been
the sheer number of events which have been captured in detail by the near-
constant view across a broad range of extreme ultraviolet (EUV) passbands.
Much of the flare literature revolves around in-depth case studies of
individual flares, or small collections of similar flares. Here we conduct an
assessment of post-flare rain in a large pool of events. Rain has been
reported for years as a feature of strong flares, but the frequency and
statistical properties of it has not, to our knowledge, ever been quantified.
While not the brightest signal of the elusive magnetic dynamics involved in
flares, post-flare rain provides a critical reservoir which we can use to
measure how much of the flare’s energy is actually deposited in the plasma
still trapped in the corona. Energy which is radiated away as the
condensations form is energy that is no longer contained in the magnetic
field, powering a CME or an SEP event (Reames (2013), Drake & Swisdak (2012)).
We review all of the X-class flares that occurred between 2011 and 2019, as
well as a large, representative sample of M- and C-class flares from the same
period (we detail the collection and analysis of these in Section 2). This
survey allows us to build a strong basis for general statements about energy
budgeting and magnetic field dynamics in more energetic flares.

##  2 Data Collection and Analysis

###  2.1 Data Selection

We queried the Heliophysics Events Knowledgebase (HEK; accessible at
https://www.lmsal.com/hek/index.html) for flares detected by the GOES mission
between 2011 January 1 and 2019 December 31. We selected this date range to
cover most of Solar Cycle 24 for which we have routine daily SDO coverage.
Solar Cycle 24 began in December 2008, but SDO finished commissioning and
began its main science mission in May 2010. During this time period, there
were 49 X-class flares, 732 M-class flares, and 7681 C-class flares. Since
this study was conducted manually, a subset of the M and C-class flares were
selected; all of the X-class flares were analyzed. We chose flares to create a
pool that covered roughly equal numbers of flares in all relevant
characteristics: class (i.e., M), intra-class magnitude (i.e., M5), and
distribution across the solar cycle (equal numbers from each year covered by
the study, when possible). All the flares in the data set were selected for
their adherence to these characteristics, and prior to analysis; a flare was
not viewed until it had been selected for the study. If enough flares in a
category were unable to be analyzed due to data dropouts, eclipsing of SDO, or
similar reasons such that it would skew the distribution of the final data
set, more flares were chosen from the list.

The final selection included all 49 X-class flares, as well as 78 M-class and
118 C-class flares. 5 of the X-class flares had data dropouts which precluded
analysis, so the statistics here are on the 44 which had available AIA data. A
few other events had truncated rain durations because of data dropouts or
subsequent flares in the same location; these are all noted as such in the
full dataset, available on Zenodo at https://doi.org/10.5281/zenodo.7028211.
From the GOES flare list, we obtained the flare start, peak, and end times;
the location on the Sun in heliographic coordinates; and the source NOAA
active region designation.

###  2.2 Observational Analysis and Challenges

The analysis itself was primarily conducted via visual inspection of SDO AIA
data. Movies were generated with JHelioviewer (Müller et al. (2017)), showing
eight hours of data at a time beginning with the start time of the flare, and
171 and 304 Å data were evaluated for signs of cool condensations (see Figure
1 for a still-frame of one such movie) within flare arcades. Observationally,
arcade loops appear as relatively short, bright loops with footpoints rooted
at (indeed, defining) the flare ribbons. They generally form perpendicular to
the solar surface, and are almost always symmetric about their apex. These
were the observational hallmarks we used, which include characteristics that
make them readily identifiable both on-disk and off-limb.

When post-flare rain was present up until the end of the 8-hour movie, the
inspection time period was extended (this occurred primarily for X-class
flares, as can be seen from Figure 2). 304 Å is the coolest main SDO EUV
channel, and the one most commonly used to observe coronal rain in emission;
we chose 171 Å for its cool, well-defined loop structures in active regions,
where condensations can often appear clearly outlined as both emission and
absorption features. Hotter lines also sometimes show rain in absorption, but
the attendant diffuse glow from the active region obfuscates individual loops,
and so these were not widely used in our analysis. Once post-flare arcade rain
was established to be observed in 304, the duration of the event was observed
and recorded, in accordance with the definition we outlined in Section 1.

There are, of course, pitfalls to attempting to observe small, cool structures
in 304, particularly when the flare occurs on disk. One of these challenges is
that 304 frequently saturates early in a flare, and some detail is lost. As
can be seen from Figure 1, we adjusted the maximum and minimum values for the
color table from the default scaling that is commonly applied to AIA data;
while this did not entirely eliminate the early saturation during the
impulsive phase, it did frequently resolve structure significantly earlier
than would normally be observed with the standard scaling rules. The early
saturation is a well-known problem which cannot be overcome without
introducing another instrument, which would have significantly restricted our
pool of events.

Another common observational challenge relates to identifying cool
condensations when the flare occurred on-disk. Since 304 Å observations are
dominated on-disk by chromospheric and transition region emission, rain can be
difficult to pick out. However, this problem affects quiescent rain more than
post-flare arcade rain; when observing the flares, the activated source region
was bright enough that the cool blobs were easily viewed the majority of the
time. For cases where the observation was somewhat ambiguous, we searched for
transient dark structures (e.g., blobs, streaks, or lines) running
perpendicular to the flare ribbons in both 171 and 304. Since 171 is the next
coolest channel to 304, it was often possible to see a brightening form and
then dim in 171, followed by a cospatial signal in 304. For the handful of
highly ambiguous cases, we used running difference movies to capture
precipitating condensation motions that were otherwise masked by more random
activity in the region.

Figure 1: Left: Screen capture of animation (Movie 1, 5 seconds in duration,
available online showing 108 minutes of flare arcade evolution) showing the
flare arcades from an M6.5 flare in SDO AIA 171 Å. An example cool
condensation is indicated in absorption with the blue arrow, and in emission
with a black arrow. Right: The same flare as seen in Movie 1 in 304 Å, with
adjusted scaling to resolve the condensations. An example post-flare arcade
rain blob is indicated in emission with a white arrow, which is not seen in
the left-hand panel, while the black arrow corresponds to an analogous
condensation signature as that indicated by the black arrow in 171.

The last but most pervasive challenge when observing medium-to-large flares
did not involve problems with observing the condensations themselves, but
rather how to define when flare-related rain transitioned to the typical
active region quiescent rain. For many of the flares analyzed here, the arcade
loops did not disappear after the flare; they either began raining again, or
continued raining without pause for tens of hours after all other
observational signatures of the flare had ended. An example of this can be
seen in the movie associated with Figure 1. As we discussed in Section 1, the
most logical, physically-motivated dividing line was the point at which the
arcade reached its apparent maximum height or the point when the flare ribbons
ceased evolving (these being physically equivalent), which both imply that the
flare reconnection forming the arcade has ended. Within a few minutes, the
last arcade-coupled rain should fall, and subsequent rain should be motivated
by more quiescent loop conditions.

Due to the limitations inherent in remote EUV observations with a resolution
lower than that of the smallest condensations, there are likely to be smaller
rain events that are unresolved or lost in the background signal. There are
also line-of-sight constraints which make it difficult to determine when the
coronal rain has reached the chromosphere. None of the uncertainties presented
in this section are unique to this study, but are indeed shared by all remote
solar observations. While it is possible to wait for a prime alignment to
observe when studying long-duration structures, the impulsive nature of solar
flares often means taking measurements at imperfect angles. This study,
therefore, represents a uniformly-analyzed body of flares which were chosen
for the best possible statistical coverage, whose condensations were
identified by eye and were categorized with a particularly restrictive
definition in order to provide the most physically meaningful results.

We also determined the eruptivity of each solar flare event (whether or not a
CME was released) through HEK following flare events of interest. HEK was
queried for 1-hour-long spans beginning at the flare start time, and we then
assessed the returned CME or CMEs to determine whether the location correlated
with the flare within the respective time frame. CME candidates for these
flares were visually identified and confirmed through the SOHO LASCO CME
Catalog.111This CME catalog is generated and maintained at the CDAW Data
Center by NASA and The Catholic University of America in cooperation with the
Naval Research Laboratory. SOHO is a project of international cooperation
between ESA and NASA.

##  3 Statistical Results

Table 1 summarizes the proportions of flares which did and did not host post-
flare arcade rain. Only two X-class flares out of the 44 did not produce
coronal post-flare arcade rain; one of these, on 2012 March 5 at 02:30 UTC,
was a large and unusual flare near the limb that simply did not form an
arcade, but which was followed several hours later by a second flare with a
classic arcade and several hours of post-flare rain (the second flare was not
part of our analysis). This pair of flares is worthy of further study, to
determine how they were related, and what magnetic changes would lead to their
unique characteristics. The other (2013 November 8, 4:20 UTC) was in a small
confined jet on the border of NOAA active region 11890, which did not exhibit
an arcade large enough to be able to detect individual condensations. The
preponderance of M-class flare arcades also exhibited post-flare arcade rain,
although there is a slight drop of rain frequency between the X and M classes.
C-class flares break this trend by showing a much lower propensity for
condensation formation; only $26\%$ of these flares in the data set hosted
post-flare arcade rain. Furthermore, this tendency showed an uneven
distribution across the C flare magnitudes: the bottom-right panel of Figure 3
shows that there is significantly more rain observed in flares of intensity C5
or higher. As we discuss in the next section, this implies a minimum energy
deposition required for condensations to form.

Table 1: Total Flares Analyzed by Post-Flare Arcade Rain Occurrence |  Post-Flare Arcade Rain |  No Post-Flare Arcade Rain
---|---|---
X |  42 (95%) |  2 (5%)
M |  56 (72%) |  22 (28%)
C |  31 (26%) |  87 (74%)
Totals |  129 (54%) |  111 (46%)

Figure 2 shows the overall correlation between GOES class and the duration of
post-flare rain events. X-class flares have a broad distribution between the
<1 to 6 hour bins; the average rain duration was 2 hours 51 minutes. M-class
flares show a strong peak of 58 in the <1-3 hour range (average value 1 hour
14 minutes), which tails off precipitately for longer post-flare arcade rain
events, while C-class flares have 87 events with no post-flare arcade rain and
30 with only short post-flare arcade rain “storms” (average value 15 minutes);
there was only one C-class flare with post-flare arcade rain slightly over 3
hours in length, and none longer. Figure 3 shows these correlations in more
detail by graphing the time spans of post-flare arcade rain events for each
class separately. The time bins there are uniform between X and M-classes, but
are broken down into smaller time spans for the C-class flares. The
correlation between flare class and post-flare arcade rain duration appears
clearly even with a data set that spans a large representative sample of the M
and C-class catalogs, so we consider the relationship strong and significant.

Figure 2: Graph summarizing the duration of post-flare rain events by GOES
class. X-class flares have a wide distribution of flare rain durations, while
M-class flares’ post-flare arcade rain most often lies in the <1-3 hour range.
C-class post-flare arcade rain never exceeded 3 hours in duration in the
sample studied here (the bins are not end-inclusive; for a more detailed
breakdown of C-class rain times, see Figure 3d).

(a) (b)

(c) (d)

Figure 3: Top left: histogram presenting the percentage frequency of X-class
flares hosting post-flare arcade rain of duration ranging from under an hour
to over 6

hours. Top right: analogous histogram for M-class flares. Bottom left:
analogous histogram for C-class flares; please note that these are subdivided
into shorter time bins only ranging between <1 hour and 4 hours, due to the
absence of any longer post-flare arcade rain events detected in the survey.
Bottom right: the percentages of flares which post-flare arcade rain by
strength within the C-class designation, showing the strong change in post-
flare arcade rain occurrence around C5.

Given the very wide range of post-flare arcade rain durations which we found
in the sample, we investigated whether there was a relationship between the
length of the flare itself, as measured by the GOES X-ray intensity, and the
duration of the subsequent post-flare arcade rain. Figure 4 shows the results
of this investigation; we applied a linear regression to the full set of each
class of raining flare and assessed the fit via the $r^{2}$ measure. None of
the data shows a statistically significant relationship, although the M-class
flares show the least-poor correlation, with an $r^{2}$ value of 0.2486.
X-class flares tended to have long post-flare arcade rain events, regardless
of the length of the flare, while C-class flares had a wide flare duration
spread but very short rain durations. There has been some debate about the
actual duration of C-class flares, and several background-subtraction methods
have been put forward (discussed recently in Reep & Knizhnik, 2019). Given the
very low rates of flare rain across all C-class flares, however, it is
unlikely that background subtracting these low-energy flares would
significantly improve the correlation.

Figure 4: Graph of flare duration, as determined by the GOES X-ray data, vs.
post-flare arcade rain duration. Linear fits for each class show that all are
generally a positive relationship, but none are statistically significant.
M-class flares have the highest correlation, with an r-squared value of
approximately 0.25.

Figure 5: This graph presents the proportions of X, M, and C class flares
which are associated with CMEs and post-flare rain. There are few statistical
conclusions that can be drawn between the raining and eruptive categories, due
to the low rates of non-raining flares in X- and M-class flares, and the low
rates of eruptive flares in C-class flares.

Figure 6: Top: a histogram summarizing the relationship between post-flare
arcade rain occurrence and the Hale class of the source active region (the
bins along the x-axis), color-coded by GOES flare class. Bottom: an analogous
histogram showing eruptivity of flares and the source region’s Hale class.
Both show strong preference for $\beta$ type active regions.

Figure 7: Top: This graph shows the relationship between the proportion of
post-flare arcade rain occurrences (by GOES class) and the solar cycle as
plotted by the annual sunspot number. Bottom left: Sunspot number versus the
percentage of flares in each class that produce rain. Linear fit lines have
been plotted and their respective $r^{2}$ values are shown on the left side of
the graph. All of are reasonably good fits, ranging from about 0.55 (C-class)
to 0.74 (M-class). These may prove to be useful tracers for the solar cycle.
Bottom left: This graph shows the residuals for all three classes of flare;
the distribution of all three are random, showing that the linear fit is a
good choice for these relationships.

Next, we explored a potential connection between the eruptivity of a flare and
whether it produced post-flare arcade rain, since both are well-known to
correlate separately with GOES class. Figure 5 shows a nested pie chart for
all three classes of flares, indicating the proportions of the flares which
fell into the classifications CME/post-flare arcade rain, CME/no post-flare
arcade rain, no CME/post-flare arcade rain, and no CME/no post-flare arcade
rain. As has already been well-established by Kawabata et al. (2018) and
references therein, most X (here $80\%$) and M-class flares ($71\%$) are
eruptive; however, only $13\%$ of the C-class flares were associated with
CMEs.

To try to establish whether there was a correlation between CMEs and post-
flare arcade rain events, we compared the proportions of post-flare arcade
rain/CME and no post-flare arcade rain/CME with no post-flare arcade rain/no
CME and post-flare arcade rain/no CME. The results of this are found in Table
2. There is a statistically insignificant correlation between eruptivity and
post-flare arcade rain occurrence for C-class flares (Fisher’s exact test
gives a value of 0.204). M-class flares show a nearly statistically-
significant relationship between rain presence and eruptivity, with a Fisher
test value of 0.084; X-class flares have the poorest relationship, with a
Fisher test value of 0.704. We chose the Fisher exact test due to the
proportionately small number of non-eruptive flares in X- and M-classes in
general, but would still caution against drawing definitive conclusions given
the low numbers in the non-eruptive category. While intriguing, further study
would be required to show that these relationships are more than simple
correlations.

We also compared the Mt. Wilson Hale class of each source active region across
the flare classes and against the flares’ eruptivity. The results of this
investigation are shown in Figure 6; there are a few flares of each class in
$\alpha$ type regions, but the overwhelming number of flares of all types
occur in the more complex $\beta$ designation – particularly
$\beta\gamma\delta$ spots, which are commonly known to be the most active
(Benz (2017) and references therein). We saw very similar results when the
Hale class was graphed against flare eruptivity, seen in the right graph of
Figure 6.

Finally, to determine whether the presence of post-flare rain followed the
solar cycle, we plotted the post-flare arcade rain-positive flares in each
class by year as a percent of all flares of that class, against the SILSO
annual sunspot number, obtained from the Royal Observatory of Belgium,
Brussels (https://wwwbis.sidc.be/silso/datafiles). The results can be seen in
Figure 7. All three classes of flares track the solar cycle in rough terms,
with reasonably good fits for positive linear relationships to the solar
cycle, as shown in the bottom left graph of 7. X-class flares have an $r^{2}$
value of 0.67, while M- and C-class flares exhibit values of 0.74 and 0.55,
respectively. The bottom right graph shows that the residuals for each class
are fairly randomly distributed, which means that the linear fit is a
reasonable one for this relationship. Therefore, we conclude that X- and
M-class flares show promise as indicators of solar cycle strength. While the
span covered here did not quite capture a full solar cycle, these data show
moderate levels of post-flare arcade rain before the peak, increasing peaks
that match the sunspot number peak, and then dying away almost entirely across
the years of solar minimum. We discuss this finding in more detail in the next
section.

Table 2: Post-Flare Arcade Rain and CME Occurrence | Post-Flare Arcade Rain | No Post-Flare Arcade Rain
---|---|---
| No CME | CME | No CME | CME
X-class | 7 | 35 | 0 | 2
M-class | 14 | 42 | 9 | 13
C-class | 26 | 5 | 76 | 11

##  4 Implications

The purpose of this paper is to establish a benchmark set of statistics for a
well-known phenomenon which has nevertheless been under-reported. The main
findings of this statistical study can be summarized as follows:

  1. 1.

Post-flare arcade rain is prevalent to the point of ubiquity in X-class flares
(as previously reported, see Benz, 2017), and very prevalent across the
M-class flare population. However, flare duration and post-flare arcade rain
duration are _not_ well-correlated. One notable subset we identified, confined
jets, did not have identifiable post-flare arcade rain despite relatively
strong flares. We suspect that this is a simple matter of insufficient spatial
resolution, since there is plenty of work showing condensations in larger jets
with the same dynamics (Raouafi et al. (2016), Kumar et al. (2019), Mason et
al. (2019)). There is also a distinct, direct relationship between GOES class
and post-flare arcade rain duration, suggesting that more sets of loops
reconnect across the neutral line to form arcades in stronger flares. Post-
flare arcade rain occurs, though with markedly reduced frequency, in C-class
flares, making up only approximately half of those analyzed.

  2. 2.

X- and M-class post-flare arcade rain may be a useful tracker of the solar
cycle. While rain in all flares tracked the solar cycle to a moderate extent,
the highest correlation occurred for M-class flares, followed by X-class
flares. If this correlation is shown to extend across more than one solar
cycle, the rate of raining X- and M-class flares may be useful at the early
stages of a solar cycle to predict the strength of solar maximum, and later in
a cycle to tell how quickly solar minimum is approaching (with the caveat that
weaker flares are more common, particularly as active regions weaken near the
tail of the cycle). Further research is planned to corroborate this finding.

  3. 3.

The presence or absence of post-flare rain is not consistently correlated with
whether the flare was eruptive or not. As pointed out in the discussion
concerning Table 2, X- and C-class flares do not show any relationship between
eruptivity and rain presence; M-class flares do, although it is not quite
statistically significant. More work is required to settle whether the weak
correlation is truly due to the post-flare arcade rain or simply the fact that
the majority of X and M-class flares are eruptive and the majority of C-class
flares are not, producing significantly smaller values in one category in each
case. Regardless, this lack of strong, consistent correlation is perhaps
unsurprising; while the two are connected by the process of reconnection,
which influences both the creation of the CME and the amount of energy
injected into the arcade loops, the dynamics which drive the arcade after its
formation (radiative cooling and conduction) are physically isolated from and
occur at different time scales from those driving CME creation (magnetic
reconnection). By all accounts, flare reconnection occurs in a very short time
period. Once the reconnection has happened, the loop cools in relative
isolation, disconnected from the presence or absence of a CME evolving above
it. Certainly, if one considers the pool of all flares, the strongest are more
likely to host CMEs and more likely to host post-flare arcade rain, and our
findings support this. Within a class, though, the reconnection should be
roughly comparable across all of the flares, and there we do not see a strong
correlation between CME and post-flare arcade rain occurrence. It is entirely
possible for, as an example, an M-class arcade loop to simply cool non-
catastrophically and drain (with a CME or without), and it is also possible
for an M-class arcade loop to cool catastrophically and form post-flare arcade
rain (with a CME or without).

  4. 4.

Post-flare rain can act as a proxy for flare energy release. The near-
universal prevalence of post-flare arcade rain in X and M-class flares,
coupled with its distinct drop in occurrence around C5, point to a direct
energy correlation for condensation formation in arcades. Further research is
required to better define the low-end energy cutoff and other factors which
may affect it. However, as individual post-flare arcade rain condensation
geometry estimates have been carried out before (Antolin et al. (2015)), and
radiative cooling rates are well-established, these observations of post-flare
arcade rain presence and duration can be utilized as a proxy to aid in
calculating well-constrained estimates of the energy released during
impulsive-phase flare reconnection.

  5. 5.

Post-flare arcades can persist for days after the flare, eventually
transitioning from post-flare rain to quiescent active region condensations.
While most C-class flares did not produce arcades, and most X-class flares
produced arcades that faded after a few hours, many M-class and some X-class
flares produced arcades that continued to be visible in 171 Å for several days
after the flare, well after the ribbons had ceased evolving. At this point,
the pronounced region of higher emission at the loops’ tops dissipated, and
the loops proceeded to produce more erratic condensations that were no longer
coordinated with neighboring loops. This shift was subtle but common, and
points to a shift from post-reconnection condensation formation to the
ubiquitous quiescent active region thermal nonequilibrium condensations that
have been widely discussed. To the authors’ knowledge, this has not been
previously reported; the life span of such loops, which do not appear to
reconnect again readily after the flare, may be useful for studying active
region aging.

##  5 Acknowledgements

The authors would like to thank the reviewer for constructive suggestions that
greatly improved the manuscript, and Jeff Reep for suggesting this paper topic
several years ago. EIM was supported during the paper’s early development by
an appointment to the NASA Postdoctoral Program at the Goddard Space Flight
Center, administered by Universities Space Research Association under contract
with NASA, and later by NSF’s Solar-Terrestrial Program (Award No.
AGS-1923377). KK contributed to the research described in this paper through
an internship arranged at the Goddard Space Flight Center through the Naval
Academy’s Summer Internship Program. Both authors gratefully acknowledge
collaboration with the NASA GSFC Internal Scientist Funding Model program
“Magnetic Energy Buildup and Explosive Release in the Solar Atmosphere”.

## References

  * Andrews (2003) Andrews, M. D. 2003, Sol. Phys., 218, 261, doi: 10.1023/B:SOLA.0000013039.69550.bf
  * Antiochos (1980) Antiochos, S. K. 1980, The Astrophysical Journal, 241, 385, doi: 10.1086/158351
  * Antolin et al. (2015) Antolin, P., Vissers, G., Pereira, T. M. D., Rouppe van der Voort, L., & Scullion, E. 2015, Astrophys. J., 806, 81, doi: 10.1088/0004-637X/806/1/81
  * Aschwanden et al. (1985) Aschwanden, M. J., Wiehl, H. J., Benz, A. O., & Kane, S. R. 1985, Sol. Phys., 97, 159, doi: 10.1007/BF00152985
  * Benz (2017) Benz, A. O. 2017, Living Reviews in Solar Physics, 14, 2, doi: 10.1007/s41116-016-0004-3
  * Cane et al. (2002) Cane, H. V., Erickson, W. C., & Prestage, N. P. 2002, Journal of Geophysical Research (Space Physics), 107, 1315, doi: 10.1029/2001JA000320
  * Cane & Lario (2006) Cane, H. V., & Lario, D. 2006, Space Sci. Rev., 123, 45, doi: 10.1007/s11214-006-9011-3
  * Carrington (1859) Carrington, R. C. 1859, MNRAS, 20, 13, doi: 10.1093/mnras/20.1.13
  * Cliver et al. (2012) Cliver, E. W., Ling, A. G., Belov, A., & Yashiro, S. 2012, Astrophys. J.l, 756, L29, doi: 10.1088/2041-8205/756/2/L29
  * Drake & Swisdak (2012) Drake, J. F., & Swisdak, M. 2012, Space Sci. Rev., 172, 227, doi: 10.1007/s11214-012-9903-3
  * Feynman & Hundhausen (1994) Feynman, J., & Hundhausen, A. J. 1994, J. Geophys. Res., 99, 8451, doi: 10.1029/94JA00202
  * Gosling (1993) Gosling, J. T. 1993, J. Geophys. Res., 98, 18937, doi: 10.1029/93JA01896
  * Hodgson (1859) Hodgson, R. 1859, MNRAS, 20, 15, doi: 10.1093/mnras/20.1.15
  * Holman (2016) Holman, G. D. 2016, Journal of Geophysical Research (Space Physics), 121, 11,667, doi: 10.1002/2016JA022651
  * Kawabata et al. (2018) Kawabata, Y., Iida, Y., Doi, T., et al. 2018, Astrophys. J., 869, 99, doi: 10.3847/1538-4357/aaebfc
  * Klimchuk (2019) Klimchuk, J. A. 2019, Sol. Phys., 294, 173, doi: 10.1007/s11207-019-1562-z
  * Kumar et al. (2019) Kumar, P., Karpen, J. T., Antiochos, S. K., et al. 2019, ApJ, 873, 93, doi: 10.3847/1538-4357/ab04af
  * Lemen et al. (2011) Lemen, J. R., Title, A. M., Akin, D. J., et al. 2011, in The Solar Dynamics Observatory (New York, NY: Springer US), 17–40, doi: 10.1007/978-1-4614-3673-7_3
  * Lin & Forbes (2000) Lin, J., & Forbes, T. G. 2000, J. Geophys. Res., 105, 2375, doi: 10.1029/1999JA900477
  * Mason et al. (2019) Mason, E. I., Antiochos, S. K., & Viall, N. M. 2019, Astrophys. J.l, 874, L33, doi: 10.3847/2041-8213/ab0c5d
  * Mason et al. (2021) Mason, E. I., Antiochos, S. K., & Vourlidas, A. 2021, ApJ, 914, L8, doi: 10.3847/2041-8213/ac0259
  * Mittal & Narain (2010) Mittal, N., & Narain, U. 2010, Journal of Atmospheric and Solar-Terrestrial Physics, 72, 643, doi: 10.1016/j.jastp.2010.03.011
  * Müller et al. (2017) Müller, D., Nicula, B., Felix, S., et al. 2017, A&A, 606, A10, doi: 10.1051/0004-6361/201730893
  * Nelson & Melrose (1985) Nelson, G., & Melrose, D. 1985, Solar Radiophysics: Studies of emission from the sun at metre wavelengths, 333
  * Parker (1953) Parker, E. N. 1953, ApJ, 117, 431, doi: 10.1086/145707
  * Pesnell et al. (2011) Pesnell, W. D., Thompson, B. J., & Chamberlin, P. C. 2011, in The Solar Dynamics Observatory (New York, NY: Springer US), 3–15, doi: 10.1007/978-1-4614-3673-7_2
  * Priest (1978) Priest, E. R. 1978, Sol. Phys., 58, 57, doi: 10.1007/BF00152555
  * Raouafi et al. (2016) Raouafi, N. E., Patsourakos, S., Pariat, E., et al. 2016, Space Science Reviews, 201, 1, doi: 10.1007/s11214-016-0260-5
  * Reames (1990) Reames, D. V. 1990, Astrophys. J.s, 73, 235, doi: 10.1086/191456
  * Reames (2013) —. 2013, Space Sci. Rev., 175, 53, doi: 10.1007/s11214-013-9958-9
  * Reames (2017) —. 2017, Solar Energetic Particles, Vol. 932, doi: 10.1007/978-3-319-50871-9
  * Reep et al. (2020) Reep, J. W., Antolin, P., & Bradshaw, S. J. 2020, The Astrophysical Journal, 890, 100, doi: 10.3847/1538-4357/ab6bdc
  * Reep & Knizhnik (2019) Reep, J. W., & Knizhnik, K. J. 2019, ApJ, 874, 157, doi: 10.3847/1538-4357/ab0ae7
  * Ruan et al. (2021) Ruan, W., Zhou, Y., & Keppens, R. 2021, ApJ, 920, L15, doi: 10.3847/2041-8213/ac27b0
  * Sheeley et al. (1983) Sheeley, N. R., J., Howard, R. A., Koomen, M. J., & Michels, D. J. 1983, Astrophys. J., 272, 349, doi: 10.1086/161298
  * Sun et al. (2015) Sun, X., Bobra, M. G., Hoeksema, J. T., et al. 2015, ApJ, 804, L28, doi: 10.1088/2041-8205/804/2/L28
  * Švestka (2001) Švestka, Z. 2001, Space Sci. Rev., 95, 135, doi: 10.1023/A:1005225208925
  * Wild & Smerd (1972) Wild, J. P., & Smerd, S. F. 1972, Annual Review of Astron. and Astrophys., 10, 159, doi: 10.1146/annurev.aa.10.090172.001111
  * Yashiro et al. (2005) Yashiro, S., Gopalswamy, N., Akiyama, S., Michalek, G., & Howard, R. A. 2005, Journal of Geophysical Research (Space Physics), 110, A12S05, doi: 10.1029/2005JA011151
  * Zuccarello et al. (2017) Zuccarello, F. P., Chandra, R., Schmieder, B., Aulanier, G., & Joshi, R. 2017, A&A, 601, A26, doi: 10.1051/0004-6361/201629836