#  Supernovae Time Profiles as a Probe of New Physics at Neutrino Telescopes

Jeff Lazar<EMAIL_ADDRESS>Department of Physics & Laboratory for
Particle Physics and Cosmology, Harvard University, Cambridge, MA 02138, USA
Department of Physics & Wisconsin IceCube Particle Astrophysics Center,
University of Wisconsin-Madison, Madison, WI 53706, USA  Ying-Ying Li
<EMAIL_ADDRESS>Peng Huanwu Center for Fundamental Theory, Hefei, Anhui
230026, China  Interdisciplinary Center for Theoretical Study, University of
Science and Technology of China, Hefei, Anhui 230026, China  Carlos A.
Argüelles<EMAIL_ADDRESS>Department of Physics & Laboratory for
Particle Physics and Cosmology, Harvard University, Cambridge, MA 02138, USA
Vedran Brdar<EMAIL_ADDRESS>Department of Physics, Oklahoma State
University, Stillwater, OK, 74078, USA

###### Abstract

Neutrino telescopes, including IceCube, can detect galactic supernova events
by observing the collective rise in photomultiplier count rates with a sub-
second time resolution. Leveraging precise timing, we demonstrate the ability
of neutrino telescopes to explore new weakly coupled states emitted from
supernovae and subsequently decaying to neutrinos. Our approach utilizes
publicly available packages, ASTERIA and SNEWPY, for simulating detector
responses and parametrizing neutrino fluxes originating from Standard Model
and new physics. We present results for two beyond the Standard Model
scenarios and introduce the tool developed for testing a diverse range of new
physics models.

††preprint: USTC-ICTS/PCFT-24-07

Introduction—The Standard Model (SM) is a remarkable but incomplete theory.
Puzzles such as the non-vanishing neutrino mass, the origin of observed
matter-antimatter asymmetry, and the nature of dark matter, among others, seek
explanations beyond the SM (BSM). Since there are also no firmly established
clues on the energy scale at which the new particles should appear, numerous
BSM searches across a wide range of scales have been performed: from collider
searches at high energies [1] to studies of cosmic microwave background at
temperatures near absolute zero [2].

Supernovae (SNe) release most of their energy in neutrinos, offering unique
opportunities to test BSM physics via neutrino interactions. The duration of
the burst [3, 4, 5] and inferred total energy of neutrinos [6, 7, 8] from SN
1987A—the only SN observed in neutrinos—have already given tentative
constraints on the energies that could be carried by light BSM particles
produced in the interior of a star. Indeed, many different BSM scenarios have
been constrained in this way, see e.g. Refs. [9, 10, 11, 12, 13, 14, 15, 16,
17, 18]. If, additionally, BSM particles emitted from SNe can decay to
neutrinos en route to Earth, even stronger limits can be set, as shown, for
instance, in Ref. [19] for Majoron-like bosons and in Ref. [20] for a
realization featuring a neutrino magnetic moment portal. In contrast to SM
neutrinos, BSM states produced in the SN core can freely stream out without
further interactions. Consequently, they exit with higher energies than
neutrinos, typically around $\mathcal{O}(100)$ MeV. Provided such particles
decay to neutrinos, strong constraints can be set from the fact that neutrinos
of $\mathcal{O}(100)$ MeV were not recorded from SN 1987A [19, 20]. Such
strategies were utilized for neutrino experiments operating during the SN
1987A event. Sensitivity projections for DUNE [21] and Hyper-Kamiokande [22]
in light of anticipated galactic SN event were also performed.

Figure 1: The decay geometry for Majorons produced in a SN. Both the
additional distance travelled relative to the direct line of sight and the
potential non-relativistic speed of the Majoron can induce a delayed signal
relative to the SM signal.

The IceCube Neutrino Observatory has pioneered the detection of neutrinos at
$\mathcal{O}$(PeV) energies [23, 24] and has identified several specific
astrophysical neutrino sources [24, 25, 26]. IceCube is also expected to be
able to detect the next galactic SN [27]. In fact, a search for MeV neutrinos
from optically obscured galactic SNe was recently performed in [28] as well as
the search for temporal correlation of MeV neutrino events with fast radio
bursts [29]. In both cases, the observable is a collective rise in all
photomultiplier rates on top of the background noise in a certain time window,
and it turns out that IceCube is sensitive to intervals as short as
$\mathcal{O}(0.01)$ seconds. In this Letter, we will demonstrate that such
precise time resolution will play a crucial role in testing new physics from
SNe.

BSM states can give rise to off-time signals in two ways. First, if the BSM
particles are non-relativistic, they will produce a delayed signal. However,
it is also possible that new states are light, implying that they may even
exit the SN at earlier times than certain flavors of SM neutrinos. Such an
early BSM signal, which to the best of our knowledge has not been previously
considered, can, for instance, be induced by light BSM states produced from
$\nu_{e}$ scattering inside the SN; they would then subsequently decay to
$\bar{\nu}_{e}$, which chiefly induces the signal at IceCube. This would
happen at early stages of an SN, around the neutralization burst, when
$\nu_{e}$ are present, and $\bar{\nu}_{e}$ have not been produced yet. Such
novel timing patterns generically exist in models with new weakly interacting
particles, such as the aforementioned Majoron model [19] and the neutrino
magnetic moment portal [18, 30].

In this work, we will show that, by using the timing measurements, neutrino
telescopes will be able to provide powerful constraints on the BSM parameter
space in association with forthcoming galactic SN events.

The Model—First, we consider the Majoron model, studied in Ref. [19], where SN
1987A constraints were derived. The relevant part of the Lagrangian reads

$\displaystyle\mathcal{L}\supset-
g_{\alpha\beta}\nu_{\alpha}\nu_{\beta}\phi+\text{h.c.}-m_{\phi}\phi\phi^{*}\,,$
(1)

where $\phi$ is the (pseudo)scalar, $m_{\phi}$ is its mass, and $g$
parametrizes interaction strength between $\phi$ and neutrinos. We assume
flavor universal interaction and hence $g_{\alpha\beta}\equiv
g_{\phi}\,\delta_{\alpha\beta}$. Majorons are produced from (anti)neutrino
coalescence in the star and subsequently decay to a pair of (anti)neutrinos,
resulting in BSM (anti)neutrino flux emitted from the star. The total decay
width is given by $\Gamma_{\phi}=3g^{2}m_{\phi}/16\pi$.

We will set SN bounce time $t_{\rm bounce}$ as $t=0$. Using data from the
simulation of $8.8M_{\odot}$ progenitor star [8] and including effects of
neutrino oscillations assuming normal mass ordering [31], we calculated
standard neutrino fluxes at Earth. Following Ref. [19], we also calculated the
flux of emitted Majorons. At $t\leq 0.05\,\mathrm{s}$, Majorons are mainly
produced through $\nu_{e}$ and $\nu_{x}$ coalescence where $\nu_{x}$ includes
(anti)neutrinos of $\tau$ and $\mu$ flavor. After $0.05\,\mathrm{s}$, the flux
of Majorons decreases with the total flux of neutrinos of all flavors. Thus,
we found that the Majoron flux peaks around $t\sim 0.05\,\mathrm{s}$. We
further cross-checked the agreement with Ref. [19] by reproducing their bound
from SN energy loss.

