# Supernova Explosions of the Lowest-Mass Massive Star Progenitors

Tianshu Wang Department of Astrophysical Sciences, Princeton University,
Princeton, NJ 08544  Adam Burrows Department of Astrophysical Sciences,
Princeton University, Princeton NJ 08544 and the Institute for Advanced Study,
1 Einstein Drive, Princeton NJ 08540 Tianshu Wang<EMAIL_ADDRESS>
(Accepted XXX. Received YYY)

###### Abstract

We here focus on the behavior of supernovae that technically explode in 1D
(spherical symmetry). When simulated in 3D, however, the outcomes of
representative progenitors of this class are quite different in almost all
relevant quantities. In 3D, the explosion energies can be two to ten times
higher, and there are correspondingly large differences in the 56Ni yields.
These differences between the 3D and 1D simulations reflect in part the
relative delay to explosion of the latter and in the former the presence of
proto-neutron star convection that boosts the driving neutrino luminosities by
as much as $\sim$50% at later times. In addition, we find that the ejecta in
3D models are more neutron-rich, resulting in significant weak r-process and
48Ca yields. Furthermore, we find that in 3D the core is an interesting,
though subdominant, source of acoustic power. In summary, we find that though
a model might be found theoretically to explode in 1D, one must perform
supernova simulations in 3D to capture most of the associated observables. The
differences between 1D and 3D models are just too large to ignore.

(stars:) supernovae: general – (stars:) neutron – (stars:) – hydrodynamics

##  1 Introduction

The emerging theory of core-collapse supernova explosions (CCSNe) suggests
that there is most often a delay in explosion after the bounce of the
Chandrasekhar-like core of most massive-star CCSN progenitors (Janka, 2012;
Burrows & Vartanyan, 2021). During this delay, turbulence driven by neutrino
heating behind the shock in the so-called gain region (Bethe & Wilson, 1985)
builds and becomes more vigorous over time, while the mass accretion rate
through the stalled shock slowly ebbs. At some critical point (Burrows &
Goshy, 1993), either when the silicon-oxygen shell and its associated density
jump is accreted and/or when the mantle between the accretion shock and the
inner core otherwise becomes unstable and can’t maintain quasi-hydrostatic
equilibrium (Wang et al., 2022), the shock is launched. The basic agencies of
explosion are continuing neutrino heating and the turbulent hydrodynamic
stress behind the shock and for most progenitors the turbulence is crucial to
the explosion (Burrows et al., 1995; Janka & Mueller, 1996; Wang et al., 2022;
Burrows et al., 2024). The accumulated energy in the post-shock ejecta, when
corrected for the binding energy of the outer stellar mantle, starts negative
and builds over time, reaching its asymptotic positive value sometime between
$\sim$0.5 and $\sim$10 seconds after bounce. The timing of explosion and its
subsequent development depend centrally on the mass density profile of the
Chandrasekhar-like core at the time of the collapse instability, frequently
associated with a “compactness” parameter (O’Connor & Ott, 2011) 111Defined as
$\frac{M/M_{\odot}}{R(M)/1000\text{km}}$, we have generally set $M$ equal to
1.75 $M_{\odot}$. However, for the compact structures encountered at the
lowest supernova progenitor masses, calculating the compactness at such large
interior masses makes little sense. Hence, on Table 1 we also provide the
initial central density and entropy, quantities better able to distinguish
such models. which increases roughly, though not monotonically (Sukhbold et
al., 2016, 2018), with progenitor mass. The lower compactness progenitors
generally seem to explode with lower energies, lower 56Ni yields, lower recoil
kicks, and more spherically. When higher compactness models explode, they do
so with greater explosion energies, greater 56Ni yields, greater
asphericities, and larger kicks (though the latter with significant scatter;
Burrows et al. (2023)). While there is not yet quantitative agreement among
supernova theorists performing sophisticated 3D simulations (Lentz et al.,
2015; Roberts et al., 2016; Takiwaki et al., 2016; Müller et al., 2017; Ott et
al., 2018; Glas et al., 2019; Müller et al., 2019; Burrows et al., 2019, 2020;
Müller & Varma, 2020; Powell & Müller, 2020; Stockinger et al., 2020; Kuroda
et al., 2020; Bollig et al., 2021; Vartanyan et al., 2022, 2023; Burrows et
al., 2024), there is overall qualitative agreement on the specifics
articulated here. Burrows et al. (2024) summarize the systematic correlations
between observables, compactness, and progenitor mass seen in their
comprehensive suite of state-of-the-art 3D simulations.

However, the 8.8 $M_{\odot}$ (n8.8) solar-metallicity progenitor of Nomoto
(1984) and Nomoto & Hashimoto (1988) and both the 8.1 $M_{\odot}$ (u8.1;
$10^{-4}$ solar metallicity) and 9.6 $M_{\odot}$ (z9.6; zero-metallicity)
models (A. Heger, private communication), with the lowest compactnesses
published, do explode in 1D (Kitaura et al., 2006; Burrows et al., 2007; Janka
et al., 2008; Fischer et al., 2010; Hüdepohl et al., 2010; Melson et al.,
2015). Though this is the realm of the so-called “ONe,” or “electron-capture”
progenitors (Nomoto, 1984; Nomoto & Hashimoto, 1988; Leung et al., 2020;
Limongi et al., 2024), the extremely-low metallicity models u8.1 and z9.6 in
fact have iron cores at the time of collapse; however, all such models
(including n8.8) have NSE (nuclear statistical equilibrium) core compositions
at bounce. For such models that explode in 1D the turbulent stress and 3D
effects are not crucial factors in the inauguration of explosion. This is
because their rarified envelopes lead to a post-bounce mass accretion rate
that plummets immediately after bounce, and with it the accretion ram
pressure, while the early post-bounce neutrino luminosity (most of it sourced
by the diffusive component from the shocked proto-neutron star (PNS) inner
mantle) is adequate to power a prompt, neutrino-driven explosion222For more
massive progenitors with shallower initial density profiles, post-bounce
accretion contributes significantly to the driving neutrino luminosities
themselves., though with low asymptotic supernova energy, and such explosions
transition into a neutrino-driven wind very early on (Wang & Burrows, 2023a).

Figure 1: Progenitor mass density profiles as a function of interior mass
coordinate (left) and radius (right) for the extremely low compactness, low or
zero metallicity, models of this study (u8.1, z9.6), compared with the
corresponding profiles of low-mass solar-metallicity progenitors (s9.0a,
s9.0b, s9.25, s9.5) and an example of a higher-mass progenitor (s16). The 8.8
solar-metallicity progenitor (Nomoto, 1984; Nomoto & Hashimoto, 1988) is also
plotted. The s9.0a, s9.0b, s9.25, s9.5, and s16 comparison models are from
Sukhbold et al. (2016) and Sukhbold et al. (2018) and the extreme compactness
models (u8.1 and z9.6) that explode in 1D are from A. Heger (private
communication). The s9.0 model does not seem to explode in 1D. See text for a
discussion.

Nomoto (1984) introduced the idea that ONe cores (such as his 8.8 $M_{\odot}$
model) achieve the effective Chandrasekhar mass and start to collapse not by
the photodissociation of iron-peak nuclei, but by electron capture, and that
there is burning on infall. However, what type of inner structures the stars
in the $\sim$8$-$9 $M_{\odot}$ progenitor mass range (and for a spread of
metallicities) actually possess at the end of their lives has not been
definitively determined. This is in part due to the many thermal pulsations
such cores experience and their electron degeneracy. The latter can lead to
many flashes that are difficult to simulate, and the former significantly
retard progress to collapse in extant stellar evolution codes. The result is
not only ambiguity, but that there are very few published models in this mass
range available to supernova modelers. Limongi et al. (2024) and Leung et al.
(2020) discuss these issues (see also Jones et al. (2019)) and Limongi et al.
(2024) suggest that solar-metallicity stars in the ZAMS mass range
$\sim$8.2$-$9.2 $M_{\odot}$ will be electron-capture supernovae. However, due
to myriads of intervening thermal pulses at the end of their calculations they
have yet actually to reach the mass density threshold for electron capture on
20Ne.

Nevertheless, despite the remaining uncertainties concerning the core
structures of massive stars in the 8$-$9 $M_{\odot}$ ZAMS mass range, the
existence of theoretical evolutionary models that lead to explosions in 1D
(spherical) enables a study of the effects of 3D through a comparison of 1D
and 3D simulations. This is the subject of this paper. In addition, the
Salpeter mass function indicates that this range of ZAMS mass space occupies a
disproportionate fraction of progenitor space, though remaining uncertainties
in stellar evolution leave the actual fraction a subject of conjecture. Should
all stars between 8 and 9 $M_{\odot}$ technically explode in 1D, then this
would constitute $\sim$15% of massive star deaths and nearly 15% of core-
collapse supernovae.

Figure 2: The temporal evolution of the average shock radius (left) and the
explosion energy (right, in Bethes $\equiv$1051 ergs) for 3D models (solid)
and 1D models (dashed, for only the u8.1, z9.6, and s9.0 progenitors). Note
that the 1D s9.0 model does not explode, while the corresponding 1D models for
the u8.1 and z9.6 models do. Note also that the 3D u8.1 ($10^{-4}$ solar) and
z9.6 (zero metallicity) models explode with a greater asymptotic energy than
the 3D solar-metallicity models s9.0, s9.25, and s9.5. Finally, we call the
reader’s attention to the significant relative delay in the explosion of the
1D u8.1 model and its very low asymptotic explosion energy. See text in §2 for
a discussion.

Figure 3: The temporal evolution of the neutrino luminosities (Top) for the
u8.1 (orange) and z9.6 (blue) models in both 1D and 3D. The net neutrino
heating rate evolution (note the log10 time scale) in the gain region is
depicted in the bottom left. PNS convection enhances the neutrino luminosities
of all species at later times and leads to higher net heating rates (bottom
left). This in turn elevates the neutrino-driven wind mass loss rate (bottom
right) in 3D vis à vis 1D. The higher neutrino luminosities in 1D at early
times are due to the later explosion times of such models and the larger
consequent accretion rates during this more extended pre-explosion phase. This
effect is much larger for the u8.1 model than for the z9.6 model. Figure 4:
The temporal evolution of the kick speeds for the 3D models compared in this
study. The recoil kicks for the extreme models u8.1 and z9.6 are quite low and
are due predominantly to the neutrino contribution. See text in §2 and Burrows
et al. (2023) for a discussion of the broader context.

We begin in §2 by describing the distinctive hydrodynamics of 1D, prompt,
neutrino-driven explosions and compare with the same models in 3D. We also
summarize their explosion energies, the recoil kick speeds, and the
gravitational mass range of their residues. In §3, we discuss sound wave
generation and the characteristics of the acoustic component and in §4 we
provide the distinctive nucleosynthesis. The latter includes the 56Ni, 57Ni,
and 44Ti, and 48Ca yields, as well as their possible weak r-process
production. We conclude with a general discussion of the results and thoughts
concerning future theoretical directions.

##  2 General Hydrodynamics and Results

The radiation/hydrodynamic code we use for this study is Fornax (Skinner et
al., 2019; Vartanyan et al., 2019b; Burrows et al., 2020, 2024). For the 1D
simulations, we use a radial grid with 1024 zones and for the 3D simulations
we use a 1024($r$)$\times$128($\theta$)$\times$256($\phi$) grid. The number of
neutrino energy groups is twelve (for $\nu_{e}$, $\bar{\nu}_{e}$, and
“$\nu_{\mu}$” neutrinos, where the latter bundles $\nu_{\mu}$,
$\bar{\nu}_{\mu}$, ${\nu}_{\tau}$, and $\bar{\nu}_{\tau}$ neutrinos).We
employed the SFHo equation of state (EOS) for all simulations (Steiner et al.,
2013).

We use the u8.1 ($10^{-4}$ solar) and z9.6 (zero metallicity) progenitors from
A. Heger (private communication) as example models that explode in 1D. These
progenitors, other than the Nomoto 8.8 $M_{\odot}$ model, have the lowest
published compactness. We compare these low-mass models with slightly higher-
compactness models that explode only in 3D (s9.0a, s9.0b, s9.25, and s9.5, all
of which have solar metallicity (Sukhbold et al., 2016)). Models s9.0a and
s9.0b differ only in the imposition of perturbations on infall in the former,
with consequences predominantly for the nucleosynthesis (see Wang & Burrows
(2023a, b)). All these comparison models have already been published (Wang &
Burrows, 2023a, b; Burrows et al., 2024). The outer boundary radius is 20,000
kilometers (km) for models u8.1 and z9.6, 30,000 km for models s9.0a, s9.0b
and s9.5, and 45,000 km for model s9.25.