Very weakly coupled Majorons immediately stream out after being produced
inside the star. On their way to Earth, Majorons of energy $E_{\phi}$ will
travel at the speed of $\beta$ for a distance of $L_{1}$, then decay to
(anti)neutrinos of energy $E_{\nu}$ at an emission angle
$\cos\alpha=(2E_{\phi}E_{\nu}-m^{2}_{\phi})/(2E_{\nu}E_{\phi}\beta)$. Daughter
(anti)neutrinos will travel for a distance $L_{2}$ and reach Earth at the
angle $\theta$ satisfying $D_{\rm SN}\sin\theta=L_{1}\sin\alpha$ for a SN that
is $D_{\rm SN}$ away; see the decay geometry in Fig. 1. The time delay,
$\delta t$, relative to a relativistic particle moving along the line-of-sight
path from a SN to reach the Earth is given by [32]

$\displaystyle\delta
t=\left(\frac{L_{1}}{\beta}-L_{1}\right)+(L_{1}+L_{2}-D_{\rm SN})\,,$ (2)

where we have split the contribution to $\delta t$ into two parts: the first
part comes from the slow-moving Majorons before their decay, and the second
part is due to the detour from the line-of-sight path at the given emission
angle $\alpha$.

Figure 2: Flux of $\bar{\nu}_{e}$ from SM production and two Majoron
hypotheses at Earth. We choose the parameters $g_{\phi}=10^{-10.2}$ and
$g_{\phi}=10^{-11.8}$ for $m_{\phi}=1\,\mathrm{MeV}$ and
$m_{\phi}=200\,\mathrm{MeV}$, respectively.

We will consider a SN event that happens in the galaxy at a distance $D_{\rm
SN}=10\,\mathrm{kpc}$, which is not unlikely [33, 34]. To get the
characteristic value of $\delta t$, we can take $L_{1}=L_{\phi}$ with
$L_{\phi}=(E_{\phi}/m_{\phi})\Gamma^{-1}_{\phi}\beta$ being the decay length
of $\phi$. As for the parameter regions we considered, $(D_{\rm SN}-D_{\rm
SN}\cos\theta)\ll 10^{-3}\,\mathrm{s}$, $\delta t$ can be approximated as

$\displaystyle\delta t$ $\displaystyle\simeq
L_{\phi}\left(\frac{1}{\beta}-1\right)+L_{\phi}(1-\cos\alpha)=\frac{8\pi}{3E_{\nu}g_{\phi}^{2}}\,,$
(3)

which is roughly independent of Majoron energy $E_{\phi}$.

Figure 3: Timing profile of hits and test statistic as a function of the
timing window. In the left panel, we show the hit rate from the SM neutrino
flux which peaks at around 0.1 s, coinciding with the peak emission of
$\bar{\nu}_{e}$ from the SN. The shaded regions correspond to the timing
windows that give peak sensitivity to the particular model. These are the same
benchmarks in the parameter space as in Fig. 2. In the right panel, we show
the test statistic as a function of the starting time and length of the window
for the heavy Majoron benchmark ($m_{\phi}=200$ MeV). The relatively broad
distribution of high test-statistic values shows that this analysis technique
does not require extremely precise reconstruction of the
$t_{\mathrm{bounce}}$.

We can obtain the flux of daughter $\nu$ by considering decays that occur only
at $L_{1}\geq R^{\nu}_{\rm SN}$. As the smallest decay distance beyond which
the daughter neutrino can escape the explosion unperturbed is roughly the
radius of the core, we take $R^{\nu}_{\rm SN}=30\,\mathrm{km}$. We will also
limit our analysis to a maximal time window of $100\,\mathrm{s}$ after SN
bounce. This time window is larger than typical values of time delay arising
in the model.

As $\bar{\nu}_{e}$ contribute chiefly to the signals at IceCube via inverse
beta decay (IBD), we show in Fig. 2 the differential fluxes of daughter
$\bar{\nu}_{e}$ at different times for two benchmark points in the Majoron
model, in comparison to the SM $\bar{\nu}_{e}$ flux. We expect the observed
signal timing distribution at IceCube to follow these patterns of
$\bar{\nu}_{e}$ flux.

From Fig. 2, we can observe that the energy of $\bar{\nu}_{e}$ from Majoron
decay extends to $100\,\mathrm{MeV}$ and above for both benchmarks, as
expected. Consider the light Majoron case in Fig. 2. Majorons, generating
neutrinos with $E_{\bar{\nu}_{e}}\geq 150\,\mathrm{MeV}$, are nearly
relativistic with the emission angle $\alpha$ close to zero. The time delay
can be estimated from the first term in Eq. 3 given by $\delta t\sim
L_{\phi}(1/\beta-1)\lesssim 0.01\,\mathrm{s}$. With such a negligible time
delay, the time dependence of the $\bar{\nu}_{e}$ flux is mainly inherited
from that of the Majorons upon production, with larger fluxes at $t\leq
0.05\,\mathrm{s}$. Notice that such a $\bar{\nu}_{e}$ flux would arrive at the
detector earlier than the peak of the SM $\bar{\nu}_{e}$ flux, which is around
$t\simeq 0.1\,\mathrm{s}$ (see Fig. 2), and it potentially leads to early
signals at IceCube. Following Eq. 3, the time delay for $\bar{\nu}_{e}$ with
energy $E_{\bar{\nu}_{e}}\sim 10\,\mathrm{MeV}$ is typically $\gtrsim
0.1\,\mathrm{s}$. These $\bar{\nu}_{e}$ could be produced from Majorons with
low energy of $\mathcal{O}(10\,\mathrm{MeV})$ or with high energies. For the
low-energy Majoron case, the time delay is mainly due to slow-moving Majorons,
while for the high-energy case, it is mainly from the detour. Such sizable
time delays will shift the $\bar{\nu}_{e}$ peak flux to later times, in
comparison to that of their mother Majorons, as exhibited in Fig. 2.

For the heavy Majoron case, the resulting $\bar{\nu}_{e}$ flux is larger at
later times, manifesting the time delay from slowly moving Majorons. For the
heavy Majoron benchmark in Fig. 2, we can estimate that $\delta t\gtrsim
10\,\mathrm{s}$ for $E_{\bar{\nu}_{e}}\sim 200\,\mathrm{MeV}$ which is delayed
compared to the peak of the SM $\bar{\nu}_{e}$ flux. Compared to the light
Majoron case, such a large time delay is a combination of a more slowly moving
Majoron and a larger detour due to its large emission angle.

As $L_{1}\lesssim L_{\phi}$, the time spread of the flux is roughly given by
the aforementioned characteristic value of time delay, which is
$0.1\,\mathrm{s}$ ($\gtrsim 10\,\mathrm{s}$) for the light (heavy) benchmark
case.