Figure 1 depicts the initial density profiles of all models mentioned above,
together with the solar-metallicity 8.8 $M_{\odot}$ progenitor Nomoto &
Hashimoto (1988) that also explodes in 1D (labeled as n8.8) and a higher-mass
comparison progenitor (16 $M_{\odot}$ with solar metallicity, labeled as s16).
It is clear that models that explode in 1D (z9.6, u8.1, and n8.8) have lower
densities in their mantles and have the steepest mantle density gradients.
Since the density decreases rapidly at interior masses of 1.3$-$1.5
$M_{\odot}$, calculating the compactness at 1.75 $M_{\odot}$, as we normally
do, may not capture the extreme contrasts between core and mantle for those
models that explode in 1D. Therefore, it may not be as useful a discriminant.
Though for those models the compactness is extremely low, using the central
density or entropy at collapse may be better (Temaj et al., 2023) and we
provide these quantities in summary Table 1. However, for now we will not
focus on specific discriminating indices. Rather, we seek to explore the
differences between 1D and 3D developments when the 1D models explode and
leave to future work the detailed classification of the core structures of the
least massive CCSN progenitor stars.

Figure 2 portrays the mean shock radius versus time after bounce for the 1D
(dashed) and 3D (solid) models (left panel) and also provides for comparison
the corresponding development of the 3D models s9.0a, s9.0b, s9.25, and s9.5.
We note that though the 1D u8.1 model explodes, that explosion is delayed
significantly. Its explosion is marginal, with an energy of only $\sim$1049
ergs (see Table 1). This is a direct consequence of its slightly thicker
envelope than found for the z9.6 model (Figure 1). A slightly enhanced initial
mantle density translates into a slightly greater mass accretion rate
($\dot{M}$) and accretion ram pressure after bounce. These delay the explosion
until the accretion rate subsides. Since they are not aided by important 3D
effects, such as neutrino-driven turbulence, in overcoming the ram pressure,
one-dimensional models suffer greater delays. For many progenitor models, such
an eventual decrease in accretion rate is accompanied by a corresponding
decrease in driving luminosity, with the result that explosion is not assured
(Wang et al., 2022), and this is why the s9.0-1D model in Figure 2, and most
other 1D models, aren’t seen to explode. However, for the models that explode
in 1D, the diffusive component of the luminosity provides a floor, so that the
rapid decrease in $\dot{M}$ is not accompanied by a decrease in the driving
neutrino luminosities. The result is explosion, though the asymptotic
explosion energy (Figure 2, right panel), contingent as it is upon the
integrated mantle neutrino heating rate, is generically low.

As the right-hand-side of Figure 2 demonstrates (and Table 1 summarizes), the
final explosion energies for the z9.6 model are $\sim$0.125 Bethes (1D) and
$\sim$0.22 Bethes (3D), while for the u8.1 model they are $\sim$0.017 Bethes
(1D) and $\sim$0.17 Bethes (3D)333One Bethe (B) is 1051 ergs.. Interestingly,
the explosion energies in 3D for the u8.1 and z9.6 models are slightly higher
than the explosion energies for the s9.0, s9.25, and s9.5 progenitors, despite
their larger compactnesses and neutrino luminosities. It is the greater delay
to explosion and the associated wasting of the heat deposited during a longer
stall to explosion of the s9.0, s9.25, and s9.5 progenitors (despite their
slightly higher luminosities at the same post-bounce times) that results in
slightly lower asymptotic explosion energies for them. This provides an
“asterisk” to the conclusion that explosion energy tracks compactness (Burrows
et al., 2024) and is likely a metallicity effect $-$ published lower
metallicity progenitors have higher compactnesses upon collapse (Sukhbold et
al., 2016).

Figure 3 portrays the evolution of the neutrino luminosities of species
$\nu_{e}$, $\bar{\nu}_{e}$, and “$\nu_{\mu}$” of models z9.6 and u8.1 in both
1D and 3D. We note that for the u8.1 model the $\nu_{e}$ and $\bar{\nu}_{e}$
luminosities in 1D start larger than the corresponding luminosities for the 3D
model. This is due to the significant delay to explosion of this 1D model and
the resultant larger mass accretion rate (and corresponding accretion
luminosity) at this earlier phase compared with that found in 3D, which
explodes much earlier. As Figure 1 shows, the envelope of the u8.1 model is
slightly thicker than that for the z9.6 model. As seen in Figure 2, this
results in a quicker explosion in 1D for the latter than the former.

However, at slightly later times, the luminosities of the 3D models overtake
those for the 1D models for both the z9.6 and u8.1 progenitors and continue
higher thereafter. This is a direct consequence of proto-neutron star (PNS)
convection in the residual core and is a purely multi-D effect. Driven by
negative lepton gradients in the PNS (Dessart et al., 2006; Nagakura et al.,
2019), such convection dredges heat from the core, is found during the first
second after bounce in these models between about 15 and 30 km, and continues
unabated for the duration of the PNS cooling phase.

The left lower panel of Figure 3 shows the net neutrino heating rate in the
gain regions for the 1D and 3D models of the u8.1 and z9.6 progenitors. Due to
PNS convection and a slight physical expansion in the shock radius (and,
hence, the gain volume) in 3D due to turbulence behind the initially
decelerating shock, the neutrino heating rates in 3D are generally larger. As
a result, the explosion energies in 3D are larger than in 1D. This is
particularly the case for the u8.1 model due to the significant delay in its
1D explosion. The difference in the wind mass loss rates between the 1D and 3D
models (bottom right of Figure 3) at later times is also quite pronounced,
amounting to approximately a factor of three in $\dot{M}$. This is a major
difference in the behavior of 3D models and emphasizes that, though these
models explode in 1D, one must perform 3D simulations to capture their true
behavior.

Figure 4 shows the evolutions of the PNS recoil kick imparted to the 3D u8.1
and z9.6 models, and compares them to similar curves for the s9.0a,b, s9.25,
and s9.5 models. A much more comprehensive discussion of this general topic
can be found in Burrows et al. (2023). We see in Figure 4 that the kicks are
low, near $\sim$40 km s-1 for the u8.1 model and at $\sim$13 km s-1 for the
z9.6 model. Not shown on this plot is that these modest kicks are due mostly
to anisotropic neutrino emissions (Vartanyan et al., 2019a). The ejecta are
mostly launched spherically and do not contribute to the same degree. That
this is likely the case for the lowest mass progenitors was discussed and
motivated in Burrows et al. (2023).

We provide the final neutron star baryon and gravitational masses for the
residual neutron stars in Table 1, along with a variety of other summary model
characteristics and results. We are using the SFHo EOS and note that the
gravitational mass in particular is contingent upon this. Nevertheless, we
find gravitational masses for the u8.1 and z9.6 neutron stars near a low value
of $\sim$1.25 $M_{\odot}$. Note that the gravitational mass for the residue of
the s9.0 model is expected to be slightly lower still, as is the neutron star
mass for the residues of accretion-induced collapse (Dessart et al., 2007).

For models that explode in 1D (and their 3D counterparts), the hydrodynamics
transitions quickly into a wind phase (Duncan et al., 1986; Burrows, 1987;
Burrows et al., 1995; Qian & Woosley, 1996; Wang & Burrows, 2023a). In fact,
for them the power in the wind contributes substantially to the total energy
of the blast. Figure 5 shows on the same plot the evolution of the various
powers and luminosities associated with the 3D u8.1 and z9.6 models. Included
are the sum of the $\nu_{e}$ and $\bar{\nu}_{e}$ luminosities, the peak of the
convective luminosity in the PNS ($L_{\rm PNS}$), the heating rate in the gain
region ($\dot{Q}$) (repeated from Figure 3), the mechanical luminosity of the
ejecta ($L_{\rm mech}$, which clearly emerges as a wind after $\sim$0.3
seconds), and the acoustic power that emerges from the core into the exploding
mantle (see §3). Note that the plot starts at 0.1 seconds after bounce.

To calculate these quantities we have used the formulae:

$\begin{split}L_{\rm mech}&=r^{2}\int d\Omega\rho
v_{r}(\frac{1}{2}(v_{r}^{2}+v_{\theta}^{2}+v_{\phi}^{2})+e+\Phi+P/\rho)\\\
L_{\rm PNS}&=r^{2}\int d\Omega\rho v_{\rm turb}(\frac{1}{2}(v_{\rm
turb}^{2}+v_{\theta}^{2}+v_{\phi}^{2})+e+\Phi+P/\rho)\\\ L_{\rm
acoustic}&=r^{2}\int d\Omega(P-\langle P\rangle_{t})(v_{r}-\langle
v_{r}\rangle_{t})\,,\end{split}$ (1)

where $\Phi$ is the gravitational potential, $v_{\rm
turb}=v_{r}-\frac{1}{4\pi}\int d\Omega v_{r}$, and $\langle A\rangle_{t}$ is
the time-averaged value of $A$. Intuitively, $L_{\rm mech}$ is the flux of
total energy, $L_{\rm PNS}$ is the energy flux carried by convection in the
PNS (where the radial speed, $v_{\rm turb}$ is the speed relative to the mean
radial motion), and $L_{\rm acoustic}$ is the energy carried by acoustic waves
in the frame of the radial motion. In our calculation of the acoustic power,
we calculate the time average using a 20-ms sliding window (50 Hz). This
filters out the influence of the expanding shock. Changing this frequency cut
to any value between 25 and 80 Hz doesn’t influence the measured acoustic
power (especially at late times), since the acoustic component generally has
frequencies above 80 Hz (§3).

Figure 5 is quite illuminating. First (as suggested in the bottom right panel
of Figure 3), the mechanical power of the late-time wind is approximately
seven and four times higher in 3D than 1D for the z9.6 and u8.1 progenitors,
respectively. This is a direct consequence of the PNS-convection-boosted
neutrino luminosities and heating rates. Second, $\dot{Q}$ and $L_{\rm mech}$
in both 1D and 3D soon follow in lock step in an approximate power-law
relation, with power 1.8 to 1.6 ($L_{mech}\propto\dot{Q}^{1.8-1.6}$), only
gradually/secularly decreasing with time. Third, while during the first second
after bounce the neutrino luminosity clearly exceeds $L_{\rm PNS}$ (early on
by an order of magnitude), the two approach one another at later times as the
accretion phase abates and the PNS settles into its quasi-static Kelvin-
Helmholtz phase. The fluctuations in $L_{\rm PNS}$ are clearly seen in Figure
5444We note that the 3D calculations we have performed over the last $\sim$5
years using Fornax follow 1) to the center, 2) to late times, and 3) in 3D the
full Kelvin-Helmholtz contraction phase of neutron star birth. This is the
only context in stellar evolution theory for which almost an entire stellar
evolutionary phase, with nested convective and radiative zones, has been
followed in 3D..