Detector Response and Statistical Treatment—The IceCube Neutrino Observatory
comprises 5,160 light-detecting digital optical modules (DOMs) buried in a
cubic kilometer of the deep, transparent Antarctic ice sheet [35]. The DOMs
are arranged on 86 strings of 60 DOMs each, with an interstring distance of
125 m (70 m for the eight centermost strings). These DOMs detect the photons
emitted by the charged byproducts of a neutrino interaction in or near the
instrumented volume. This enables IceCube to resolve individual neutrinos with
energies ranging between a few GeV and a few PeV. The neutrinos produced by
SNe are far below this threshold and cannot be individually resolved; however,
the immense number of neutrinos produced in a SN increases the single-photon
rate of the detector. This dramatic increase in the rate of photons can be
distinguished from the background caused by dark noise and radioactive
activity in the DOM glass to enable the detection of galactic SN events [36].

New physics will affect the development of a SN and distort the temporal
structure of the photon signal seen in the IceCube detector. Thus, one may
look for excess events in certain time intervals as evidence of this new
physics. In this work, we simulate the light curve produced by a standard SN
with a progenitor mass of $8.8~{}M_{\odot}$ and those produced by different
BSM scenarios mentioned previously; see the left-hand panel of Fig. 3 for
light curves of the same two benchmarks from Fig. 2.

Since more than 93% of the photons detected in IceCube originate in IBD
interactions [37], the SM-only event rate peaks at around 0.1 s. BSM
scenarios, on the other hand, can give rise to an early or delayed signal
depending on the particulars of the model. In essence, Majorons can be
generated before the peak of the SM $\bar{\nu}_{e}$ flux, and as light
Majorons can travel nearly relativistically, they can escape the SNe with
negligible time delay and then decay to an all-flavor flux of neutrinos,
producing early $\bar{\nu}_{e}$ fluxes before the SM one. The higher-mass
Majorons, on the other hand, acquire a time delay from traveling sub-luminally
and at a large detour, thus producing a later signal with sizable time delay.

We use the ASTERIA [38] package to simulate the detector response to the
different SN scenarios and the thermal and radioactive noise in the DOMs. This
package simulates light yields from coherent $\nu_{e}$ scattering off
electrons in the ice, IBD of $\bar{\nu}_{e}$ on nuclei in the ice, and
charged- and neutral-current interactions for $\nu_{\alpha}$. In addition to
photons produced in neutrino interactions, it also simulates photons from
thermal noise in the IceCube photomultipliers (PMTs) and from radioactive
decay in the pressure glass housing. In order to interface with ASTERIA, we
use the SNEWPY [39] package and, in particular, the ParametrizedFlux object to
parametrize the flux with the “pinching” factor [40, 39].

Since each DOM only sees photons originating from within a sphere of radius
$\sim 5.2$ m [37], each module can be treated independently if the intermodule
spacing is $\gtrsim 10.4$ m. This criteria will be met in IceCube-Gen2 [41],
and thus, we should expect the total number of hits to scale like the
effective photocathode area, i.e.:

$A_{\mathrm{eff}}^{\mathrm{PC}}=\sum_{i}\varepsilon_{i}A_{i}^{\mathrm{PC}},$

where $\varepsilon_{i}$ and $A_{i}^{\mathrm{PC}}$ are the quantum efficiency
and photocathode area of a PMT and the sum runs over all PMTs in the detector.
Using this method, we rescale the results from the simulation from IceCube to
IceCube-Gen2 in order to estimate IceCube-Gen2’s sensitivity.

With the number of hits in the detector from noise, SM-only scenarios, and BSM
scenarios as a function of time, we can then quantify the probability of
seeing a certain number of hits in a given time window with

$-2\Delta\mathrm{LLH}=2\left[N_{\mathrm{exp.}}-N_{\mathrm{obs.}}+N_{\mathrm{obs.}}\log\left(\frac{N_{\mathrm{obs.}}}{N_{\mathrm{exp.}}}\right)\right],$
(4)

where $N_{\mathrm{obs.}}$ is the number of photons seen in the detector in the
given time and $N_{\mathrm{exp.}}$ is the number of photons expected from a
particular BSM hypothesis. For each BSM hypothesis, we select the time range
that maximizes the test statistic and use that maximal test statistic value.
The optimal ranges for each physics hypothesis are shown by the shaded regions
in the left panel of Fig. 3, which is around $0.1\,\mathrm{s}$
($10\,\mathrm{s}$) for the light (heavy) benchmark case, consistent with its
time spread of the produced $\bar{\nu}_{e}$ flux.

In the right-hand panel of Fig. 3, we show the test statistic as a function of
the time at which the time range begins, $t_{\mathrm{start}}$, and the
duration of the time range, $\Delta t$, for an example BSM hypothesis. The
distribution of test-statistic values is not sharply peaked in time. Thus,
this analysis retains much of its sensitivity even with errors on the
reconstruction of
$t_{\mathrm{bounce}}\sim\mathcal{O}\left(10^{-3}~{}\mathrm{s}\right)$ [42].

Results and Discussions—

Figure 4: Exclusion sensitivities for the Majoron case. The lines on this plot
show the exclusion sensitivity for IceCube, IceCube Gen-2, Hyper-Kamiokande
and DUNE at $2\sigma$ CL. Additionally, we also show the reaches by SN cooling
and by the lack of observation of high-energy neutrinos from SN1987A as black
lines. The sensitivities shown by the dashed line are obtained following [20]
with time-integrated data.

Assuming that IceCube does not observe an excess neutrino above the expected
detector backgrounds and SM neutrinos, we obtain the expected $2\sigma$
exclusion limits by finding the coupling where Eq. 4 is equal to 3.841. We
show these in Fig. 4 for the Majoron case considering both IceCube and
IceCube-Gen2. To illustrate the broad physics application of the analysis, we
also show results for the sterile neutrino case with magnetic portal in the
Supplemental Materials. Current limits [19] from the energy loss requirement,
and non-observation of high-energy SN neutrinos from SN 1987A are also shown
in Fig. 4 for the Majoron case. A future SN explosions would enable IceCube to
improve current limits by one order of magnitude. Moreover, while current
limits from time-integrated observables follow the behavior of
$g_{\phi}\propto m^{-1}_{\phi}$ as Majoron production rate is proportional to
$(g_{\phi}m_{\phi})^{2}$, IceCube can provide a stronger limit than this by
focusing on the time window, especially before the peak of SM $\bar{\nu}_{e}$
fluxes, to reduce SM condemnations.

We also present the estimated reaches of DUNE and Hyper-Kamiokande in Fig. 4
by looking for time-integrated high-energy neutrinos following Refs. [19, 20].
While the constraints from IceCube for the intermediate mass region are
comparable to that from DUNE and Hyper-Kamiokande, we observe that for both
low-mass region $m_{\phi}\lesssim 10\,\mathrm{MeV}$ and high-mass region
$m_{\phi}\gtrsim 200\,\mathrm{MeV}$, IceCube can provide stronger constraints.
This is the outcome we expect. For light $\phi$ case with negligible time
delay, IceCube would observe early signals in a time window separated from
that when most of the SM flux contribution appears, as we show in the left
panel of Fig. 3, thus enhancing the reaches. On the high mass end, neutrino
signals will arrive much later than the SM case, where a time window at late
time can be considered to reduce the standard neutrino contamination.

There are existing uncertainties from different modeling of SN explosions (see
e.g. Ref. [43]) as well as discrepancies between SN 1987A neutrino data and
modeling. Nevertheless, we point out that even the uncertainties of around
$2\sigma$ will only change our exclusion limits by a small factor. We make the
code used in this work publicly available so that the impact of SN modeling
uncertainties on the sensitivity of neutrino telescopes to BSM scenarios can
be further explored.

Summary— In this Letter, we have demonstrated that neutrino telescopes such as
IceCube can provide unprecedented probes of new physics, such as Majoron-like
bosons and neutrino magnetic moment portal, by measuring SN time profiles to a
window as short as $\mathcal{O}(0.01)$ s. Specifically, signatures outside the
SM neutrino time window, including both early and delayed signals, can be
effectively probed by neutrino telescopes. This analysis can be extended to
test other BSM realizations including, but not limited to, sterile neutrinos
with a mixing portal. Additionally, the code used to probe these models is
publicly available on GitHub.

Acknowledgements—We want to thank Thomas Janka and Daniel Kresse for providing
the data from Garching supernovae simulations—available at https://wwwmpa.mpa-
garching.mpg.de/ccsnarchive/ in machine-readable form—and also for stimulating
discussions. We also thank Jackapan Pairin for their help designing figures in
this work. Y.-Y. L is supported by the NSF of China through Grant No.
12305107, 12247103. CAA are supported by the Faculty of Arts and Sciences of
Harvard University, the National Science Foundation (NSF), the NSF AI
Institute for Artificial Intelligence and Fundamental Interactions, the
Research Corporation for Science Advancement, and the David & Lucile Packard
Foundation. CAA and JL were partially supported by the Alfred P. Sloan
Foundation for part of this work. JL is supported by the Belgian American
Educational Foundation.

## References

  * Nath  _et al._ [2010] P. Nath _et al._ , Nucl. Phys. B Proc. Suppl. 200-202, 185 (2010), arXiv:1001.2693 [hep-ph] .
  * Baumann  _et al._ [2016] D. Baumann, D. Green, J. Meyers, and B. Wallisch, JCAP 01, 007 (2016), arXiv:1508.06342 [astro-ph.CO] .
  * Hirata  _et al._ [1987] K. Hirata _et al._ (Kamiokande-II), Phys. Rev. Lett. 58, 1490 (1987).
  * Bionta  _et al._ [1987] R. M. Bionta _et al._ , Phys. Rev. Lett. 58, 1494 (1987).
  * Alexeyev  _et al._ [1988] E. Alexeyev, L. Alexeyeva, I. Krivosheina, and V. Volchenko, Physics Letters B 205, 209 (1988).
  * Loredo and Lamb [2002] T. J. Loredo and D. Q. Lamb, Phys. Rev. D 65, 063002 (2002), arXiv:astro-ph/0107260 .
  * Pagliaroli  _et al._ [2009] G. Pagliaroli, F. Vissani, M. L. Costantini, and A. Ianni, Astropart. Phys. 31, 163 (2009), arXiv:0810.0466 [astro-ph] .
  * Hüdepohl  _et al._ [2010] L. Hüdepohl, B. Müller, H. Janka, A. Marek, and G. G. Raffelt, Physical Review Letters 104, 251101 (2010), arXiv:0912.0260 [astro-ph.SR] .
  * Raffelt and Zhou [2011] G. G. Raffelt and S. Zhou, Phys. Rev. D 83, 093014 (2011), arXiv:1102.5124 [hep-ph] .
  * Argüelles  _et al._ [2019] C. A. Argüelles, V. Brdar, and J. Kopp, Phys. Rev. D 99, 043012 (2019), arXiv:1605.00654 [hep-ph] .
  * Suliga  _et al._ [2020] A. M. Suliga, I. Tamborra, and M.-R. Wu, JCAP 08, 018 (2020), arXiv:2004.11389 [astro-ph.HE] .
  * Lucente and Carenza [2021] G. Lucente and P. Carenza, Phys. Rev. D 104, 103007 (2021), arXiv:2107.12393 [hep-ph] .
  * Caputo  _et al._ [2022a] A. Caputo, G. Raffelt, and E. Vitagliano, JCAP 08, 045 (2022a), arXiv:2204.11862 [astro-ph.SR] .
  * Caputo  _et al._ [2022b] A. Caputo, G. Raffelt, and E. Vitagliano, Phys. Rev. D 105, 035022 (2022b), arXiv:2109.03244 [hep-ph] .
  * Chu  _et al._ [2019] X. Chu, J.-L. Kuo, J. Pradler, and L. Semmelrock, Phys. Rev. D 100, 083002 (2019).
  * DeRocco  _et al._ [2019] W. DeRocco, P. W. Graham, D. Kasen, G. Marques-Tavares, and S. Rajendran, JHEP 02, 171 (2019), arXiv:1901.08596 [hep-ph] .
  * Kazanas  _et al._ [2014] D. Kazanas, R. N. Mohapatra, S. Nussinov, V. L. Teplitz, and Y. Zhang, Nucl. Phys. B 890, 17 (2014), arXiv:1410.0221 [hep-ph] .
  * Magill  _et al._ [2018] G. Magill, R. Plestid, M. Pospelov, and Y.-D. Tsai, Phys. Rev. D 98, 115015 (2018), arXiv:1803.03262 [hep-ph] .
  * Fiorillo  _et al._ [2023] D. F. G. Fiorillo, G. G. Raffelt, and E. Vitagliano, Phys. Rev. Lett. 131, 021001 (2023), arXiv:2209.11773 [hep-ph] .
  * Brdar  _et al._ [2023] V. Brdar, A. de Gouvêa, Y.-Y. Li, and P. A. N. Machado, Phys. Rev. D 107, 073005 (2023), arXiv:2302.10965 [hep-ph] .
  * Acciarri  _et al._ [2015] R. Acciarri _et al._ (DUNE), (2015), arXiv:1512.06148 [physics.ins-det] .
  * Abe  _et al._ [2018] K. Abe _et al._ (Hyper-Kamiokande), (2018), arXiv:1805.04163 [physics.ins-det] .
  * Aartsen  _et al._ [2014] M. G. Aartsen _et al._ (IceCube), Phys. Rev. Lett. 113, 101101 (2014), arXiv:1405.5303 [astro-ph.HE] .
  * Aartsen  _et al._ [2018] M. G. Aartsen _et al._ (IceCube), Science 361, 147 (2018), arXiv:1807.08794 [astro-ph.HE] .
  * Abbasi  _et al._ [2023] R. Abbasi _et al._ (IceCube), Science 380, adc9818 (2023), arXiv:2307.04427 [astro-ph.HE] .
  * Abbasi  _et al._ [2022] R. Abbasi _et al._ (IceCube), Science 378, 538 (2022), arXiv:2211.09972 [astro-ph.HE] .
  * [27] L. K. (for the IceCube Collaboration), Journal of Physics: Conference Series 309, 012029 (2011).
  * Ackermann  _et al._ [2023] M. Ackermann _et al._ (IceCube), (2023), arXiv:2308.01172 [astro-ph.HE] .
  * Aartsen  _et al._ [2020] M. G. Aartsen _et al._ (IceCube), Astrophys. J. 890, 111 (2020), arXiv:1908.09997 [astro-ph.HE] .
  * Brdar  _et al._ [2021] V. Brdar, A. Greljo, J. Kopp, and T. Opferkuch, JCAP 01, 039 (2021), arXiv:2007.15563 [hep-ph] .
  * Dighe and Smirnov [2000] A. S. Dighe and A. Y. Smirnov, Phys. Rev. D 62, 033007 (2000), arXiv:hep-ph/9907423 .
  * Jaeckel  _et al._ [2018] J. Jaeckel, P. C. Malta, and J. Redondo, Phys. Rev. D 98, 055032 (2018), arXiv:1702.02964 [hep-ph] .
  * Reed [2005] B. C. Reed, Astron. J. 130, 1652 (2005), arXiv:astro-ph/0506708 .
  * Rozwadowska  _et al._ [2021] K. Rozwadowska, F. Vissani, and E. Cappellaro, New Astron. 83, 101498 (2021), arXiv:2009.03438 [astro-ph.HE] .
  * Aartsen  _et al._ [2017] M. G. Aartsen _et al._ (IceCube), JINST 12, P03012 (2017), arXiv:1612.05093 [astro-ph.IM] .
  * Griswold and BenZvi [2023] S. Griswold and S. BenZvi (IceCube), in  _38th International Cosmic Ray Conference_ (2023) arXiv:2308.01843 [astro-ph.HE] .
  * Abbasi  _et al._ [2011] R. Abbasi _et al._ (IceCube), Astron. Astrophys. 535, A109 (2011), [Erratum: Astron.Astrophys. 563, C1 (2014)], arXiv:1108.0171 [astro-ph.HE] .
  * Griswold  _et al._ [2020] S. Griswold, S. BenZvi, N. Uberoi, and R. Cross, “Icecubeopensource/asteria: v1.0.0,”  (2020).
  * Baxter  _et al._ [2022] A. L. Baxter _et al._ (SNEWS), Astrophys. J. 925, 107 (2022), arXiv:2109.08188 [astro-ph.IM] .
  * Keil  _et al._ [2003] M. T. Keil, G. G. Raffelt, and H.-T. Janka, Astrophys. J. 590, 971 (2003), arXiv:astro-ph/0208035 .
  * Aartsen  _et al._ [2021] M. G. Aartsen _et al._ (IceCube-Gen2), J. Phys. G 48, 060501 (2021), arXiv:2008.04323 [astro-ph.HE] .
  * Halzen and Raffelt [2009] F. Halzen and G. G. Raffelt, Physical Review D 80 (2009), 10.1103/physrevd.80.087301.
  * Li  _et al._ [2023] S. W. Li, J. F. Beacom, L. F. Roberts, and F. Capozzi, (2023), arXiv:2306.08024 [astro-ph.HE] .
  * Caputo  _et al._ [2022c] A. Caputo, H.-T. Janka, G. Raffelt, and E. Vitagliano, Phys. Rev. Lett. 128, 221103 (2022c).
  * Chauhan  _et al._ [2024] G. Chauhan, S. Horiuchi, P. Huber, and I. M. Shoemaker, (2024), arXiv:2402.01624 [hep-ph] .
  * Plestid [2021] R. Plestid, Phys. Rev. D 104, 075027 (2021), arXiv:2010.04193 [hep-ph] .