Figure 5: The neutrino ($L_{\nu_{e}}+L_{\bar{\nu}_{e}}$), mechanical, PNS convective, and acoustic luminosities of the z9.6 and u8.1 models in the neutrino-driven wind phase. The net neutrino heating rate ($\dot{Q}$) for is also shown. The mechanical luminosity is measured at 500 km, the convective luminosity is the maximum value in the PNS convective region, and the acoustic luminosity is measured just above the PNS radius defined at 1011 g cm-3. Note that the plots start at 0.1 seconds after bounce. Also shown for comparison is $L_{\rm mech}$ for the 1D models (dashed) at later times, well into the wind phase. $M_{\text{ZAMS}}$ | Dimension | Metallicity | $\rho_{\rm c,ini}$ | $S_{\rm c,ini}$ | $t_{\rm pb,max}$ | $M_{\rm b,final}$ | $M_{\rm g,final}$ | $E_{\rm exp}$ |  56Ni | $|v_{\rm kick}^{\rm total}|$
---|---|---|---|---|---|---|---|---|---|---
[$M_{\odot}$] |  | Solar Units | [g cm-3] | [$k_{b}$/baryon] | [s] | [$M_{\odot}$] | [$M_{\odot}$] | [B] | [$10^{-2}M_{\odot}$] | [km s-1]
u8.1 | 3D | $10^{-4}$ | $8.61\times 10^{9}$ | 0.842 | 0.844 | 1.363 | 1.250 | 0.171 | 0.572 | 37.8
u8.1 | 1D | $10^{-4}$ | $8.61\times 10^{9}$ | 0.842 | 2.715 | 1.386 | 1.269 | 0.017 | 0.164 | N/A
z9.6 | 3D | $0$ | $5.80\times 10^{9}$ | 0.888 | 0.721 | 1.356 | 1.245 | 0.207 | 0.566 | 13.0
z9.6 | 1D | $0$ | $5.80\times 10^{9}$ | 0.888 | 0.775 | 1.365 | 1.252 | 0.129 | 0.922 | N/A
s9.0(a) | 3D | 1 | $1.63\times 10^{10}$ | 0.532 | 1.775 | 1.347 | 1.237 | 0.111 | 0.168 | 120.7
s9.0(b) | 3D | 1 | $1.63\times 10^{10}$ | 0.532 | 1.950 | 1.348 | 1.238 | 0.0945 | 0.612 | 78.6
s9.25 | 3D | 1 | $1.87\times 10^{10}$ | 0.540 | 3.532 | 1.378 | 1.263 | 0.124 | 1.04 | 140.1
s9.5 | 3D | 1 | $1.78\times 10^{10}$ | 0.540 | 2.375 | 1.397 | 1.278 | 0.142 | 1.47 | 208.6
Table 1: Summary of simulation outcome data and initial model characteristics.
$\rho_{\rm c,ini}$ and $S_{\rm c,ini}$ are the central density and entropy of
the progenitor model (when the maximum infall speed reaches 950 km $s^{-1}$).
$t_{\rm pb,max}$ is the post-bounce duration of the simulation. We stop the
simulations when the shock radii reach about 20,000 km. The u8.1-1D model
explodes later than the other models, and with lower energy, with the result
that it has a significantly lower shock velocity. Thus, this simulation was
run much longer than the others. $M_{\rm b,final}$ is the PNS baryon mass at
the end of each simulation, while $M_{\rm g,final}$ is the derived
gravitational mass of the cold neutron star. $E_{\rm exp}$ is the explosion
energy and $|v_{\rm kick}^{\rm total}|$ is the PNS recoil kick speed at the
end of each simulation. For the lowest mass CCSNe, 56Ni production (here in
units of 10-2 $M_{\odot}$) is lower in the 3D simulations due to the larger
amount of neutron-rich ejecta. The u8.1 model is a special case, since its 1D
counterpart explodes so weakly that the total ejecta mass is much lower than
in 3D. Therefore, u8.1-1D 56Ni production is lower even if the 56Ni mass
fraction is higher than in 3D.

##  3 Sound Generation and Acoustic Power

Nevins & Roberts (2023) have speculated that the nucleosynthetic yields of
such lowest-compactness progenitors could be affected by this acoustic
component. Gossan et al. (2020) have suggested that acoustic power could help
drive the supernova explosion. Figure 6 depicts snapshots at two different
times of the divergence of the velocity in the inner 1000 km of the explosions
of the u8.1 and z9.6 3D models. Roughly concentric waves of velocity
perturbations are clearly seen. Movies of these waves show that these
perturbations propagate at speeds of $c_{s}+v$, where $c_{s}$ is the local
speed of sound and $v$ is the wind speed. After inner turbulence settles
(within $\sim$0.2 seconds of bounce), a Fourier analysis reveals that these
sound waves have frequencies between roughly 100 and 250 Hertz (Hz). We
comment that no such waves are seen in the 1D simulations.

Figure 7 depicts space-time diagrams of the radial velocity along
$(\theta,\phi)=(\pi/2,0)$ (along the x-direction of Figure 6) for both models
u8.1 and z9.6. The colors represent the magnitude of the radial velocity. The
region of PNS convection, the epoch and realm of early turbulent convection,
and the striations that track the emergent sound waves are all encapsulated in
this plot. The black line follows the PNS radius, defined as where the mean
density equals 1011 g cm-3.

However, though sound waves are unambiguously generated (Figures 6 and 7), the
inferred luminosity is quite small. To derive this quantity we subtract out
the mean flow to obtain the pressure and speed perturbations and then use the
last (sub)equation in Eq. 1 to obtain the sonic power. Figure 5 demonstrates
that after $\sim$0.3 seconds after bounce and just outside the PNS radius (see
Figure 7), the emergent acoustic power is small, about five orders of
magnitude below $L_{\rm PNS}$. However, at earlier times our estimate of the
acoustic power emerging from deeper in gets as high as $\sim$1048-49 ergs s-1.
The turbulence in this region at these earlier times makes it quite difficult
to derive this quantity (see eq. 1 and Figure 7), since extracting the
relevant pressure and velocity perturbations in this region at this early
phase before the flow settles into a smooth, approximately spherical, wind is
messy. Nevertheless, it might be that the core generated sound power reaches
early on within three orders of magnitude of the peak mechanical luminosity in
the PNS convection region and four orders of magnitude of
$L_{\nu_{e}}+L_{\bar{\nu_{e}}}$, and that it is absorbed within a few tens of
kilometers of the region of PNS convection. Wave damping would likely be due
to non-linear breaking and/or neutrino cooling (Gossan et al., 2020), and if
the former could very marginally contribute to powering the supernova and
subsequent wind.

Nevertheless, our results challenge the notion that sound generated by tapping
PNS convection might be a major factor in powering the supernova or in
nucleosynthesis (Nevins & Roberts, 2023; Gossan et al., 2020). However, we
have yet to extract the corresponding quantities for other progenitor models.
All we can say currently is that we derive quite low acoustic powers for the
u8.1 and z9.6 models. It would, however, be illuminating to determine whether
higher spatial resolution simulations yield a similar conclusion.

Figure 6: Velocity divergence at 240 ms and 480 ms post-bounce on the x-y
plane of the 3D z9.6 and u8.1 models. The presence of sound waves generated in
the interior is clear. See text in §3 for a discussion.

Figure 7: Space-time diagrams of radial velocity along
$(\theta,\phi)=(\pi/2,0)$ (along the x-direction of Figure 6). Top: Radial
velocity of the z9.6 model. Bottom: Same as the top panel, but for the u8.1
model. The solid black line shows the PNS radius (defined as the
$\rho=10^{11}$g cm-3 isosurface).

##  4 Nucleosynthesis

The origin of the elements is the primary topic of investigation in nuclear
astrophysics and has been studied in great detail for more than seventy years.
A rich literature (Thielemann et al., 1996; Rauscher et al., 2002; Woosley et
al., 2002; Heger & Woosley, 2010; Nomoto et al., 2013; Sukhbold et al., 2016;
Limongi & Chieffi, 2018; Curtis et al., 2019) has been created that
encompasses pre-supernova stellar evolution, explosive nucleosynthesis,
$\alpha$-rich freeze-out (Woosley et al., 2002; Nomoto et al., 2013; Arcones &
Thielemann, 2023; Wang & Burrows, 2023b), post-explosion winds (Duncan et al.,
1986; Burrows, 1987; Burrows et al., 1995; Qian & Woosley, 1996; Wanajo et
al., 2018; Wang & Burrows, 2023a), and neutrino nucleosynthesis (Woosley et
al., 1990; Fröhlich et al., 2006; Pruet et al., 2006) that forms the core of
our general understanding of the sources of the isotopes of Nature. However,
recent developments (Lentz et al., 2015; Wongwathanarat et al., 2015; Harris
et al., 2017; Wanajo et al., 2018; Sieverding et al., 2020, 2023; Wang &
Burrows, 2023a, b) highlighting the differences between the traditional 1D
supernova nucleosynthesis calculations with ad hoc explosions in spherical
symmetry and more modern 3D perspectives informed by self-consistent
simulations of the explosion, complicated debris trajectories, and the
associated nucleosynthesis are emerging to deepen, if also complicate, this
narrative.

In this context, and for this paper, we have sought to determine the
differences between 1D and 3D nucleosynthetic realizations for the u8.1 and
z9.6 progenitors as representatives of those (rare) models that explode in
current theoretical state-of-the-art 1D simulations. We will not here be
exhaustive, but will emphasize only the most salient differences we find and
will defer to a later work in the context of a larger suite of progenitor
models a more detailed investigation. We use the procedures outlined in Wang &
Burrows (2023b) to perform the nucleosynthetic calculations.

The left-hand side of Figure 8 depicts the final electron fraction ($Y_{e}$)
distributions for both the 1D (dashed) and 3D (solid) simulations of the u8.1
and z9.6 models. The differences are stark. The 1D models generally explode a
bit later and less vigorously. As a consequence, $\nu_{e}$ and $\bar{\nu}_{e}$
absorption on the ejecta has time to push the $Y_{e}$ into a predominantly
proton-rich state, with very little neutron-rich matter ejected. In contrast,
the 3D models explode earlier and more energetically, with the consequence
that $\nu_{e}$ and $\bar{\nu}_{e}$ absorption does not have sufficient time to
elevate initially neutron-rich matter near the PNS surface into a very proton-
rich state. Moreover, the higher late-time neutrino luminosities in 3D (Figure
3) and the longer dwell times of a fraction of the ejecta result in a spread
in the $Y_{e}$ on the proton-rich side to even higher values. The net outcome
is a pronounced widening of the $Y_{e}$ distributions of the 3D models vis à
vis the 1D models, with significant consequences for the emergent
nucleosynthesis. At the same time, the asymptotic entropy distributions
(right-hand panel of Figure 8) are not dramatically different between the 1D
and 3D realizations. There are differences, however, in the entropy
distributions of the u8.1 and z9.6 models themselves.

Figure 9 shows one of the major consequences. On the left hand panel is the
mass fraction distribution of the ejecta versus atomic mass number ($A$). On
its right-hand panel is the corresponding production factor (see Wang &
Burrows (2023b)). For the lower $A$s, most of the nucleosynthesis is
accomplished prior to the supernova during stellar evolution, and the models
are quite similar. However, for the higher $A$s on these plots there is a
pronounced difference in the yields. The 1D models do not generate much near
the first peak of the r-process (a “weak” r-process), while the 3D models
clearly do, and this is a direct result of the presence of significant
neutron-rich ejecta in them. This conclusion echoes a similar conclusion in
Wang & Burrows (2023a) concerning the nucleosynthetic differences between the
s9.0a and s9.0b models. The former (with imposed initial perturbations that
kick start turbulence) experienced a more rapid explosion that left enough
neutron-rich ejecta to result in significant weak r-process production555We
note that Wang & Burrows (2023a) did not find that imposed perturbations made
much of a nucleosynthetic difference for those models with a greater post-
bounce delay before explosion.. These results (re)emphasize the potential
importance of the lower mass CCSN progenitors studied in Wang & Burrows
(2023a) and Wang & Burrows (2023b) as contributing sources of the first peak
of the r-process. We provide the corresponding 56Ni yields in Table 1.

Curiously, the 48Ca yields for the 3D u8.1 and z9.6 models are
$\sim$1.58$\times$10-5 $M_{\odot}$ and $\sim$3.31$\times$10-5 $M_{\odot}$,
respectively. This reflects the neutron-richness (Figure 8) of a large
fraction of their ejecta and emphasizes, as do the results for the s9.0a
model, the possibility that this most neutron-rich of stable nuclei might
originate in the explosions of the lowest mass massive stars. The
corresponding 44Ti yields are $\sim$3.5$\times$10-6 $M_{\odot}$ and
$\sim$2.9$\times$10-6 $M_{\odot}$, respectively, with the values still growing
at the end of the simulation (Wang and Burrows, in preparation). This
emphasizes again the important role played by the neutrino-driven winds in
44Ti production, as discussed in Sieverding et al. (2023); Wang & Burrows
(2023b).

For 57Ni, the yields are $\sim$1.5$\times$10-4 $M_{\odot}$ and
$\sim$1.4$\times$10-4 $M_{\odot}$ for 3D models u8.1 and z9.6, respectively,
resulting in 56Ni to 57Ni ratios of $\sim$40\. This is near the solar value of
$\sim$43\. Moreover, the stable Ni (58Ni, 60Ni, and 62Ni) to 56Ni ratio of the
ejecta is around $\sim$1.0$-$1.5 for both 3D models, or 15$-$30 relative to
the solar ratio of $\sim$0.057 (Karrer et al., 2007). If the u8.1 and z9.6
models had had solar metallicity envelopes, the ratio of the aggregate nickel
to iron ejecta (after 56Ni decay) would have been $\sim$5$-$6, and this may be
what we might expect for a solar-metallicity model that theoretically would
have exploded in 1D. However, such a conclusion must await a proper simulation
of the relevant progenitor.