##  Appendix A Supplemental Methods and Tables

##  Appendix B Appendix A: Dipole Magnetic Moment Portal

SUPPL. FIG. 1: Exclusion sensitivities for the magnetic moment case. The lines
on this plot show the exclusion sensitivity for IceCube, IceCube-Gen2, and the
limits computed in [20] for Hyper-Kamiokande and DUNE. The sensitivities shown
by the dashed line have been previously computed [20] and use time-integrated
data. The shaded regions are excluded by the combination of energy loss
requirement, non-observation of photon and neutrino signals from SN1987A
studied in [20] and by constraints on the energy release from SN explosion
(SNIIP explosion) [44, 45].

Here we also present results for active-to-sterile neutrino transition
magnetic moment described by [18, 30, 20]

$\displaystyle\mathcal{L}\supset\sum_{\alpha}d_{\alpha}\bar{N}\sigma_{\mu\nu}\nu^{\alpha}F^{\mu\nu}-\frac{M_{N}}{2}\bar{N}^{c}N+\text{h.c.}\,,$
(A1)

where $\nu^{\alpha}$ and $N$ represent active and sterile neutrinos,
respectively. Further, $F^{\mu\nu}$ is the field strength tensor of the
electromagnetic field, $d_{\alpha}$ is the dimensionful coefficient of this
dimension-5 term and $M_{N}$ is sterile neutrino mass. We assume flavor
universal interaction $d_{\alpha}\equiv d$. We consider the two production
channels for $N$ inside the SN: $\nu e^{-}\to Ne^{-}$ at lower energies and
$\nu\gamma\to N$ for larger active neutrino energies. Both processes occur due
to the interaction term in Eq. A1; after $N$ are produced, they decay to
active neutrinos and photons which is again realized through the same term in
the Lagrangian, with the decay width for $N\to\nu\gamma$ given by
$\Gamma_{N}=6d^{2}M_{N}^{3}/4\pi$ [46]. The limits from IceCube are shown in
SUPPL. Fig. 1, together with previously excluded regions by energy loss, non-
observations of photons and neutrinos from SN1987A [20], and constraints on
the energy release from low-energy SN explosions (SNIIP) [44, 45]. We also
point out that by neglecting the $N$ production channel $\nu p^{+}\to Np^{+}$
inside SN, the limits we obtained are conservative relative to the SNIIP
region taken from [45]. We notice that the time delay near the exclusion
boundary is typically of $\Delta t\sim\mathcal{O}(1)$ sec, which is longer
than the standard $\bar{\nu}_{e}$ peak luminosity time scale and this is why
IceCube does not improve over the limits computed in [20] for DUNE and Hyper-
Kamiokande.


# Universal fitness dynamics through an adaptive resource utilization model

Joshua A. Dijksman  Physical Chemistry and Soft Matter, Wageningen University,
Wageningen, The Netherlands

###### Abstract

The fitness of a species determines its abundance and survival in an
ecosystem. At the same time, species take up resources for growth, so their
abundance affects the availability of resources in an ecosystem. We show here
that such species-resource coupling can be used to assign a quantitative
metric for fitness to each species. This fitness metric also allows for the
modeling of drift in species composition, and hence ecosystem evolution
through speciation and adaptation. Our results provide a foundation for an
entirely computational exploration of evolutionary ecosystem dynamics on any
length or time scale. For example, we can evolve ecosystem dynamics even by
initiating dynamics out of a single primordial ancestor and show that there
exists a well defined ecosystem-averaged fitness dynamics that is resilient
against resource shocks.