Figure 8: Electron fraction and entropy distribution of the 1D and 3D models
u8.1 and z9.6. The differences in the entropy distributions of the 1D and 3D
models are not large, but, due to differences in the explosion energies and
mantle densities of the u8.1 and z9.6 models, differences in the associated
maximum entropies emerge. Moreover, there is a tendency for the peak entropy
to be slightly larger in 3D. Note that the more vigorous and earlier explosion
in 3D leads to large differences in the $Y_{e}$ distributions with the
corresponding 1D results. In particular, the 3D ejecta contain a wider range
of $Y_{e}$s and much more neutron-rich matter.

Figure 9: Ejecta abundances (left) and production factors (right) for the 1D
(dashed) and 3D (solid) u8.1 (orange) and z9.6 (blue) progenitors. Note that
only the 3D models manifest significant “weak” r-process production near
$A=90$.

##  5 Conclusions

In this paper we have focused on the behavior of the cores of representative
massive stars that technically explode in spherical symmetry. According to
current theory, this is a small subclass of supernova progenitors with very
tenuous mantles that at solar-metallicity might have ONe cores at the onset of
collapse and inhabit the $\sim$8$-$9 $M_{\odot}$ mass range (Nomoto, 1984;
Limongi et al., 2024). These would be so-called “electron-capture” supernovae.
The classic such model is the solar-metallicity 8.8 $M_{\odot}$ model of
Nomoto (1984), which was shown originally to explode in 1D by Kitaura et al.
(2006) and in 2D by Burrows et al. (2007). However, as the behavior of models
u8.1 and z9.6 demonstrates, the mass range at extremely low metallicity for
which models theoretically explode in 1D may be broader and such progenitors
may in fact have iron cores at the onset of collapse.

When simulated in 3D with a modern supernova simulation code, we find that the
observable outcomes of representative progenitors of this class (for this
paper chosen to be the u8.1 and z9.6 models of A. Heger) are quite different
than found in 1D in almost all relevant quantities. In particular, we find for
the u8.1 and z9.6 progenitors that in 3D the explosion energies are 0.171 B
and 0.207 B and the 56Ni yields are 0.572$\times$10-2 $M_{\odot}$ and
0.566$\times$10-2 $M_{\odot}$, while in 1D the corresponding quantities are
0.017 B and 0.129 B and 0.164$\times$10-2 $M_{\odot}$ and 0.922$\times$10-2
$M_{\odot}$ (see Table 1). The large differences in the explosion energy and
56Ni yield between the 3D and 1D simulations of the u8.1 progenitor reflect in
part the significant delay to explosion of the latter, associated with its
slightly denser mantle (Figure 1). However, the explosion energies are
significantly larger in 3D in both progenitors, a consequence not only of the
shorter delay to explosion (most prominant for the u8.1 model), but also of
the presence in the core of PNS convection (suppressed in 1D). Such lepton-
driven convection (Dessart et al., 2006; Nagakura et al., 2020) boosts the
driving neutrino luminosities by as much as $\sim$50% at later times (Figure
3) 666We note that Melson et al. (2015) performed a 3D simulation of the same
z9.6 model and obtained an explosion energy of 0.08$-$0.1 B by the end of
their simulation (taken to 0.4 seconds after bounce). Stockinger et al. (2020)
also simulated the z9.6 model in 3D and obtained an explosion energy of
$\sim$0.086 B. These explosion energies are a factor of $\sim$two lower than
what we obtain and the detailed origin of this difference is unknown. Both
Melson et al. (2015) and Stockinger et al. (2020) used the LS220 nuclear EOS
(not the SFHo EOS), but we think it unlikely that this alone could explain the
difference. Slight differences in the neutrino-matter coupling rates might
translate non-linearly into a quantitative difference in the asymptotic
energy. Whatever the reason, it would be useful to determine the origin of
this quantitative difference.. The range of final gravitational masses is
narrow, and for the 3D models is $\sim$1.25 $M_{\odot}$. This number reflects
not only the initial stellar model core structure, but the nuclear equation of
state (SFHo) employed in our simulations.

We find that the distribution of the electron fraction ($Y_{e}$) of the debris
is significantly broadened for the 3D models with respect to their 1D
counterparts. The more neutron-rich 3D ejecta contain significant weak
r-process and 48Ca yields that may be important features of this lower end of
the CCSN progenitor mass range. This has also been seen in the solar
metallicity s9.0a model (Wang & Burrows, 2023b, a). In addition, we find that
the stable nickel to iron ratio of these models is $\sim$20, and is much
higher than that measured for the Crab Nebula (T. Temim, private
communication). However, the corresponding ratio for our lower-mass/solar-
metallicity/iron-core progenitor models s9.0a and s11 do comport with the Crab
measurements, though the ratios for most of the rest of our recent 3D
simulation suite (Burrows et al., 2024) are roughly solar. What this may mean
will be discussed in Wang & Burrows (in preparation).

With such tenuous mantles, neutrino-driven mechanical winds emerge quickly and
are a major factor in the overall supernova energy. Importantly, after a few
tenths of a second after bounce, $L_{\rm mech}$ is $\sim$0.25 and $\sim$0.15
times smaller in 1D than 3D for the u8.1 and z9.6 progenitors, respectively.
We also find that $L_{mech}\propto\dot{Q}^{1.8-1.6}$, where $\dot{Q}$ is the
total neutrino energy deposition rate at the base of the wind.

Coupled with the results for a much broader spectrum of progenitors (Burrows
et al., 2024), the results we witness in this paper support the idea that
there is a continuity in dynamical behavior and outcomes across the entire
progenitor continuum and that the lower mass realm is not qualitatively
distinct in terms of the explosion mechanism or dynamics. However, there may
be an interesting stellar metallicity dependence of the explosion energy at
the low mass end, with lower metallicity explosions for these progenitors
having $\sim$twice the explosion energy at the same low ZAMS mass. This
conclusion, however, has yet to be made robust and will require further study.

Interestingly, we find that the cores of the 3D models (as opposed to 1D
models) radiate acoustic power that may play a modest, though quite sub-
dominant, role in supernova and wind dynamics (Gossan et al., 2020; Nevins &
Roberts, 2023). This component seems sourced in the fluctuations associated
with inner PNS convection and has a frequency range between $\sim$100 and
$\sim$250 Hz. Whether acoustic power plays an important role anywhere across
the supernova progenitor continuum is unlikely, but this has yet to be
investigated in adequate detail.

In conclusion, we have found that though a model might be found theoretically
to explode in 1D (spherical symmetry), one must perform modern supernova
simulations in 3D to capture most of the associated observables. The
differences between 1D and 3D models are just too large to ignore.

## Acknowledgments

We thank Tea Temim, David Vartanyan, Chris White, and Matt Coleman for
previous conversations and/or ongoing collaborations and Alex Heger for
sharing his u8.1 and z9.6 progenitor models. We acknowledge support from the
U. S. Department of Energy Office of Science and the Office of Advanced
Scientific Computing Research via the Scientific Discovery through Advanced
Computing (SciDAC4) program and Grant DE-SC0018297 (subaward 00009650) and
support from the U. S. National Science Foundation (NSF) under Grants
AST-1714267 and PHY-1804048 (the latter via the Max-Planck/Princeton Center
(MPPC) for Plasma Physics). Some of the models were simulated on the Frontera
cluster (under awards AST20020 and AST21003), and this research is part of the
Frontera computing project at the Texas Advanced Computing Center (Stanzione
et al., 2020). Frontera is made possible by NSF award OAC-1818253.
Additionally, a generous award of computer time was provided by the INCITE
program, enabling this research to use resources of the Argonne Leadership
Computing Facility, a DOE Office of Science User Facility supported under
Contract DE-AC02-06CH11357. Finally, the authors acknowledge computational
resources provided by the high-performance computer center at Princeton
University, which is jointly supported by the Princeton Institute for
Computational Science and Engineering (PICSciE) and the Princeton University
Office of Information Technology, and our continuing allocation at the
National Energy Research Scientific Computing Center (NERSC), which is
supported by the Office of Science of the U. S. Department of Energy under
contract DE-AC03-76SF00098.

## Data Availability

The summary numerical data underlying this article will be shared upon
reasonable request to either author.

## References

  * Arcones & Thielemann (2023) Arcones, A., & Thielemann, F.-K. 2023, A&A Rev., 31, 1, doi: 10.1007/s00159-022-00146-x
  * Bethe & Wilson (1985) Bethe, H. A., & Wilson, J. R. 1985, ApJ, 295, 14, doi: 10.1086/163343
  * Bollig et al. (2021) Bollig, R., Yadav, N., Kresse, D., et al. 2021, The Astrophysical Journal, 915, 28, doi: 10.3847/1538-4357/abf82e
  * Burrows (1987) Burrows, A. 1987, ApJ, 318, L57, doi: 10.1086/184937
  * Burrows et al. (2007) Burrows, A., Dessart, L., & Livne, E. 2007, in American Institute of Physics Conference Series, Vol. 937, Supernova 1987A: 20 Years After: Supernovae and Gamma-Ray Bursters, ed. S. Immler, K. Weiler, & R. McCray, 370–380, doi: 10.1063/1.3682931
  * Burrows & Goshy (1993) Burrows, A., & Goshy, J. 1993, ApJ, 416, L75, doi: 10.1086/187074
  * Burrows et al. (1995) Burrows, A., Hayes, J., & Fryxell, B. A. 1995, ApJ, 450, 830, doi: 10.1086/176188
  * Burrows et al. (2019) Burrows, A., Radice, D., & Vartanyan, D. 2019, MNRAS, 485, 3153, doi: 10.1093/mnras/stz543
  * Burrows et al. (2020) Burrows, A., Radice, D., Vartanyan, D., et al. 2020, Monthly Notices of the Royal Astronomical Society, 491, 2715–2735, doi: 10.1093/mnras/stz3223
  * Burrows & Vartanyan (2021) Burrows, A., & Vartanyan, D. 2021, Nature, 589, 29, doi: 10.1038/s41586-020-03059-w
  * Burrows et al. (2024) Burrows, A., Wang, T., & Vartanyan, D. 2024, arXiv e-prints, arXiv:2401.06840, doi: 10.48550/arXiv.2401.06840
  * Burrows et al. (2023) Burrows, A., Wang, T., Vartanyan, D., & Coleman, M. S. B. 2023, arXiv e-prints, arXiv:2311.12109, doi: 10.48550/arXiv.2311.12109
  * Curtis et al. (2019) Curtis, S., Ebinger, K., Fröhlich, C., et al. 2019, ApJ, 870, 2, doi: 10.3847/1538-4357/aae7d2
  * Dessart et al. (2006) Dessart, L., Burrows, A., Livne, E., & Ott, C. D. 2006, ApJ, 645, 534, doi: 10.1086/504068
  * Dessart et al. (2007) —. 2007, ApJ, 669, 585, doi: 10.1086/521701
  * Duncan et al. (1986) Duncan, R. C., Shapiro, S. L., & Wasserman, I. 1986, ApJ, 309, 141, doi: 10.1086/164587
  * Fischer et al. (2010) Fischer, T., Whitehouse, S. C., Mezzacappa, A., Thielemann, F. K., & Liebendörfer, M. 2010, A&A, 517, A80, doi: 10.1051/0004-6361/200913106
  * Fröhlich et al. (2006) Fröhlich, C., Martínez-Pinedo, G., Liebendörfer, M., et al. 2006, Phys. Rev. Lett., 96, 142502, doi: 10.1103/PhysRevLett.96.142502
  * Glas et al. (2019) Glas, R., Just, O., Janka, H. T., & Obergaulinger, M. 2019, ApJ, 873, 45, doi: 10.3847/1538-4357/ab0423
  * Gossan et al. (2020) Gossan, S. E., Fuller, J., & Roberts, L. F. 2020, MNRAS, 491, 5376, doi: 10.1093/mnras/stz3243
  * Harris et al. (2017) Harris, J. A., Hix, W. R., Chertkow, M. A., et al. 2017, ApJ, 843, 2, doi: 10.3847/1538-4357/aa76de
  * Heger & Woosley (2010) Heger, A., & Woosley, S. E. 2010, ApJ, 724, 341, doi: 10.1088/0004-637X/724/1/341
  * Hüdepohl et al. (2010) Hüdepohl, L., Müller, B., Janka, H. T., Marek, A., & Raffelt, G. G. 2010, Phys. Rev. Lett., 104, 251101, doi: 10.1103/PhysRevLett.104.251101
  * Janka (2012) Janka, H.-T. 2012, Annual Review of Nuclear and Particle Science, 62, 407, doi: 10.1146/annurev-nucl-102711-094901
  * Janka & Mueller (1996) Janka, H. T., & Mueller, E. 1996, A&A, 306, 167
  * Janka et al. (2008) Janka, H. T., Müller, B., Kitaura, F. S., & Buras, R. 2008, A&A, 485, 199, doi: 10.1051/0004-6361:20079334
  * Jones et al. (2019) Jones, S., Röpke, F. K., Fryer, C., et al. 2019, A&A, 622, A74, doi: 10.1051/0004-6361/201834381
  * Karrer et al. (2007) Karrer, R., Bochsler, P., Giammanco, C., et al. 2007, Space Sci. Rev., 130, 317, doi: 10.1007/s11214-007-9220-4
  * Kitaura et al. (2006) Kitaura, F. S., Janka, H.-T., & Hillebrandt, W. 2006, A&A, 450, 345, doi: 10.1051/0004-6361:20054703
  * Kuroda et al. (2020) Kuroda, T., Arcones, A., Takiwaki, T., & Kotake, K. 2020, The Astrophysical Journal, 896, 102, doi: 10.3847/1538-4357/ab9308
  * Lentz et al. (2015) Lentz, E. J., Bruenn, S. W., Hix, W. R., et al. 2015, ApJ, 807, L31, doi: 10.1088/2041-8205/807/2/L31
  * Leung et al. (2020) Leung, S.-C., Nomoto, K., & Suzuki, T. 2020, ApJ, 889, 34, doi: 10.3847/1538-4357/ab5d2f
  * Limongi & Chieffi (2018) Limongi, M., & Chieffi, A. 2018, ApJS, 237, 13, doi: 10.3847/1538-4365/aacb24
  * Limongi et al. (2024) Limongi, M., Roberti, L., Chieffi, A., & Nomoto, K. 2024, ApJS, 270, 29, doi: 10.3847/1538-4365/ad12c1
  * Melson et al. (2015) Melson, T., Janka, H.-T., & Marek, A. 2015, ApJ, 801, L24, doi: 10.1088/2041-8205/801/2/L24
  * Müller et al. (2017) Müller, B., Melson, T., Heger, A., & Janka, H.-T. 2017, Monthly Notices of the Royal Astronomical Society, 472, 491–513, doi: 10.1093/mnras/stx1962
  * Müller & Varma (2020) Müller, B., & Varma, V. 2020, arXiv:2007.04775 [astro-ph].  http://arxiv.org/abs/2007.04775
  * Müller et al. (2019) Müller, B., Tauris, T. M., Heger, A., et al. 2019, MNRAS, 484, 3307, doi: 10.1093/mnras/stz216
  * Nagakura et al. (2019) Nagakura, H., Burrows, A., Radice, D., & Vartanyan, D. 2019, Monthly Notices of the Royal Astronomical Society, 490, 4622–4637, doi: 10.1093/mnras/stz2730
  * Nagakura et al. (2020) —. 2020, Monthly Notices of the Royal Astronomical Society, 492, 5764–5779, doi: 10.1093/mnras/staa261
  * Nevins & Roberts (2023) Nevins, B., & Roberts, L. F. 2023, MNRAS, 520, 3986, doi: 10.1093/mnras/stad372
  * Nomoto (1984) Nomoto, K. 1984, ApJ, 277, 791, doi: 10.1086/161749
  * Nomoto & Hashimoto (1988) Nomoto, K., & Hashimoto, M. 1988, Phys. Rep., 163, 13, doi: 10.1016/0370-1573(88)90032-4
  * Nomoto et al. (2013) Nomoto, K., Kobayashi, C., & Tominaga, N. 2013, ARA&A, 51, 457, doi: 10.1146/annurev-astro-082812-140956
  * O’Connor & Ott (2011) O’Connor, E., & Ott, C. D. 2011, ApJ, 730, 70, doi: 10.1088/0004-637X/730/2/70
  * Ott et al. (2018) Ott, C. D., Roberts, L. F., da Silva Schneider, A., et al. 2018, ApJ, 855, L3, doi: 10.3847/2041-8213/aaa967
  * Powell & Müller (2020) Powell, J., & Müller, B. 2020, Monthly Notices of the Royal Astronomical Society, 494, 4665–4675, doi: 10.1093/mnras/staa1048
  * Pruet et al. (2006) Pruet, J., Hoffman, R. D., Woosley, S. E., Janka, H. T., & Buras, R. 2006, ApJ, 644, 1028, doi: 10.1086/503891
  * Qian & Woosley (1996) Qian, Y. Z., & Woosley, S. E. 1996, ApJ, 471, 331, doi: 10.1086/177973
  * Rauscher et al. (2002) Rauscher, T., Heger, A., Hoffman, R. D., & Woosley, S. E. 2002, ApJ, 576, 323, doi: 10.1086/341728
  * Roberts et al. (2016) Roberts, L. F., Ott, C. D., Haas, R., et al. 2016, ApJ, 831, 98, doi: 10.3847/0004-637X/831/1/98
  * Sieverding et al. (2023) Sieverding, A., Kresse, D., & Janka, H.-T. 2023, arXiv e-prints, arXiv:2308.09659, doi: 10.48550/arXiv.2308.09659
  * Sieverding et al. (2020) Sieverding, A., Müller, B., & Qian, Y. Z. 2020, ApJ, 904, 163, doi: 10.3847/1538-4357/abc61b
  * Skinner et al. (2019) Skinner, M. A., Dolence, J. C., Burrows, A., Radice, D., & Vartanyan, D. 2019, The Astrophysical Journal Supplement Series, 241, 7, doi: 10.3847/1538-4365/ab007f
  * Stanzione et al. (2020) Stanzione, D., West, J., Evans, R. T., et al. 2020, in Practice and Experience in Advanced Research Computing, PEARC ’20 (New York, NY, USA: Association for Computing Machinery), 106–111, doi: 10.1145/3311790.3396656
  * Steiner et al. (2013) Steiner, A. W., Hempel, M., & Fischer, T. 2013, ApJ, 774, 17, doi: 10.1088/0004-637X/774/1/17
  * Stockinger et al. (2020) Stockinger, G., Janka, H. T., Kresse, D., et al. 2020, MNRAS, 496, 2039, doi: 10.1093/mnras/staa1691
  * Sukhbold et al. (2016) Sukhbold, T., Ertl, T., Woosley, S. E., Brown, J. M., & Janka, H. T. 2016, ApJ, 821, 38, doi: 10.3847/0004-637X/821/1/38
  * Sukhbold et al. (2018) Sukhbold, T., Woosley, S. E., & Heger, A. 2018, ApJ, 860, 93, doi: 10.3847/1538-4357/aac2da
  * Takiwaki et al. (2016) Takiwaki, T., Kotake, K., & Suwa, Y. 2016, MNRAS, 461, L112, doi: 10.1093/mnrasl/slw105
  * Temaj et al. (2023) Temaj, D., Schneider, F. R. N., Laplace, E., Wei, D., & Podsiadlowski, P. 2023, arXiv e-prints, arXiv:2311.05701, doi: 10.48550/arXiv.2311.05701
  * Thielemann et al. (1996) Thielemann, F.-K., Nomoto, K., & Hashimoto, M.-A. 1996, ApJ, 460, 408, doi: 10.1086/176980
  * Vartanyan et al. (2019a) Vartanyan, D., Burrows, A., & Radice, D. 2019a, MNRAS, 489, 2227, doi: 10.1093/mnras/stz2307
  * Vartanyan et al. (2019b) Vartanyan, D., Burrows, A., Radice, D., Skinner, M. A., & Dolence, J. 2019b, MNRAS, 482, 351, doi: 10.1093/mnras/sty2585
  * Vartanyan et al. (2023) Vartanyan, D., Burrows, A., Wang, T., Coleman, M. S. B., & White, C. J. 2023, Phys. Rev. D, 107, 103015, doi: 10.1103/PhysRevD.107.103015
  * Vartanyan et al. (2022) Vartanyan, D., Coleman, M. S. B., & Burrows, A. 2022, MNRAS, 510, 4689, doi: 10.1093/mnras/stab3702
  * Wanajo et al. (2018) Wanajo, S., Müller, B., Janka, H.-T., & Heger, A. 2018, ApJ, 852, 40, doi: 10.3847/1538-4357/aa9d97
  * Wang & Burrows (2023a) Wang, T., & Burrows, A. 2023a, ApJ, 954, 114, doi: 10.3847/1538-4357/ace7b2
  * Wang & Burrows (2023b) —. 2023b, arXiv e-prints, arXiv:2311.03446, doi: 10.48550/arXiv.2311.03446
  * Wang et al. (2022) Wang, T., Vartanyan, D., Burrows, A., & Coleman, M. S. B. 2022, MNRAS, 517, 543, doi: 10.1093/mnras/stac2691
  * Wongwathanarat et al. (2015) Wongwathanarat, A., Müller, E., & Janka, H. T. 2015, A&A, 577, A48, doi: 10.1051/0004-6361/201425025
  * Woosley et al. (1990) Woosley, S. E., Hartmann, D. H., Hoffman, R. D., & Haxton, W. C. 1990, ApJ, 356, 272, doi: 10.1086/168839
  * Woosley et al. (2002) Woosley, S. E., Heger, A., & Weaver, T. A. 2002, Reviews of Modern Physics, 74, 1015, doi: 10.1103/RevModPhys.74.1015


# Entropic extensivity and large deviations in the presence of strong
correlations

Ugur Tirnakli1<EMAIL_ADDRESS>Mauricio Marques2
<EMAIL_ADDRESS>Constantino Tsallis2,3,4<EMAIL_ADDRESS>1Department of
Physics, Faculty of Science, Ege University, 35100 Izmir, Turkey
2Centro Brasileiro de Pesquisas Fisicas and National Institute of Science and
Technology for Complex Systems
Rua Xavier Sigaud 150, Rio de Janeiro 22290-180, Brazil
3 Santa Fe Institute, 1399 Hyde Park Road, Santa Fe, New Mexico 87501, USA
4 Complexity Science Hub Vienna, Josefstädter Strasse 39, 1080 Vienna, Austria

###### Abstract

The standard Large Deviation Theory (LDT) mirrors the Boltzmann-Gibbs (BG)
factor which describes the thermal equilibrium of short-range Hamiltonian
systems, the velocity distribution of which is Maxwellian. It is generically
applicable to systems satisfying the Central Limit Theorem (CLT), among
others. When we focus instead on stationary states of typical complex systems
(e.g., classical long-range Hamiltonian systems), both the CLT and LDT need to
be generalized. We focus here on a scale-invariant stochastic process
involving strongly-correlated exchangeable random variables which, through the
Laplace-de Finetti theorem, is known to yield a long-tailed $Q$-Gaussian
$N\to\infty$ attractor in the space of distributions ($1<Q<3)$. We present
strong numerical indications that the corresponding LDT probability
distribution is given by
$P(N,z)=P_{0}\,e_{q}^{-r_{q}(z)N}=P_{0}[1-(1-q)r_{q}(z)N]^{1/(1-q)}$ with
$q=2-1/Q\in(1,5/3)$. The rate function $r_{q}(z)$ seemingly equals the
relative nonadditive $q_{r}$-entropy per particle, with
$q_{r}\simeq\frac{7}{10}+\frac{6}{10}\frac{1}{Q-1}$, thus exhibiting a
singularity at $Q=1$ and recovering the BG value $q_{r}=1$ in the $Q\to 3$
limit. Let us emphasize that the extensivity of $r_{q}(z)N$ appears to be
verified, consistently with what is expected, from the Legendre structure of
thermodynamics, for a total entropy. The present analysis of a relatively
simple model somewhat mirroring spin-1/2 long-range-interacting ferromagnets
(e.g., with strongly anisotropic XY coupling) might be helpful for a deeper
understanding of nonequilibrium systems with global correlations and other
complex systems.

##  I Introduction