An ecosystem is a set of species that reproduce and interact by competing for
resources that support their growth. A single ecosystems can span a wide range
of length and time scales Azaele et al. (2016); Darwin already eloquently
referred to this in his “tangled bank” remark Darwin (1999). Ecosystems also
emerge on any length scale: the competition between bacterial colonies in a
Petri dish appears to have commonalities with the dynamics of competing
multinationals in the global economy. The universal nature of the emergence,
adaptation and extinction dynamics in such ecosystems suggests that ecosystem
dynamics should be amenable to simple modeling with only few ingredients that
are independent of specific physical mechanisms Nowak (2006); Azaele et al.
(2016). It seems reasonable to assume that any set of entities that have found
a thermodynamically allowed method of reproduction England (2013) is subjected
to such ecosystem dynamics. What is the simplest quantitative description that
displays all the salient features of evolution?

Here we provide a simple stochastic dynamical system from which all Darwinian
features of evolution emerge. Evolutionary dynamics modeling has settled on
describing population dynamics as a Malthusian exponential growth function
$\dot{\textbf{n}}(t)=\textbf{n}(t)f(\textbf{{n}})$ with n the set of species’
population sizes and the dot denotes a time derivative. Here $f(\textbf{n})$
is a growth rate determining function that depends on the population size and
coupling constants that specify inter-species competition, preying efficiency
Wangersky (1978); Cressman and Tao (2014). This general approach is
tremendously successful even in capturing quantitative experimental
observations of low dimensional systems Korolev et al. (2011). However,
describing the dynamics of larger ecosystems with many evolving species is
challenging, as high dimensional systems quickly lose their numerical and
analytical tractability, even without incorporating the additional complexity
of the evolution of species. At the same time, $f(\textbf{n})$ can be
considered as the _fitness_ of a species that determines its relative
abundance. We quantify fitness within the context of the ecosystem properties
by considering that species require resources to grow: in modeling ecosystem
dynamics, resource dynamics should be intricately coupled to population
dynamics.

We provide this coupling by _defining_ every species $j$ with a _strategy
vector_ $\textbf{s}_{j}$ that couples the growth dynamics
$\dot{\textbf{n}}(t)=\textbf{n}(t)f(\textbf{{n}})$ to a dynamic resource bath
vector $\textbf{r}(t)$ that represents the amount of available resources.
Here, $\textbf{s}_{j}$ describes which fraction of each resource component
$r_{i}$ is used by every species $j$ at every time step. The (time-dependent)
fitness is then naturally captured by the alignment
$\textbf{s}_{j}\cdot\textbf{r}$, and thus also time dependent growth rate of
each species. We can thus write
$\dot{\textbf{n}}_{j}=\textbf{n}_{j}f(\textbf{n},\textbf{s}_{j}\cdot\textbf{r})$
while introducing resource time dependence via a function
$\dot{\textbf{r}}=g(\textbf{n},\textbf{s})$. We will see that the perspective
provided here yields a computational framework to capture all the qualitative
phenomenology of evolution including species emergence, extinction and
adaptation.
_Strategy vector —_ Central to the modeling is the concept of a strategy
vector $\textbf{s}_{j}$ that every species $j\in\\{1\ldots k\\}$ has in order
to harvest resources. These resources are characterized by a time dependent
vector $\textbf{r}=\\{r_{1},r_{2},\ldots,r_{l}\\}$. The concept of a resource
component $r_{i}$ should be considered extremely general: it can refer to a
specific molecule, chemical energy influx or even a certain amount of habitat.
The strategy vector serves two purposes: (i) $\textbf{s}_{j}$ quantifies how
much of each resource every individual would like to take out of the resource
bath. We assume that resource consumption is proportional to $\textbf{s}_{j}$
and the population size of each species $n_{j}$. The time dynamics of every
resource component $r_{i}$ is then described by:

$\frac{\partial r_{i}}{\partial t}=-\beta\sum_{j=1\ldots
k}s_{ij}n_{j}+\gamma.$ (1)

Here $\beta$ is a timescale, and $\gamma$ is a resource replenishment factor.
(ii) $\textbf{s}_{j}$ quantifies fitness as it determines how well species $j$
is adjusted to the current resource environment. If the preferred resource
intake of the species $\textbf{s}_{j}$ is similar to the composition of the
resource environment r, the growth rate should be maximal; in case
$\textbf{s}_{j}$ and r are not aligned, the species should perform poorly. The
simplest way of incorporating species-resource dependency is then to make the
population growth factor proportional to the normalized alignment
$\frac{\textbf{s}_{j}\cdot\textbf{r}}{\|\textbf{r}\|}\equiv\Theta_{j}$ while
we always set $\|\textbf{s}_{j}\|=1$. It is now natural to write for
$n_{j}(t)$ that:

$\frac{\partial n_{j}}{\partial t}=\left(\alpha\Theta_{j}-\delta\right)n_{j}.$
(2)

Writing explicitly that $\textbf{s}_{j}\cdot\textbf{r}\equiv\sum_{i=1\ldots
l}s_{ij}r_{i}$ makes clear that $s_{ij}$ is the resource utilization
coefficient of species $j$ for resource $i$. $\alpha$ is again a time
constant; $\delta$ sets the population decay rate. The coupled set of Eqns 1 &
2 reduce to a pseudo Lotka-Volterra (LV) model Wangersky (1978) in case of one
species and one resource ($k,l=1$): predator dynamics is then described by Eq.
1 and prey dynamics by Eqn. 2. In Eqns 1 & 2, the most utilized resources also
get depleted with an increasing population size. The model thus naturally
embeds a logistic growth limitation in its dynamics, but also inter-species
competition for resources. Indeed, Eqns 1 & 2 are a simplified version of the
MacArthur equations MacArthur (1970); Chesson (1990); Haygood (2002) and will
be referred to as such from now on. We interpret the coupling matrix $s_{ij}$
however much more specifically: it is essential to see how $\textbf{s}_{j}$
here serves as the _definition_ of species $j$. This generic quantitative
definition allows us define species fitness, but also to implement species
adaptation. We can use $\textbf{s}_{j}$ to specify how a species evolves with
a simple Monte Carlo evolutionary model: we can spawn a new species $k+1$ by
taking any existing species definition $\textbf{s}_{j}$ and add a noise vector
$\eta$. The noise vector represents a small shift or mutation _amplitude_ in
the species resource utilization composition. New species can be generated
stochastically, for example by generating a new species from species $j$ when
a random number drawn from a normal distribution is larger than $\nu$ standard
deviations. $\nu$ then becomes a mutation (inverse) _rate_. Note that we
explicitly refrain from calling this “genetic drift” to avoid giving the
impression that our model can only be applied to organisms with a genetic
make-up. Indeed, by “phenotypically” defining our species solely in terms of
$\textbf{s}_{j}$, a natural link to genetic variation is lost. However, the
MacArthur equations as here interpreted and amended with the implementation of
evolution do show a remarkable ability to reproduce behavior that can be
interpreted in any evolutionary context.