Boltzmann-Gibbs (BG) statistical mechanics yields various important relations.
Still, it is fair to consider the Maxwellian distribution of velocities and
the exponential distribution of energies (BG weight or BG factor) as its most
important fingerprints statphys1 ; statphys2 . These facts mirror the Central
Limit Theorem (CLT) CLT1 ; CLT2  which leads, when the number $N$ of involved
random variables increases indefinitely, to convergence towards Gaussian
distributions, and the Large Deviation Theory (LDT) LDT1 ; Peliti ; LDT2 ;
Lacomte1 ; Touchette2009 ; Lacomte2  which characterizes the speed at which
Gaussians are approached while $N$ increases. To be more precise, the BG
distribution $p_{BG}$ associated with a many-body Hamiltonian ${\cal H}_{N}$
at thermal equilibrium is given by $p_{BG}\propto e^{-\beta{\cal H}_{N}}$
whenever ${\cal H}_{N}$ includes short-range interactions or no interactions
at all. We may then write that $p_{BG}\propto e^{-[\beta{\cal H}_{N}/N]N}$,
where, consistently with thermodynamics, $[\beta{\cal H}_{N}/N]$ is an
intensive quantity. The corresponding LDT statement for a binary stochastic
system with $N$ random variables yielding $n$ times say 0, and $(N-n)$ times
say 1 concerns the probability $P_{N}(n/N>z)\in[0,1]$ of the random variable
$n/N$ taking values larger than a fixed value $z\in\Re$ for increasingly large
values of $N$. Under the hypothesis of probabilistic independence, or similar
settings, we expect $P_{N}(n/N>z)\approx e^{-r_{1}(z)N}$, where the rate
function $r_{1}$ equals a BG relative entropy per particle. Therefore
$r_{1}(z)N$ plays the role of a thermodynamic total entropy which,
consistently with the Legendre structure of classical thermodynamics, is
extensive, i.e., $r_{1}(z)N\propto N$ ($N\gg 1)$.

Within nonextensive statistical mechanics ($q$-statistics for
short)Tsallis1988 ; Santos1997 ; Abe2000 ; WilkWlodarczyk2000 ; Beck2001 ;
BeckCohen2003 ; HanelThurnerGellMann2014 ; EncisoTempesta2017 ;
JizbaKorbelLavickaProcksSvdovaBeck2018 ; JizbaKorbel2019 ;
SuyariMatsuzoeScarfone2020 ; Tsallisbook , we typically tackle with long-
range-interacting Hamiltonian systems, among other strongly correlated ones.
The associated distributions of velocities appear to be $Q$-Gaussians with
$Q>1$ (see, for instance, AnteneodoTsallis1998 ; CirtoAssisTsallis2013 ;
CirtoRodriguezNobreTsallis2018  for the $\alpha$-XY ferromagnet,
RodriguezNobreTsallis2019  for the $\alpha$-Heisenberg ferromagnet, and
ChristodoulidiTsallisBountis2014 ; ChristodoulidiBountisTsallisDrossos2016 ;
BagchiTsallis2016 ; BagchiTsallis2017  for the $\alpha$-Fermi-Pasta-Ulam
model), with $Q$ approaching unity when the range of the interactions
approaches the short-range regime. These facts are to be associated with a
$Q$-Central Limit Theorem ($Q$-CLT) which leads, when $N\to\infty$, to a
convergence on a $Q$-Gaussian distribution. Sufficient conditions for the
$Q$-CLT to hold are already available UmarovTsallisSteinberg2008  (see also
UmarovTsallisGellMannSteinberg ; umarov ) but the necessary conditions for a
$Q$-CLT still remain as a challenge.

Within $q$-statistics we expect, for the total energy of a long-range-
interacting system at its stationary, or quasi-stationary, state, to be super-
extensive, hence, not proportional to $N$, as it is the case for short-range-
interacting systems. More precisely, we expect $p_{q}\propto
e_{q}^{-\beta_{q}{\cal H}_{N}}$, where ${\cal H}_{N}$ is a super-extensive
Hamiltonian, $\beta_{q}$ playing the role of an inverse effective temperature
($\beta_{1}\equiv\beta$ corresponds to the usual inverse kinetic temperature);
$\beta_{q}$ generically differs from $\beta$ (in extreme cases, by orders of
magnitude Andradeetal2010 ). We remind that
$e_{q}^{z}\equiv[1+(1-q)z]^{1/(1-q)}$ with $e_{1}^{z}=e^{z}$. For say two-body
(attractive) interactions decaying like $1/(distance)^{\alpha}$
($\alpha\in[0,\infty)$) within a $d$-dimensional system, we may rewrite
$p_{q}\propto e_{q}^{-[(\beta_{q}\tilde{N})({\cal H}_{N}/N\tilde{N})]N}$ where
$\tilde{N}\equiv\frac{N^{1-\alpha/d}-1}{1-\alpha/d}$ is, for $N$ increasingly
large, constant for $\alpha/d>1$ (short-range), increases like
$N^{1-\alpha/d}$ for $0\leq\alpha/d<1$ (long-range), and increases like $\ln
N$ for $\alpha/d=1$.

Let us emphasize at this point that both $(\beta_{q}\tilde{N})$ and $({\cal
H}_{N}/N\tilde{N})$ are intensive quantities. Indeed, let us illustrate these
facts by focusing on the following paradigmatic Gibbs thermodynamical energy:

$\displaystyle G(V,T,p,\mu,H,\dots)=U(V,T,p,\mu,H,\dots)$ $\displaystyle-
TS(V,T,p,\mu,H,\dots)+pV-\mu N(V,T,p,\mu,H,\dots)$ $\displaystyle-
HM(V,T,p,\mu,H,\dots)-\cdots\,,$ (1)

where $T,p,\mu,H$ are the temperature, pressure, chemical potential, external
magnetic field, and $U,S,V,N,M$ are the internal energy, entropy, volume,
number of particles (in turn proportional to the number of degrees of
freedom), magnetization. By dividing both sides by $N\tilde{N}$ we obtain

$\displaystyle\frac{G(V,T,p,\mu,H,\dots)}{N\tilde{N}}=\frac{U(V,T,p,\mu,H,\dots)}{N\tilde{N}}$
$\displaystyle-\frac{T}{\tilde{N}}\frac{S(V,T,p,\mu,H,\dots)}{N}+\frac{p}{\tilde{N}}\frac{V}{N}-\frac{\mu}{\tilde{N}}$
$\displaystyle-\frac{H}{\tilde{N}}\frac{M(V,T,p,\mu,H,\dots)}{N}-\cdots\,.$
(2)

It has been profusely verified in the literature (see Tsallisbook ;
TsallisCirto2013  and references therein) that all the quantities
$G/(N\tilde{N})$, $U/(N\tilde{N})$, $T/\tilde{N}$, $S/N$, $p/\tilde{N}$,
$V/N$, $\mu/\tilde{N}$, $H/\tilde{N}$, $M/N$ are thermodynamically intensive,
in the sense that, in the $N\to\infty$ limit, they all yield finite
quantities, thus preserving the Legendre structure of classical thermodynamics
for both short- and long-range interactions.

We then identify three classes of thermodynamical variables for all values of
$\alpha/d$, namely (i) those that are expected to always be extensive
($S,V,N,M,\ldots$), i.e., scaling with $N$, (ii) those that characterize the
external conditions under which the system is placed ($T,p,\mu,H,\ldots$),
scaling with $\tilde{N}$, and (iii) those corresponding to energies ($G,U$),
scaling with $N\tilde{N}$. For short-range interactions (i.e., $\alpha/d>1$),
these three classes collapse into the traditional two (intensive and
extensive) currently indicated in the textbooks of thermodynamics.

The desirable mathematical counterpart for such systems would of course be to
have a $q$-Large Deviation Theory ($q$-LDT) with a probability corresponding
to $n/N-1/2\geq z$ given by $P(N,z)\approx e_{q}^{-r_{q}(z)N}$, where the rate
function $r_{q}(z)$ would once again equal some relative nonadditive entropy
per particle defined through Tsallis1988

$S_{q}=k\frac{1-\sum_{i}p_{i}^{q}}{q-1}\;\;(q\in\mathbb{R})\,$ (3)

with $S_{1}=S_{BG}\equiv-k\sum_{i}p_{i}\ln p_{i}$, $k$ being a conventional
positive constant (hereafter taken to be $k=1$). A more precise notation for
$r_{q}(z)$ would be $r_{q_{r}}(z)$, since, as we shall verify here below,
there is no reason for being $q_{r}=q$; in fact, in all the nontrivial cases
that we are aware of, it appears to be $q_{r}\neq q$. However, for simplicity,
we shall maintain the notation $r_{q}$.

The quantity $r_{q}(z)N$ is expected to play a role similar to that of a total
system thermodynamic entropy which, as mentioned above, should always be
extensive, i.e., $\propto N$ ($N\gg 1)$. Naturally, in order to unify all the
above situations, we expect $q=f(Q)$, $f(Q)$ being a smooth function which
satisfies $f(1)=1$, thus recovering the usual LDT.

The above $q$-LDT scenario has already been numerically verified for a purely
probabilistic model with strong correlations RuizTsallis2012 ; Touchette2012 ;
RuizTsallis2013 , as well as for diverse physical models TirnakliTsallisAy2021
.

##  II Model and results

In the present paper, we focus on a scale-invariant probabilistic model
introduced in HanelThurnerTsallis2009  and based on the Laplace-de Finetti
theorem for exchangeable stochastic processes. The random variables are binary
(Ising-like) and can be either uncorrelated ($Q=1$) or strongly correlated
($Q>1$); each of them takes the values $0$ and $1$. A specific micro-state
with $n$ values $0$ and $(N-n)$ values $1$ corresponds to $r_{n}^{N}=1/2^{N}$
for $Q=1$ and to

$r_{n}^{N}=\frac{B\Bigl{(}\frac{3-Q}{2Q-2}+n,\frac{3-Q}{2Q-2}+N-n\Bigr{)}}{B\Bigl{(}\frac{3-Q}{2Q-2},\frac{3-Q}{2Q-2}\Bigr{)}}$
(4)

for long-tailed distributions ($1<Q<3$), where
$B(x,y)=\Gamma(x)\Gamma(y)/\Gamma(x+y)$ is the Euler Beta function. In this
model, $(Q-1)$ measures the strength of the global correlations, and varies
from zero to 2. As a physical analog we may think of a classical
$d$-dimensional spin-1/2 highly anisotropic XY ferromagnet with two-body
interactions decaying as $1/(distance)^{\alpha}$. This class of systems
approach the Ising ferromagnet and typically corresponds to magnets with two-
body highly anisotropic XY interactions in the absence of one-body terms, and
also to two-body fully isotropic XY interactions but having, in addition,
highly anisotropic one-body terms. Then $Q=1$ corresponds to $\alpha/d>1$,
whereas $Q>1$ is expected to mirror systems with $0<\alpha/d<1$. An extreme
such case is the $\alpha=0$ one, where all spins interact equally strongly
with all the other spins of the system. The $Q$-Gaussian distribution would
then possibly mirror a non-Maxwellian distribution of velocities. Such
distributions have been repeatedly observed in similarly complex systems,
e.g., long-range isotropic XY and Heisenberg ferromagnets as well as the long-
range Fermi-Pasta-Ulam $\beta$-model.

In the present model there are $N!/[n!(N-n)!]$ equivalent such micro-states
($n=0,1,2,\dots,N$). Consistently, we have

$\sum_{n=0}^{N}\frac{N!}{n!(N-n)!}r_{n}^{N}=1\,.$ (5)

The quantities $r_{n}^{N}$ defined in Eq. (4) satisfy (see
HanelThurnerTsallis2009 ) the so-called Leibnitz triangle rule, i.e.,

$r_{n}^{N}+r_{n+1}^{N}=r_{n}^{N-1}\;\;(\forall N,\,\forall n)\,.$ (6)

Such sets of probabilities are also known as scale-invariant (see also
RodriguezSchwammleTsallis2008 ; RodriguezTsallis2012 ; RodriguezTsallis2014 )
since their distribution at a given scale $N$ can always be obtained by
probabilistic marginalization of higher scales (corresponding to
$N+1,N+2,\dots$). This property is rather special indeed and should not be
confused with the so-called “rule of marginalization”, which is valid for any
distribution of probabilities. An example of well defined distribution of $N$
binary random variables which generically violates relation (6) can be seen in
TsallisGellMannSato2005 . The strength of the Leibnitz triangle rule (see, for
instance Polya1945  and related works) can be illustrated by the fact that,
under this remarkable hypothesis, the entire set
$\\{r_{n}^{N}\\}\,,\forall(N,n)$, can be generically and univocally recovered
by only providing, for instance, the set $\\{r_{0}^{N}\\}\,,\forall N$. In
order to provide a more complete view on this issue, let us clarify however
that if the system satisfies, as in the present model, exchangeability of the
random variables, then the rule of marginalization implies the Leibnitz
triangle rule.

The $N\to\infty$ attractor of this model for fixed $Q$ turns out
HanelThurnerTsallis2009  to precisely be a $Q$-Gaussian, which makes it an
interesting case for checking its LDT behavior.

Let us briefly review the algorithm of HanelThurnerTsallis2009  for
constructing the distributions which, in the $N\to\infty$ limit, yield exact
$Q$-Gaussians. We define

$u_{n}^{N}\equiv\frac{(n/N-1/2)}{\sqrt{(Q-1)(n/N)(1-n/N)}}$ (7)

and also

${\tilde{u}}_{n}^{N}\equiv\frac{u_{n}^{N}}{2\,\max_{n=1,2,\dots,N-1}\\{u_{n}^{N}\\}}$
(8)

where
$\max_{n=1,2,\dots,N-1}\\{u_{n}^{N}\\}=\frac{\frac{1}{2}-\frac{1}{N}}{\sqrt{(Q-1)(N-1)/N^{2}}}$.
We also define the discrete width

$du_{n}^{N}\equiv\frac{[(n/N)(1-n/N)]^{-3/2}}{4(N+1)\sqrt{Q-1}}$ (9)

from which it follows the (un-normalized) distribution

$F_{n}^{N}=(du_{n}^{N})^{-1}\frac{N!}{n!(N-n)!)}r_{n}^{N}\,,$ (10)

where $r_{n}^{N}$ is given by Eq. (4) and, after normalization, we have

${\tilde{F}}_{n}^{N}\equiv\frac{F_{n}^{N}}{\sum_{n=1}^{N-1}F_{n}^{N}}\,.$ (11)

The values $n=0$ and $n=N$ are excluded from the sum because these values map
to infinity in Eq. (9). ${\tilde{F}}_{n}^{N}$ is the distribution which, for
$N\to\infty$, is attracted by a $Q$-Gaussian. In Fig. 1a we have represented
the data $\\{(N{\tilde{u}}_{n}^{N},{\tilde{F}}_{n}^{N})\\}$, whereas in Fig.
1b the same data have been plotted with respect to ${\tilde{u}}$ so that the
distribution can be given in the region $[-1/2,1/2]$. Through the rescaling
$\bar{u}\equiv\sqrt{2N}{\tilde{u}}$, one can easily obtain the normalized
attractor
${\bar{F}}_{n}^{N}(\bar{u})=\sqrt{2N}{\tilde{F}}_{n}^{N}(\tilde{u})=e_{Q}^{-{\bar{u}}^{2}}$.
The discrete width in Eq. (9) is related to the nonequidistance observed
between the points in Fig. 1.

Now, in the LDT realm, we focus on the probability $P(N,z)\in[0,1]$ which is
defined as that whose values of $F_{n}^{N}$ correspond to $n/N>1/2+z$. More
precisely, it is the sum of all values whose $\tilde{u}>z$ (see shadowed areas
in Fig. 1b). It is clear that $P(N,0)=1/2$ and $P(N,1/2)=0$. In view of what
has been discussed above, we expect to numerically verify that

$P(N,z)=P_{0}(Q,z)\,e_{q}^{-r_{q}(Q,z)\,N}\,,$ (12)

with $P_{0}(Q,0)=1/2$, $P_{0}(Q,1/2)=0$, and $r_{q}(Q,0)=0$.

By optimally fitting, with respect to $(q,r_{q},P_{0})$, the data with Eq.
(12), we have heuristically found $q$ to be given by the simple composition of
the well known additive duality ($q\to 2-q$) on top of the multiplicative
duality ($q\to 1/q$) (see, for instance, TsallisGellMannSato2005 ;
MoyanoTsallisGellMann2006 ), namely

$q=2-\frac{1}{Q}\;\;\;(1\leq Q<3)\,,$ (13)

or, equivalently,

$\frac{1}{q-1}=\frac{1}{Q-1}+1\,.$ (14)

For uncorrelated binary variables, we have (from RuizTsallis2012  with
$z=x-1/2$) the following BG relative entropy (with regard to equal
probabilities)

$\displaystyle r_{1}(z)$ $\displaystyle=$ $\displaystyle\ln
2+\frac{1+2z}{2}\ln\frac{1+2z}{2}+\frac{1-2z}{2}\ln\frac{1-2z}{2}$ (15)
$\displaystyle\sim$ $\displaystyle 2z^{2}+\frac{4}{3}z^{4}\;\;\;(z\to 0)\,.$
(16)

hence $r_{1}(z)\in[0,\ln 2]$. More generally, for strongly correlated binary
variables, we have (from RuizTsallis2012 )

$\displaystyle r_{q}(z)$ $\displaystyle=$
$\displaystyle\frac{1}{q_{r}-1}\Bigl{\\{}\frac{1}{2}[(1+2z)^{q_{r}}+(1-2z)^{q_{r}}]-1\Bigr{\\}}$
(17) $\displaystyle\sim$ $\displaystyle
2q_{r}z^{2}+\frac{2}{3}(3-q_{r})(2-q_{r})q_{r}z^{4}\;\;(z\to 0),$ (18)

where we have used $p=1/2+z$ and $1-p=1/2-z$ in expression (3), adopting as a
reference distribution the equal probability case. It follows that
$r_{q}(z)\in[0,\ln_{2-q_{r}}2]$, and we verify that $r_{q}(z)$ recovers
$r_{1}(z)$ for $q_{r}\to 1$.

Figure 1:  ${\tilde{F}}_{n}^{N}$ distributions are given for some
representative values of $N$. Dashed lines are the corresponding $Q$-Gaussians
($\propto e_{Q}^{-2N{\tilde{u}}^{2}}$). The sums of the values of all points
equals unity. Notice a relevant point, namely that the abscissa values of
these points are not equidistant. (a) The distribution is represented as a
function of $N{\tilde{u}}\in[-N/2,N/2]$. (b) The distribution is represented
as a function of $\tilde{u}\in[-1/2,1/2]$. An arbitrary value $z\in[0,1/2]$ is
indicated as well.

Figure 2:  $P(N,z)$ for $Q=1.5$ in log-log (a), $q$-log (b) and ratio (c)
representations. Note that the only distribution which provides, in all
scales, straight lines in a $\ln_{q}x$ versus $x$ representation is the
$q$-exponential function. It is worthy to mention here that, in our
calculations, $z$ values lies in $[0.035,0.17]$. Calculations for $z$ even
larger overcome our present computational capacity. To numerically process the
generically chosen $z$ value, we attribute to ${\tilde{F}}_{n}^{N}$ a 5-point
cubic (polynomial) interpolation using five neighboring data points (we
specifically used interp1d class in scipy.interpolate of Python). We have
checked for all the values of $z$ but have illustrated in (c) with only one
representative example.

Figure 3:  (a) $q_{r}$ versus $Q$ is plotted. Dots are obtained through
optimization of the overall fitting of $P(N,z)$ with regard to $(P_{0},q_{r})$
for typical values of $Q$ and various values for $z\in[0.035,0.17]$. The
optimization procedure uses the scipy.optimize module in Python with method
SLSQP. The dashed line is the conjecture given in Eq. 19. For each value of
$z$, we obtain, through an optimized fitting, a value of $q_{r}$. For the
whole set of values of $z$, we then obtain the maximal and minimal values of
$q_{r}$, which determine the upper and lower error bars. The averaged values
are indicated by points. Inset: $P_{0}$ versus $z$ for all $Q$ values seen in
the figure; the dashed line is given by Eq. 20 with the illustrative value
$u=0.4$. (b) $r_{q}$ values calculated from Eq. 17 are plotted as a function
of $z^{2}$ for a representative $Q$ value. The shaded region indicates the
interval consistent with the error bars of $q_{r}$ for this $Q$ value.
Calculations for $z^{2}$ even larger overcome our present computational
capacity.

Through the optimized fitting, we heuristically found the following relation

$q_{r}=\frac{7}{10}+\frac{6}{10}\frac{1}{Q-1}\;\;\;(1<Q<3).$ (19)

This result was obtained by making $z$ to typically vary up to $0.17$ for
diverse values of $Q$, which guarantees the verification of the dominant term
in Eq. (18). If we could numerically check up to $z=1/2$, we could guarantee
the full expression (17), but this remains out of our present computational
capacity.

The numerical determination of $P_{0}(Q,z)$ is much harder than that of
$(q,q_{r})$. However, just as a simple indication, we have compared the
numerical data with

$P_{0}(Q,z)=1/4-az^{u}+a(1/2-z)^{u}\;\;\;(0\leq z\leq 1/2))\,,$ (20)

where $a=2^{u}/4$ in order to satisfy the conditions $P_{0}(Q,0)=1/2$,
$P_{0}(Q,1/2)=0$. This particular form was heuristically proposed as a simple
illustration; it satisfies $P_{0}(Q,z)=P_{0}(Q,1/2-z)$, thus exhibiting an
inflexion point at $z=1/4$ which was suggested by numerical exploration (see
the Inset of Fig. 3 (a)). Let us however emphasize here that the precise
numerical values of $P_{0}$ (as well as its unknown exact analytical
expression) have no particular relevance. They play in fact a rather minor
role in the conjecture (12), in complete similarity with the corresponding
pre-factor in the standard LDT.

As can be seen in Figs. 2 and 3, strong numerical evidence supports conjecture
in Eq. (12) with relations (13), (17), (19) and (20). Let us stress that the
rate function $r_{q}$ yields an extensive (relative) total entropy $r_{q}N$
for all values of $Q\in[1,3)$, as mandated by the Legendre structure of
thermodynamics. However, the corresponding entropic index is not $q(Q)$ but a
different one, namely $q_{r}(Q)$. It remains as an open question whether this
is the generic case, or rather an exception associated to the present specific
model. The fact that three different entropic indices emerge, namely
$(Q,q,q_{r})$, is not particularly surprising, given the Moebius algebra which
characterizes $q$-statistics (see TsallisGellMannSato2005 , GazeauTsallis2019
and references therein). Another point which deserves emphasis concerns the
fact that infinitely many classes of distributions exist which asymptotically
are power-laws, for $N>>1$. However, only one of these infinite functional
forms is the $q$-exponential. All these power-laws definitively differ in the
non-asymptotic region, i.e., for relatively low values of $N$, say $N\sim
50,100$. This is the reason for which special numerical attention has been
here devoted to that region of $N$. The evidence that has been achieved
strongly points in favor of precisely the $q$-exponential form.

##  III Final remarks

At this point, let us conclude by reminding that our aim is to approach,
within a more general context, the fingerprints of Boltzmann-Gibbs statistical
mechanics, namely the Maxwellian distribution of velocities and the BG
exponential weight for the energies. Indeed, in the realm of $q$-statistics
based on nonadditive entropies, a $Q$-Gaussian distribution emerges for the
velocities and a $q$-exponential weight emerges for the energies, with $Q\geq
q\geq 1$, the equalities holding precisely for the BG theory. These
generalizations should respectively mirror corresponding generalizations of
the classical Central Limit Theorem and the Large Deviation Theory. This
scenario has been successfully verified for a purely probabilistic model (see
RuizTsallis2013  and references therein), as well as for some simple physical
models TirnakliTsallisAy2021 .

Finally, a comparison with previous results might be helpful. In Ref.
RuizTsallis2012 , a completely different probabilistic model was focused on.
In Ref. TirnakliTsallisAy2021  we have focused on four existing, physically
motivated, models. The model that we focus on here is the one introduced in
Ref. HanelThurnerTsallis2009 , where no reference at all exists to a possible
$q$-large deviation theory. In our present manuscript, it has been possible to
numerically discuss (with satisfactory precision in some cases) (i) The
$(Q,z)$-dependence of the pre-factor $P_{0}(Q,z)$ (this is the first time such
a pre-factor is focused on in the literature of complex systems); (ii) The
possible identification of the rate function $r_{q}(z)$ with a nonadditive
relative entropy whose index is $q_{r}$, definitively different from the index
$q$ (a possibility that has never been handled before) [a definitive numerical
identification of $r_{q}(z)$ with the $q_{r}$-entropy for the full range of
$z$ was not possible because our present computational capacity does not allow
us to increase $z$ up to $z=1/2$; still, the discussion for $z\simeq 0$ has
never been undertaken before]; (iii) The $Q$-dependence $q_{r}(Q)$ as
indicated in Eq. (19), which includes an unexpected singularity at $Q=1$ (a
feature coming from the specificity of the model introduced in Ref.
HanelThurnerTsallis2009 , which, as said above, neatly differs from the model
discussed in Ref. RuizTsallis2012 ). This $q_{r}(Q)$ dependence only became
accessible due to the mathematical fact that all the necessary information is
already available in the first asymptotic term, namely in the quadratic term
of $r_{q}(z)$ as a function of $z$, as exhibited in Eq. (18). To the best of
our knowledge, the above three points have never before been simultaneously
attained for any nontrivial model, and certainly not in Refs. RuizTsallis2012
and TirnakliTsallisAy2021 . Last but not least, in all four models focused on
in TirnakliTsallisAy2021 , the values of $(Q,q)$ are numerical ones, whereas
in the present paper we have obtained analytical expressions for arbitrary
real $Q>1$.