Figure 1: (a) $n_{j}(t)$ for the first (red) and all subsequently emerged
species; color indicates the spawning time, the dot the emergence of a new
species. (b) $r_{i}(t)$ for the five resources available in the evolving
ecosystem. At later times, multiple resources emerge after initial depletion.

To see this, we implement the MacArthur equations in MATLAB with a fourth
order Runge-Kutta (RK) scheme; see Supplementary Material SM  for details. We
can distinguish many different cases, depending on the (initial) number of
species $k$ and total number of resources $l$. The phenomenological richness
of the coupling embedded in the our stochastic MacArthur equations can however
already be observed for starting evolution with the most stringent starting
condition of _one_ species $k=1$; we choose an $l=5$ resource space. Note that
capturing the emergence of new interacting species from a primordial
reproducing entity is an explicit aim of the current modeling approach:
evolving an ecosystem from a single species is for example not possible in the
classical MacArthur, LV or replicator equation context.
We show typical ecosystem evolution initiated from one species $k=1$ starting
at size $n_{1}=1$ in Fig. 1; all further settings can be found in SM SM . Due
to the stochastic nature of species emergence and extinction, every
realization of ecosystem dynamics is different. However, several important
qualitative features reproduce and are visible in any example: (i) The initial
species size oscillates in time until a viable new species has emerged; in
sync, the resource dynamics is also oscillatory for the the smallest component
of $\textbf{s}_{j}\equiv\min(\textbf{s}_{1})$ as here the resource usage stays
dynamically balanced with the resource influx. (ii) the emergence of new
species affects the timescale of periodic oscillations; also new species can
make older species go extinct. (iii) later in the evolution, the population
fluctuations shift in frequency and decay in amplitude and multiple resources
become utilized. The interpretation of these trends is obvious: the randomly
selected initial species favors the survival of one resource, for which
$s_{i}$ is smallest. After this transient, the dynamics follow the $l=1$,$k=1$
system which is pseudo-LV in character and allows for periodic orbits of fixed
frequency. In this phase, the possibility of the random emergence of new
species is consequential; the emergence of new species that are $\eta$
different from their parents will suppress the dominant role of the first
species and thus limit its overuse of other resources, which makes the
remaining resources emerge again as they are always continuously replenished
at rate $\gamma$.

Figure 2: The phylogenetic tree derived from the emergence of species in the
evolution shown in Fig. 1. The vertical scale is arbitrary; the horizontal
scale depicts time. The species number is indicated in color consistent with
labeling in Fig. 1 and on the right to provide the chronology of emerged
species in the tree. The red dashed line shows the lineage of the species with
which the tree started; it went exinct around $t\sim 2.6\times 10^{4}$ steps.

_Speciation & Emerging fitness —_ Speciation is naturally embedded in the
current implementation of the stochastic MacArthur equations. Starting from
$k=1$, every species can spawn a new species at any time step. This generates
a “tree” of evolution; again the example ecosystem from Fig. 1a can be used as
a reference point. In this evolving ecosystem, many new species emerged and
went extinct. It is however not transparent from Fig. 1a,c from which existing
species the new species emerged. To visualize this, we present the
phylogenetic tree of this ecosystem in Fig. 2. The binary tree shows several
long and short branches representing the successful and unsuccessful species
in the ecosystem.
Even though the time dependent dynamics of $\textbf{r}(t)$ and the ecosystem
make-up $\textbf{n}(t)$ is intrinsically stochastic and can be very complex
depending on $\nu$, $\eta$ and $l$, we see immediately that the model
naturally embeds a fitness that grows with time. The alignment between
resource utilization and availability is a natural representative of fitness:
$\Theta_{j}$ increases with time for successive generations $j$. We can show
this already in the example ecosystem introduced in Fig. 1a,b. Only the new
species for which $\alpha\textbf{s}_{1}\cdot\textbf{r}>\delta$ will survive
the initial dynamics, but the example above shows that many newly spawned
species satisfy this condition. We take the example ecosystem from Fig. 1a and
compute the alignment $\Theta_{j}$ first for all species that were active in
the ecosystem – see Fig. 3a. We scale out the effect of arbitrary starting
conditions by showing $\Theta_{j}/\min(\textbf{s}_{j})-1$. Clearly the fitness
of the individual species is increasing with time, even though large
fluctuations are visible. It seems natural to see fitness growth originating
in the fittest species having the largest growth rates, yet this picture is
naive: a large population size also removes more resources, leading to an
overall worse alignment $\textbf{s}_{j}\cdot\textbf{r}$.

Figure 3: (a) $\Theta_{j}/\min(\textbf{s}_{1})-1$ for all species that emerged
during the evolution; color coding is the same as in (a). The filled circle
shows species emergence; extinction of an individual species is denoted as a
vertical line. (b) $\langle\overline{\Theta}/\Theta_{1}(t=0)\rangle$ for
$\eta=0.005,~{}l=5$ (red, dashed), $\eta=0.05,~{}l=5$ (blue, dashed) and
$\eta=0.01,~{}l=5,15,50,150$ (green, solid). Short arrows indicate the three
main features present in all ecosystem dynamics: i resource depletion, ii r
restoration due to new species emergence, iii overall fitness growth in larger
ecosystem balancing. Long arrows indicate the trend with $\eta$ and $\l$.

We also demonstrate that the trend of growing fitness is a natural consequence
of the stochastic MacArthur equations. We analyze aggregate dynamics of many
ecosystem repetitions in which we only vary the initial $\textbf{s}_{j}$ and,
due to the random emergence of new species that evolve with noise vector
$\eta$, the composition of new species. We define the mean ecosystem fitness
$\overline{\Theta}=\sum n_{j}\Theta_{j}/N(t)$ with $N=\sum n_{j}$, weighing
the contribution of each living species with its contemporary population size.
This mean ecosystem fitness still has an arbitrary initial amplitude
$\overline{\Theta}(t=0)$ set by the initial species composition, which we
again scale out while taking the average over 100 realizations of each
ecosystem. We run all ecosystems starting from one species, with
$\alpha,=0.005$, $\beta=0.01$, $\delta=0.1$, $\nu=3.9$; for robustness we now
test $\gamma_{i}=1-(i-1)/l$, $\eta=0.005,0.05$ and $l=5,15,50,150$. The
resulting $\langle\overline{\Theta}/\Theta_{1}(t=0)\rangle$ are shown in Fig.
3b. Remarkably, the aggregate dynamics shows a characteristic time evolution
for _all_ combinations of parameters. Several features stand out: (i) Initial
fitness rapidly decays due to total depletion of all but one resource;
$\overline{\Theta}$ reaches a first minimum where total ecosystem
disappearance ($n_{j}<0.1\forall j$) occurs in about 10% of all runs due to
the vanishing of resources. (ii) After that first minimum a maximum is reached
due to the emergence of new species who rapidly restructure the composition of
r. (iii) The second stage of evolution displays another broader minimum:
initially the emergence of more new species put pressure on r yet eventually
the enhanced diversity in s allows for the re-emergence of all resources (as
in Fig. 1b) and an overall ability to yet again improve fitness. By varying
the model ingredients, we can interpret the role of noise and dimensionality:
the ecosystem fitness dynamics is not affected by the level of noise $\eta$
apart from a shift in the inflection point where fitness growth eventually
becomes significant. This suggests that _any_ magnitude in species composition
drift overall leads eventually to the emergence of a growing fitness in the
ecosystem, hence that discrete evolutionary steps are not required for
evolution. The caveat is that with a small evolutionary drift, an ecosystem
might be wiped out before it reaches the stage of fitness growth. The
dimensionality of the resource space $l$ affects average ecosystem dynamics
only by modifying time scales and fitness amplitudes: for smaller $l$, overall
fitness drops to lower values and regains overall fitness at a lower pace,
suggesting that the expansion of species’ abilities or habitats provides a
positive feedback mechanism for ecosystem expansion, naturally accelerating
evolution.

_Adaptation —_ The resource alignment interpretation of the stochastic
MacArthur equations make their many generalizations physically meaningful and
provide clear links to real world ecosystems. Empirically the most promising
variation is to provide a time varying resource influx by introducing
$\gamma_{i}(t)$, sometimes also called a “pulse” experiment Tilman (1987). We
implement an example where we change $\gamma_{i}=1-(i-1)/l$ for $t<t_{c}$ and
$\gamma_{i}=i/l$ for $t\geq t_{c}$, which keeps the $\|\gamma\|$ constant at
all times. We choose $t_{c}=2\times 10^{4},3\times 10^{4},4\times 10^{4}$ This
abrupt change in resource influx induces adaptation by the ecosystem. Solving
the MacArthur equations with this time dependent resource influx over 100
realizations at previously defined $\alpha,\beta,\delta,\nu$ and with
$\eta=0.01,l=5$, we observe that ecosystems are able to recover from such a
resource shock; example time dynamics is shown in Fig. 4a. The fitness goes
_up_ after the resource shock, as the Malthusian pressure on resource
availability and resulting fitness decrease is temporarily halted. After a
resource shock, ecosystem dynamics converges back to the general fitness
dynamics as determined by the system parameters. The mean number of _living_
species $\langle k\rangle$ in the ecosystems exposed to a resource shock is
not affected – see Fig. 4b. The mean number of living species keeps growing
with time without showing any signatures at the resource shock moments. The
dynamics of available resources and population size are however affected,
providing clear empirical benchmarks.

Figure 4: (a) Average mean ecosystem fitness as a function of time with a
resource shock $\gamma_{i}=1-(i-1)/l\rightarrow\gamma_{i}=i/l$ applied at
times indicated by the long arrows. (b) The mean number of living species in
the ecosystems from panel (a). Clearly the resource shock does not affect the
growth dynamics of the number of species.

The stochastic MacArthur equations allow for many further generalizations, as
partly noted before MacArthur (1970); Chesson (1990); Haygood (2002). It is
obvious that $\alpha$, $\beta$, $\gamma$ and $\delta$ can be made time
dependent and non-uniform among species. Predator dynamics can be introduced
by adding another predator coupling matrix term $\sum_{i}m_{ij}n_{i}$, which
can have both positive and negative elements, when species $j$ is a predator
or prey respectively. The addition of predator matrix coupling then allows for
some species to take up natural resources while others in the ecosystem only
prey on existing species, to model for example marine ecosystems with plankton
at the base. In Eq. 1 we have only considered that species _consume_
resources; they may also _provide_ resources. The required sign change of
components of $\textbf{s}_{j}$ would allow for the emergence of mutualism.
Oxygenic photosynthesis Knoll and Nowak (2017) is one example; on a different
scale also gut microbes provide natural resources for each other Faust and
Raes (2012). The criteria for the generation of new species can be
significantly generalized, for example by creating more new species when
$\textbf{n}_{j}$ becomes larger, or by making the noise amplitude dependent on
the species size: $\eta(\textbf{n}_{j})$: in an evolutionary bottleneck or in
large populations, the composition drift may be faster or slower. Note that
replicator or LV-models can also be generalized to spawn new species in the
way described here, and the species composition in the LV model is in some
sense “defined” by the coupling matrix $m_{ij}$, to which stochastic
components can be added. The current model however offers the essential
benefit that it can initiate evolution from a single reproducing entity
coupled to a resource bath, making it more physically realistic and more
general than the LV model.

_Conclusions —_ We presented a quantitative model for evolution that
reproduces the phenomenology of species emergence, speciation, extinction and
adaptation. The model hinges on quantifying fitness by defining a species
through its resource uptake — the species _is_ the resource strategy. The
quantitative nature of the model provides us with a new path to explore basic
questions about evolution. For example: evolution must occur in small steps,
yet species emerge often in radically different phenotype. Does the current
model retain the same dynamics in the limit of $\eta\rightarrow 0$? If so, the
evolutionary noise may simply set an additional time scale for evolutionary
dynamics. How does the phylogenetic tree structure depend on the noise
amplitude and species creation threshold? The dimensionality of the resource
space is also a crucial element that requires deeper study: in the limit of
large $l$, the curious geometry of high dimensional vector spaces dictates
that the distance between two species as measured by their angle will always
be finite, suggesting that the emergence of new species becomes a natural
outcome of the quantitative structure of evolution.
_Acknowledgments_ We thank Christian Fleck for several useful comments and for
mentioning that $\textbf{s}_{ij}$ can have either sign. Discussions with
Kirsty Wan on ordinary differential equation solvers are gratefully
acknowledged.

## References

  * Azaele et al. (2016) S. Azaele, S. Suweis, J. Grilli, I. Volkov, J. R. Banavar, and A. Maritan, Reviews of Modern Physics 88 (2016).
  * Darwin (1999) C. Darwin, _On The Origin of Species, 6th Ed._ (Project Gutenberg, 1999).
  * England (2013) J. L. England, The Journal of Chemical Physics 139, 121923 (2013).
  * Nowak (2006) M. A. Nowak, _Evolutionary Dynamics_ (Harvard University Press, 2006).
  * Wangersky (1978) P. J. Wangersky, Annual Review of Ecology and Systematics 9, 189 (1978).
  * Cressman and Tao (2014) R. Cressman and Y. Tao, Proc Natl Acad Sci U S A 111 Suppl 3, 10810 (2014).
  * Korolev et al. (2011) K. S. Korolev, J. B. Xavier, D. R. Nelson, and K. R. Foster, Am Nat 178, 538 (2011).
  * MacArthur (1970) R. MacArthur, Theoretical Population Biology 1, 1 (1970).
  * Chesson (1990) P. Chesson, Theoretical Population Biology 37, 26 (1990).
  * Haygood (2002) R. Haygood, Theor Popul Biol 61, 215 (2002).
  * Tilman (1987) D. Tilman, The American Naturalist 129, 769 (1987).
  * Knoll and Nowak (2017) A. H. Knoll and M. A. Nowak, Science Advances 3 (2017).
  * Faust and Raes (2012) K. Faust and J. Raes, Nat Rev Microbiol 10, 538 (2012).
  * (14) See Supplemental Material at [URL will be inserted by publisher] for details about the numerical implementation of stochastic MacArthur equation.