Along the lines of the promising results presented herein, analytical
approaches (or very high precision numerical approaches) would naturally be
very welcome, either for the specific models studied here and elsewhere
RuizTsallis2012 ; TirnakliTsallisAy2021 , or in the ambitious form of a
$q$-generalized theorem for large deviations based on say a $Q$-generalized
central limit theorem for an important class of strongly correlated random
variables, which frequently emerges in physics, geophysics, astrophysics,
economics, among other areas. We hope that the present numerical indications
may stimulate $q$-generalizations of theorems such as the Gärtner-Ellis one
Gartner1977 ; Ellis1984 . In fact, a fascinating research field seemingly
exists at the crossroad of the large deviation theory Lanford1973 ; Ellis1985
; Oono1989 ; Ellis1995 ; Ellis1999  and nonadditive entropies such as $S_{q}$.
Since $S_{q}$ stands at the foundations of nonextensive statistical mechanics,
such contributions would naturally enable a deeper understanding of the
mathematical structure of this current generalization of the Boltzmann-Gibbs
theory.

In $q$-statistics, new concepts such as the $q$-triplet appear. This was
experimentally verified for the first time in BurlagaVinas2005  and, since
then, in a plethora of other systems. No such complex structure exists in the
Boltzmann-Gibbs theory, where, generically, a unique value for $q$ is
admissible, namely $q=1$. The various values of the indices
$(Q,q(Q),q_{r}(Q))$ that emerge in the present model, together with one more
HanelThurnerTsallis2009  ($q_{entropy}=1$, $\forall Q$, corresponding to the
index of the entropic functional $S_{q_{entropy}}(N)$ which, for the present
exchangeable stochastic system, is thermodynamically extensive, i.e.,
$S_{q_{entropy}}(N)\propto N$ for $N\to\infty$, $\forall Q$) are not yet fully
elucidated. Some preliminary understanding is available in
TsallisGellMannSato2005 ; GazeauTsallis2019 , but this rich issue still
remains as a nontrivial and intriguing open problem.

## Acknowledgments

The numerical calculations reported in this paper were partially performed at
TUBITAK ULAKBIM, High Performance and Grid Computing Center (TRUBA resources).
We acknowledge interesting remarks from E.M.F. Curado and R. Hanel, as well as
partial financial support from CNPq and Faperj (Brazilian agencies). U.T. is a
member of the Science Academy, Bilim Akademisi, Turkey and acknowledges
partial support from TUBITAK (Turkish Agency) under the Research Project
number 121F269.

## References

  * (1) F. Reif, Fundamentals of Statistical and Thermal Physics (Waveland Press, Long Grove, 2008).
  * (2) R.K. Pathria and P.D. Beale, Statistical Mechanics (Academic Press, New York, 2011).
  * (3) P. Billingsley, Convergence of Probability Measures (Wiley, New York, 1968).
  * (4) N.G. van Kampen, Stochastic Processes in Physics and Chemistry (North-Holland, Amsterdam, 1981).
  * (5) R.S. Ellis, Entropy, Large Deviations and Statistical Mechanics (Springer, Berlin, 1985).
  * (6) C. Giardina, J. Kurchan and L. Peliti, Direct evaluation of large-deviation functions, Phys. Rev. Lett. 96, 120603 (2006).
  * (7) F. den Hollander, Large Deviations (American Mathematical Society, USA, 2008).
  * (8) V. Lecomte and J. Tailleur, A numerical approach to large deviations in continuous-time, J. Stat. Mech., P03004 (2007).
  * (9) H. Touchette, The large deviation approach to statistical mechanics, Phys. Rep. 478, 1-69 (2009).
  * (10) T. Nemoto, E. G. Hidalgo and V. Lecomte, Finite-time and finite-size scalings in the evaluation of large deviation functions: Analytical study using a birth-death process, Phys. Rev. E 95, 012102 (2017).
  * (11) C. Tsallis, Possible generalization of Boltzmann-Gibbs statistics, J. Stat. Phys. 52, 479-487 (1988).
  * (12) R.J.V. dos Santos, Generalization of Shannon’s theorem for Tsallis entropy, J. Math. Phys. 38, 4104 (1997).
  * (13) S. Abe, Axioms and uniqueness theorem for Tsallis entropy, Phys. Lett. A 271, 74 (2000).
  * (14) G. Wilk and Z. Wlodarczyk, Interpretation of the nonextensivity parameter $q$ in some applications of Tsallis statistics and Levy distributions, Phys. Rev. Lett. 84, 2770 (2000).
  * (15) C. Beck, Dynamical foundations of nonextensive statistical mechanics, Phys. Rev. Lett. 87, 180601 (2001).
  * (16) C. Beck and E.G.D. Cohen, Superstatistics, Physica A 322, 267 (2003).
  * (17) R. Hanel, S. Thurner and M. Gell-Mann, How multiplicity determines entropy: derivation of the maximum entropy principle for complex systems, PNAS 111, 6905 (2014).
  * (18) A. Enciso and P. Tempesta, Uniqueness and characterization theorems for generalized entropies, J. Stat. Mech. 123101 (2017).
  * (19) P. Jizba, J. Korbel, H. Lavicka, M. Proks, V. Svoboda and C. Beck, Transitions between superstatistical regimes: validity, breakdown and applications, Physica A 493, 29 (2018).
  * (20) P. Jizba and J. Korbel, Maximum entropy principle in statistical inference: case for non-Shannonian entropies, Phys. Rev. Lett. 122, 120601 (2019).
  * (21) H. Suyari, H. Matsuzoe and A.M. Scarfone, Advantages of q-logarithm representation over q-exponential representation from the sense of scale and shift on nonlinear systems, Eur. Phys. J. Special Topics 229, 773-785 (2020).
  * (22) C. Tsallis, Introduction to Nonextensive Statistical Mechanics–Approaching a Complex World (Springer, New York, 2009).
  * (23) C. Anteneodo and C. Tsallis, Breakdown of the exponential sensitivity to the initial conditions: Role of the range of the interaction, Phys. Rev. Lett. 80, 5313 (1998).
  * (24) L.J.L. Cirto, V.R.V. Assis and C. Tsallis, Influence of the interaction range on the thermostatistics of a classical many-body system, Physica A 393, 286-296 (2014).
  * (25) L.J.L. Cirto, A. Rodriguez, F.D. Nobre and C. Tsallis, Validity and failure of the Boltzmann weight, EPL 123, 30003 (2018).
  * (26) A, Rodriguez, F.D. Nobre and C. Tsallis, $d$-Dimensional classical Heisenberg model with arbitrarily-ranged interactions: Lyapunov exponents and distributions of momenta and energies, Entropy 21, 31 (2019).
  * (27) H. Christodoulidi, C. Tsallis and T. Bountis, Fermi-Pasta-Ulam model with long-range interactions: Dynamics and thermostatistics, EPL 108, 40006 (2014).
  * (28) H. Christodoulidi, T. Bountis, C. Tsallis and L. Drossos, Dynamics and Statistics of the Fermi–Pasta–Ulam $\beta$–model with different ranges of particle interactions, JSTAT, 123206 (2016).
  * (29) D. Bagchi and C. Tsallis, Sensitivity to initial conditions of $d$-dimensional long-range-interacting quartic Fermi-Pasta-Ulam model: Universal scaling, Phys. Rev. E 93, 062213 (2016).
  * (30) D. Bagchi and C. Tsallis, Long-ranged Fermi-Pasta-Ulam systems in thermal contact: Crossover from q-statistics to Boltzmann-Gibbs statistics, Phys. Lett. A 381, 1123 (2017).
  * (31) S. Umarov, C. Tsallis, S. Steinberg, On a $q$-central limit theorem consistent with nonextensive statistical mechanics, Milan J. Math. 76, 307 (2008).
  * (32) S. Umarov, C. Tsallis, M. Gell-Mann, S. Steinberg, Generalization of symmetric $\alpha$-stable Lévy distributions for $q>1$, J. Math. Phys. 51, 033502 (2010).
  * (33) M.G. Hahn, X.X. Jiang, S. Umarov, On $q$-Gaussians and exchangeability, J. Phys. A 43, 165208 (2010).
  * (34) J.S. Andrade Jr., G.F.T. da Silva, A.A. Moreira, F.D. Nobre and E.M.F. Curado, Thermostatistics of overdamped motion of interacting particles, Phys. Rev. Lett. 105, 260601 (2010).
  * (35) C. Tsallis and L.J.L. Cirto, Black hole thermodynamical entropy, Eur. Phys. J. C 73, 2487 (2013).
  * (36) G. Ruiz and C. Tsallis, Towards a large deviation theory for strongly correlated systems, Phys. Lett. A 376, 2451 (2012).
  * (37) H. Touchette, Comment on “Towards a large deviation theory for strongly correlated systems”, Phys. Lett. A 377, 436 (2013).
  * (38) G. Ruiz and C. Tsallis, Reply to Comment on “Towards a large deviation theory for strongly correlated systems”, Phys. Lett. A 377, 49 (2013).
  * (39) U. Tirnakli, C. Tsallis and N. Ay, Approaching a large deviation theory for complex systems, Nonlinear Dynamics (2021), https://doi.org/10.1007/s11071-021-06904-3.
  * (40) R. Hanel, S. Thurner and C. Tsallis, Limit distributions of scale-invariant probabilistic models of correlated random variables with the $q$-Gaussian as an explicit example, Eur. Phys. J. B 72, 263 (2009).
  * (41) A. Rodriguez, V. Schwammle and C. Tsallis, Strictly and asymptotically scale-invariant probabilistic models of $N$ correlated binary random variables having q-Gaussians as $N\to\infty$ limiting distributions, J. Stat. Mech. P09006 (2008).
  * (42) A. Rodriguez and C. Tsallis, A dimension scale-invariant probabilistic model based on Leibniz-like pyramids, J. Math. Phys. 53, 023302 (2012).
  * (43) A. Rodriguez and C. Tsallis, Connection between Dirichlet distributions and a scale-invariant probabilistic model based on Leibniz-like pyramids, J. Stat. Mech. P12027 (2014).
  * (44) C. Tsallis, M. Gell-Mann and Y. Sato, Asymptotically scale-invariant occupancy of phase space makes the entropy $S_{q}$ extensive, Proc. Natl. Acad. Sc. USA 102, 15377 (2005).
  * (45) G. Pólya, How to Solve It, (Princeton University Press,1945).
  * (46) L.G. Moyano, C. Tsallis and M. Gell-Mann, Numerical indications of a $q$-generalised central limit theorem, Europhys. Lett. 73, 813 (2006).
  * (47) J.-P. Gazeau and C. Tsallis, Mobius transforms, cycles and $q$-triplets in statistical mechanics, Entropy 21, 1155 (2019).
  * (48) J. Gärtner, On large deviations from an invariant measure, Teor. Verojatnost. i Primenen. 22, 27-42 (1977).
  * (49) R.S. Ellis, Large deviations for a general class of random vectors, Ann. Probab. 12, 1-12 (1984).
  * (50) O.E. Lanford, Entropy and equilibrium states in classical statistical mechanics, in Statistical Mechanics and Mathematical Problems, A. Lenard, Ed., Lecture Notes in Physics 20, 1-113 (Springer, Berlin, 1973).
  * (51) R.S. Ellis, Entropy, Large Deviations, and Statistical Mechanics (Springer, New York, 1985).
  * (52) Y. Oono, Large deviation and statistical physics, Progr. Theoret. Phys. Suppl. 99,165-205 (1989).
  * (53) R.S. Ellis, An overview of the theory of large deviations and applications to statistical mechanics, Scand. Actuar. J. 1, 97-142 (1995).
  * (54) R.S. Ellis, The theory of large deviations: From Boltzmann’s 1877 calculation to equilibrium macrostates in 2D turbulence, Physica D 133 106-136 (1999).
  * (55) L.F. Burlaga and A.F. Vinas, Triangle for the entropic index $q$ of non-extensive statistical mechanic observed by Voyager 1 in the distant heliosphere, Physica A 356, 375 (2005).