# Memory-triggered supernova neutrino detection

Mainak Mukhopadhyay<EMAIL_ADDRESS>Department of Physics, Arizona State
University PO Box 871504, Tempe, AZ 85287-1504, USA.  Zidu Lin<EMAIL_ADDRESS>Department of Physics and Astronomy, University of Tennessee Knoxville, TN
37996-1200, USA  Cecilia Lunardini<EMAIL_ADDRESS>Department of
Physics, Arizona State University PO Box 871504, Tempe, AZ 85287-1504, USA

###### Abstract

We demonstrate that observations of the gravitational memory from core
collapse supernovae at future Deci-Hz interferometers enable time-triggered
searches of supernova neutrinos at Mt-scale detectors. Achieving a sensitivity
to characteristic strains of at least $\sim 10^{-25}$ at $f\simeq 0.3$ Hz –
e.g., by improving the noise of DECIGO by one order of magnitude – will allow
robust time triggers for supernovae at distances $D\sim 40-300$ Mpc, resulting
in a nearly background-free sample of $\sim 3-70$ neutrino events per Mt per
decade of operation. This sample would bridge the sensitivity gap between rare
galactic supernova bursts and the cosmological diffuse supernova neutrino
background, allowing detailed studies of the neutrino emission of supernovae
in the local Universe.

##  I Introduction

Neutrinos are major players in the emerging field of multi-messenger
astronomy. With gravitational waves (GWs) and photons, they have the potential
to probe the most extreme astrophysical phenomena in unprecedented detail.
Core collapse supernovae (CCSNe) are prime targets of multi-messenger
observations, where neutrinos dominate the energy output and carry direct
information on the extremely dense environment surrounding the collapsed core.
The $\sim 10$ s burst of neutrinos from a supernova will also allow tests of
particle physics beyond the Standard Model [1, 2, 3, 4].

Figure 1: The number of memory-triggered supernova neutrinos detected at a 1
Mt water Cherenkov detector in 30 years, as a function of the noise of the GW
detector at $f=0.3$ Hz. The vertical lines mark specific experiments
considered here. The lower and upper shaded regions refer respectively to a
homogeneous population with moderate memory strain and a mixed population
where 40% of collapses have stronger memory strain; the shading describes the
effect of varying the neutrino spectrum, see Table II.2. The dots (upper set:
NSFC and lower set: BHFC) and legends on the curves give the GW distance of
sensitivity ($r^{GW}_{max}$, see text below Eq. (2)) corresponding to the
noise on their abscissa.  Figure 2: Solid: the characteristic gravitational
memory strain $h_{c}(f)$ for the NSFC and BHFC models (thin and thick lines
respectively). The distance to the supernova is $r=1$ Mpc. Dashed: sky-
averaged noise curves for representative detectors (see fig. 1).

The detection of an individual supernova neutrino burst is exciting as well as
challenging. A statistically significant observation is possible only for
supernovae within 1-3 Mpc from Earth [5, 6], where collapses are rare,
resulting in decades of waiting time. An alternative is to search for the
Diffuse Supernova Neutrino Background (DSNB), from all the supernovae in the
universe [7, 8, 9, 10], which has a substantial cosmological component.
${\mathcal{O}}(10-100)$ DSNB neutrinos could be detected in a decade (see,
e.g., [11]), and preliminary data could be available in just a few years [12,
13, 14, 15, 16, 17].

Burst and DSNB searches lack sensitivity to the local universe, $r\sim 3-100$
Mpc, where many supernova-rich galaxies are situated. Ideas to overcome this
gap typically rely on time-triggers that would allow to identify a single
neutrino as signal instead of background. One could use either a neutrino
self-trigger — where $2-3$ neutrinos observed less than 10 s apart can be
attributed to a supernova with high confidence [5, 18] — , or the time
coincidence with the $O(10^{2})$ Hz supernova GW signal from interferometers
like LIGO-Virgo and its successors [19, 20, 21]. Both methods are still
limited to a few Mpc distance, except for the most optimistic GW models (see,
e.g., [22] and references therein) and futuristic multi-Megaton neutrino
detectors [6]111Astronomical observations of supernovae can not serve as time
triggers, due to the ${\mathcal{O}}(1)$ hour uncertainty in the time delay
between the neutrino and the electromagnetic signal from the same star..

Figure 3: Detection probabilities for a memory signal, $P^{GW}_{det}(r)$, at
three GW detectors of reference, and neutrino detection probabilities,
$P^{\nu}(1,r)$ and $P^{\nu}(2,r)$ (see eq. (4)). Shadings describe the
variations with the varying neutrino spectrum, see Table II.2. The left
(right) panel is for NSFC (BHFC). Also shown is the cumulative rate of core
collapses (vertical axis on the right). See labels on the curves for details.

.

In this paper, we propose a new time-triggered method to detect supernova
neutrinos, which is potentially sensitive to supernovae up to $\sim$ 100 Mpc.
The time trigger is the observation of the gravitational memory signal caused
by the neutrino emission itself. The memory is a non-oscillatory, permanent
distortion of the local space time due to the anisotropic emission of matter
or energy by a distant source. The memory due to neutrino emission by a
supernova at distance $r$ has characteristic strain $h_{c}\sim
10^{-23}-10^{-21}(\mathrm{10~{}kpc}/r)$ and frequencies in the Deci-Hz band,
$f\sim 0.1-3$ Hz [23, 24, 25, 26, 27, 28, 29]. The memory develops $\sim 0.1$
s from the start of the neutrino emission, thus being an ideal time-trigger.
Next generation powerful Deci-Hz GW detectors, like the Deci-hertz
Interferometer Gravitational wave Observatory (DECIGO) [30, 31, 32, 33] and
the Big Bang Observer (BBO) [31] will provide robust triggers for supernovae
at 10 Mpc and beyond [34]. These would result in a nearly pure sample of $\sim
10-100$ supernova neutrino events from the local universe within a few
decades; see our summary figure, fig.1. Here we illustrate our proposed
methodology and its physics potential.

##  II Formalism

###  II.1 Gravitational memory signals

The supernova neutrino memory strain can be expressed as [35, 36, 24]

$h^{xx}_{TT}=h(r,t)=\frac{2G}{rc^{4}}\int_{-\infty}^{t-r/c}dt^{\prime}L_{\nu}(t^{\prime})\alpha(t^{\prime})\,.$
(1)

where $c$ is the speed of light, $t$ is the time post bounce and $G$ is the
universal gravitational constant. $L_{\nu}$ is the all-flavors neutrino
luminosity and $\alpha\sim{\mathcal{O}}(10^{-3}-10^{-2})$ is the time-varying
anisotropy parameter [25, 28]222 In axisymmetric simulations [24, 25], only
the ${}^{\prime}+^{\prime}$ strain may be extracted, and the observer is
positioned such that the ${}^{\prime}+^{\prime}$ strain is maximized..
Simulations show that $\alpha(t)$ becomes non-zero within a few ms post-
collapse, during the accretion phase, and can change sign multiple times
within the first second, as a result of the dynamics of the matter near the
collapsed core. The behavior of $\alpha(t)$ at $t>1$ s, during the cooling
phase, is unknown. Following [34], we consider two phenomenological models for
the memory: the first, characterized by a weaker and shorter anisotropic
phase, is representative of a neutron-star-forming collapse (NSFC); the second
has larger and prolonged anisotropy, and could represent a black-hole-forming
collapse (BHFC). In both models, $\alpha=0$ for $t>1$ s. Maximum values of
$rh(r,t)\sim 26.5$ cm and $rh(r,t)\sim 400$ cm are obtained for the two models
respectively. In fig. 2, we show the memory characteristic strain [27],
$h_{c}(r,f)=2f|\tilde{h}(r,f)|$, where $\tilde{h}(r,f)$ is the Fourier
Transform of $h(r,t)$. Also shown are the noise curves of Deci-Hz detectors,
which are given by the quantity $h_{n}(f)=\Upsilon\sqrt{5\ f\ S_{n}(f)}$ [27],
where $S_{n}(f)$ is the power spectral noise density [37]. We choose
$\Upsilon=1,10^{-1},10^{-3}$; the first and last correspond to DECIGO and its
optimal (futuristic) realization, Ultimate DECIGO [30, 31, 32]; the middle
value represents an hypothetical intermediate case (DECIGO+ from here on).

The detectability of a memory signal is determined by the signal-to-noise
(SNR) ratio of the detector333Here the comparison with published SNR curves
has indicative character only; a signal-specific study of the detectability is
ultimately needed, and is left for future work. , which is defined as [38]

$\rho^{2}(r)=\int_{-\infty}^{\infty}d(\text{log}f)\Bigg{(}\frac{h_{c}(r,f)}{h_{n}(f)}\Bigg{)}^{2}\,.$
(2)

We compute the probability of detecting a CCSN memory, $P_{det}^{GW}$, for a
fixed false alarm probability $P_{FA}^{GW}=0.1$. This requires producing
Receiver Operating Curves (ROCs) in the plane $P_{det}^{GW}-P_{FA}^{GW}$,
which we do following the formalism in [39] for $N=3$ degrees of freedom (here
$N$ is set equal to the number of Gaussian functions used to represent
$\alpha(t)$, see [34])444In Ref. [39], the formalism of $P_{det}$ and $P_{FA}$
are presented in the context of matched filter analysis. In the search for
gravitational memory signals, we applied a filter studied in [34], which
reasonably reproduce the results from numerical simulations.. The result is
that $P_{det}^{GW}$, at a fixed $P_{FA}^{GW}$, is an increasing function of
$\rho(r)$, through which it depends on the distance, $r$. We define the GW
detector distance of sensitivity, $r^{GW}_{max}$ such that
$P^{GW}_{det}(r^{GW}_{max})=0.5$. $P_{det}^{GW}(r)$ is shown in fig. 3 for our
cases of reference. For DECIGO, and for NSFC and BHFC respectively, we have
$r^{GW}_{max}\simeq 4$ Mpc and $r^{GW}_{max}\simeq 33$ Mpc. We find
$r^{GW}_{max}\simeq 40$ Mpc and $r^{GW}_{max}\simeq 335$ Mpc for DECIGO+; for
Ultimate DECIGO, $r^{GW}_{max}>350$ Mpc for both population models.

We note in passing that, in principle, the stochastic effect of the memory
signals from cosmological supernovae contributes to the noise in a GW
detector, and therefore to $r^{GW}_{max}$. For real-time searches of transient
signals at a modern interferometer like LIGO, the noise spectral density is
measured over sliding time windows of ${\mathcal{O}}(10^{2})$ s width, leading
to a fast identification of seconds-long transients [40]. Due to the high
supernova rate ($\dot{\rho}_{SN}\sim 10^{-4}~{}\mathrm{yr^{-1}~{}Mpc^{-3}}$
locally, corresponding to $\sim 10^{7}$ core collapses per year in the visible
universe) [41, 42, 43], the individual cosmological memory signals combine to
constitute a continuum, that is best described by an integral over the cosmic
volume. Such integral represents the contribution of supernovae to the
fraction of cosmic energy density in GW, $\Omega_{GW}$ (see, e.g., [44, 45,
46, 47] for the formalism).

Following Ref. [45], we have estimated the supernova memory contribution to
$\Omega_{GW}$, and found that it affects the probability curves in fig. 3
solely for Ultimate DECIGO, and only for $r\gtrsim 300$ Mpc and for the most
optimistic memory model (BHFC curve in fig. 2, corresponding to a GW spectral
energy density $\Omega_{\mathrm{G}W}=\mathcal{O}(10^{-17})$). As will be seen
in the next section, the triggered neutrino search is limited to $r\lesssim
300$ Mpc by the background at the neutrino detector. Therefore, the stochastic
GW noise from supernovae is negligible in the present context, and will not be
considered further.

###  II.2 Neutrino signals

For neutrino detection, we consider a water Cherenkov experiment, where the
main channel of sensitivity is inverse beta decay (IBD),
$\mathrel{{\bar{\nu}}_{e}}+p\rightarrow n+e^{+}$. For the time-integrated
(over $\Delta t=10$ s) $\mathrel{{\bar{\nu}}_{e}}$ flux at Earth,
$\Phi(r,E_{\nu})$ we use analytical quasi-thermal spectra of the form given in
[48]. The average $\mathrel{{\bar{\nu}}_{e}}$ energy is varied in an interval
motivated by numerical simulations [49, 50, 51], in a way to effectively
account for neutrino oscillations. The spectrum shape parameter, $\beta$, and
the total energy in $\mathrel{{\bar{\nu}}_{e}}$ are fixed. See Table II.2 for
details.

Model |  Energy |  |  |
---|---|---|---|---
$(\times 10^{53}\text{ergs})$ |  $\beta$ |  $\langle E_{\nu}\rangle$ |  |
(in MeV) |  |  |  |  |
|  Ac. ph. |  $\bar{\nu}_{e}$ |  |  Lower |  Upper
NSFC |  $1.2$ |  $0.5$ |  $3$ |  $11$ |  $15$
BHFC |  $2$ |  $0.45$ |  $2$ |  $15$ |  $20$
Table 1: The neutrino flux parameters, from numerical simulations [49, 50,
51]. The Ac. ph. and $\mathrel{{\bar{\nu}}_{e}}$ columns refer to the all
flavor energy in the accretion phase only (which contributes to the memory
signal, see text) and to the energy in $\mathrel{{\bar{\nu}}_{e}}$ emitted
over the time-triggered interval of 10 s. $\beta$ is related to the second
momentum of the spectrum: $\beta=(2\langle E_{\nu}\rangle^{2}-\langle
E_{\nu}^{2}\rangle)/(\langle E_{\nu}^{2}\rangle-\langle E_{\nu}\rangle^{2})$.

The predicted number of events in the detector from a CCSN at distance $r$ is:

$N(r)=\int_{E^{th}_{\nu}}^{E^{max}_{\nu}}N_{p}\eta\sigma(E_{\nu})\Phi(r,E_{\nu})\
dE_{\nu}\,,$ (3)

where $N_{p}$ is the number of target protons, $\eta=0.9$ is the detection
efficiency [52, 53, 54] and $\sigma(E_{\nu})$ is the IBD cross-section [55].
We take an energy interval
$\left[E^{th}_{\nu},E^{max}_{\nu}\right]=[19.3,50]~{}\mathrm{MeV}$ to avoid
the spallation background at low energy and the atmospheric neutrino
background at high energy [53, 15, 56]. We find $N(1~{}\mathrm{Mpc})\simeq
5-12$ and $N(1~{}\mathrm{Mpc})\simeq 12-18$ for NSFC and BHFC respectively, by
varying the mean $\mathrel{{\bar{\nu}}_{e}}$ energy in the intervals given in
Table II.2.

The Poisson probability of observing $N\geq N_{min}$ neutrino events in a
detector is

$P^{\nu}(N_{min},r)=\sum_{n=N_{min}}^{\infty}\frac{N^{n}(r)}{n!}e^{-N(r)}\,.$
(4)

It is plotted for $N_{min}=1,~{}2$ in fig. 3 for the two models of reference.
As expected, $P^{\nu}(N_{min},r)$ declines rapidly at $r\gtrsim 3$ Mpc.

Figure 4:  Number of background events and of memory-triggered neutrino events
from collapses at distance $r<D$, as a function of $D$, for a Mt water
Cherenkov detector and 30 years running time. The upper to lower shaded
regions are for triggers from Ultimate DECIGO, DECIGO+ and DECIGO (the latter
is invisible in the left panel). Shadings describe the effect of varying the
neutrino spectrum, see Table II.2. Left panel: homogeneous NSFC population.
Right: mix of 60% NSFC and 40% BHFC.

##  III Memory-triggered neutrino observations

###  III.1 Event rates

To estimate the rate of memory-triggered neutrino events, we model the rate of
core collapses as a function of $r$. For $r\lesssim 11$ Mpc, we use the rates
for individual galaxies from [57]. For $r>11$ Mpc we assume a constant
volumetric rate of $R_{SN}=1.5~{}10^{-4}~{}\mathrm{Mpc^{-3}yr^{-1}}$ (the
evolution with redshift is negligible for the distances of interest here). The
cumulative rate (total rate of core collapses with $r<D$) is shown in fig. 3.

The number of memory-triggered neutrino events from all supernovae within a
distance $D$, over a detector running time $\Delta T$ can be calculated as a
sum over all the galaxies (index $j=1,2,...$), at distance $r_{j}<D$:

$N^{trig}_{\nu}(D)=\Delta
T\sum_{j,r_{j}<D}R_{j}N(r_{j})P_{det}^{GW}(r_{j})\,,$ (5)

where $R_{j}$ indicates the supernova rate in the galaxy $j$. This discrete
expression is replaced by a continuum one, involving an integral, for $D>11$
Mpc, where the cosmological supernova rate is used.

We now discuss the background of the time-triggered neutrino search. The
number of supernova memory signals observed in the time $\Delta T$ is,
$N^{trig}_{SN}(D)=\Delta T\sum_{j,r_{j}<D}R_{j}P_{det}^{GW}(r_{j})$, and the
number of expected background events is
$N^{trig}_{bckg}(D)=N^{trig}_{SN}(D)\lambda\Delta t$, where $\lambda\simeq
1313$ events/year is the background rate in the detector [53, 15, 56]. Note
that the background level is reduced by a factor
$\epsilon_{bckg=}N^{trig}_{SN}(D)\Delta t/\Delta T$ compared to an un-
triggered search555The time delay effect due to the non-zero neutrino mass can
be neglected; it is estimated to be only a fraction of a second for energies
and distances of interest here, see, e.g., [58]. .

We limit our study to neutrino events (eq. (5)) from CCSNe in the cosmic
volume with $4<D<350$ Mpc, thus accounting for the fact that a nearby
supernova ($D<4$ Mpc) is unlikely to occur in three decades time. The upper
bound on $D$ is justified because beyond it the total event rate becomes
dominated by background. Experimentally, a distance cut can be accomplished in
different ways. For NSFC, one can make a selection using estimates of $D$ from
astronomy follow ups, which will benefit from the alerts from the memory
detection and should have excellent sensitivity to supernovae in the local
universe (see, e.g., [59, 60, 61, 62, 63] for dedicated projects). In the
absence of an optical counterpart (BHFC), a similar (although less efficient)
data selection can be performed using minimal input from theoretical models,
e.g. to obtain conservative upper limits on the distances of individual
observed BHFCs via GW memory signals. In the mature stage of this search,
specifically designed data-analysis algorithms – exploiting the correlation of
multiple observables – could reduce the level of model-dependency to a
minimum.

###  III.2 Results

Our main results are in fig. 4 and fig. 1 for $\Delta T=30$ yrs and for two
scenarios: (i) a supernova population entirely comprised of NSFC; and (ii) a
mixed population with 60% NSFC and 40% BHFC. Fig. 4 shows $N^{trig}_{\nu}(D)$
as a function of $D$. We observe the (expected) trend
$N^{trig}_{\nu}(D)\propto D$ for $D\lesssim r^{GW}_{max}$ 666Recall that, in
the continuum limit, the number of supernovae scales like $D^{3}$ and the flux
dilution factor like $D^{-2}$., with a flattening of the curves at larger $D$
due to the loss of sensitivity of the GW detector. For case (i), time triggers
from DECIGO+ will result in $N^{trig}_{\nu}\sim 10-30$. For Ult. DECIGO,
$N^{trig}_{\nu}\sim 100-300$ is expected777For comparison, our estimated
number of CCSNe within 350 Mpc is $N^{trig}_{SN}\sim 1.21\times 10^{6}$.. For
the mixed population (case (ii)), results for Ult. DECIGO change only
minimally, due to the different neutrino parameters between NSFC and BHFC.
Instead, $N^{trig}_{\nu}$ increases dramatically, surpassing 100, for DECIGO+,
due to the larger distance of sensitivity to BHFC. Indeed, the number of
triggered neutrino events from collapses with $30<D<350$ Mpc is dominated by
BHFC (see also fig. 3). For this mixed population scenario, even DECIGO could
be effective, providing a few triggers of BHFC up to $D\sim 30$ Mpc, resulting
in $N^{trig}_{\nu}\sim 10$. As fig. 4 shows, in all cases the signal exceeds
the background for triggers with $r\lesssim 100$ Mpc. For Ultimate DECIGO,
even for the largest $D$ the signal is comparable to the background, and would
cause a statistically significant excess.

Our summary figure, fig. 1, shows $N^{trig}_{\nu}(350~{}{\mathrm{M}pc})$, as a
function of $h_{n}$, together with representative values of $r^{GW}_{max}$.
Roughly, we find $N^{trig}_{\nu}\propto 1/h_{n}$, for $h_{n}\gtrsim 10^{-26}$,
with a flattening at lower values of $h_{n}$, due to upper cutoff on $D$. It
appears that, even for the most conservative parameters, a ${\mathcal{O}}(10)$
noise abatement with respect to DECIGO (i.e., DECIGO+) is sufficient to obtain
a signal at a Mt scale detector in $\sim 20-30$ years.

##  IV Conclusions and Discussion

Summarizing, we have described a new multimessenger approach to core collapse
supernovae, where a time-triggered search of supernova neutrinos is enabled by
observing the gravitational memory caused by the neutrinos themselves. This
scenario could be realized a few decades from now, when powerful Deci-Hz
interferometers (noise $h_{n}\lesssim 10^{-25}$) and Mt-scale neutrino
detectors start operating. For optimistic parameters, DECIGO and
HyperKamiokande (mass $M=0.260$ Mt) might already achieve a low statistics
observation. This approach will also enable joint analyses of neutrino, GW and
light curves of CCSNe in local universe.

Our proposed method will deliver a sample of neutrino events from supernovae
in the _local_ universe, from which the main neutrino properties – i.e, the
(population-averaged) energy spectra and time profiles– will be measured.
These can then be compared to the same quantities from (1) SN 1987A, to
measure the deviation between SN1987A and an average local supernova (the same
exercise can be done for a future nearby supernova burst, if it occurs); (2)
the DSNB, to distinguish the contributions to the DSNB by CCSNe in the distant
universe and by other transients (e.g. binary mergers). The comparison between
cosmological and local contributions to the DSNB will test hypotheses of how
the supernova progenitor population evolves with the distance. Even within the
local-neutrino sample, one could test the evolution with distance, if the
latter is estimated for each supernova using multi-messenger observations
(e.g., the amplitude of the memory signal and astronomical imaging).

Correlating memory and neutrino data might reveal two distinct populations,
like those described here (NSFC and BHFC), which could be statistically
separated. For example, events having a relatively large neutrino-memory time
separation (bigger than 1 s, as black hole formation typically occur within 1
s, cutting off the neutrino luminosity [64, 65, 66]) and(or) followed by
electromagnetic (EM) signals of a CCSN could be attributed to NSFC. The
possibility to study such sub-population individually is unique of this local-
collapses neutrino sample. Additionally, our method provides a unique chance
to jointly analyze neutrino and follow-up EM signals [67, 59] from the _same_
NSFC. Although only $\approx 1$ event would be detected from a specific NSFC,
it can help to determine the time when the core of a NSFC collapses and the
shock is formed. Such estimation would be relatively precise, considering that
the neutrino burst from a NSFC only lasts for $\approx 10$ s. A supernova EM
signal is delayed relative to the neutrinos, by at least the time it takes the
shock to propagate through the envelope, typically hours. Measuring this time
delay will provide a crucial confirmation and can test the variation of the
CCSNe explosion mechanism.

To conclude, we have demonstrated that the interplay between neutrino
detectors and sub-Hz GW observatories will open a new path to studying
supernova neutrinos. Although several decades may pass before the first
results become available, the work of designing the next generation of
experiments is well under way, and we hope that our work will contribute to
its progress.

###### Acknowledgements.

We are grateful to Raffaella Margutti and Michele Zanolin for useful
discussions and comments. We also thank the anonymous referees for suggesting
valuable improvements. We acknowledge funding from the National Science
Foundation grant number PHY-2012195. MM was also supported by the Fermi
National Accelerator Laboratory (Fermilab) Award No. AWD00035045 during this
work. ZL acknowledges funding from the NSF Grant No. PHY-1554876, PHY 21-16686
and from DOE Scidac Grant DE-SC0018232.

## References

  * Raffelt and Seckel [1988] G. Raffelt and D. Seckel, Bounds on Exotic Particle Interactions from SN 1987a, Phys. Rev. Lett. 60, 1793 (1988).
  * Turner [1988] M. S. Turner, Axions from SN 1987a, Phys. Rev. Lett. 60, 1797 (1988).
  * Mayle  _et al._ [1988] R. Mayle, J. R. Wilson, J. R. Ellis, K. A. Olive, D. N. Schramm, and G. Steigman, Constraints on Axions from SN 1987a, Phys. Lett. B 203, 188 (1988).
  * Chang  _et al._ [2018] J. H. Chang, R. Essig, and S. D. McDermott, Supernova 1987A Constraints on Sub-GeV Dark Sectors, Millicharged Particles, the QCD Axion, and an Axion-like Particle, JHEP 09, 051, arXiv:1803.00993 [hep-ph] .
  * Ando  _et al._ [2005] S. Ando, J. F. Beacom, and H. Yüksel, Detection of neutrinos from supernovae in nearby galaxies, Phys. Rev. Lett. 95, 171101 (2005).
  * Kistler  _et al._ [2011] M. D. Kistler, H. Yuksel, S. Ando, J. F. Beacom, and Y. Suzuki, Core-Collapse Astrophysics with a Five-Megaton Neutrino Detector, Phys. Rev. D 83, 123008 (2011), arXiv:0810.1959 [astro-ph] .
  * Bisnovatyi-Kogan and Seidov [1982] G. S. Bisnovatyi-Kogan and Z. F. Seidov, Medium-Energy Neutrinos in the Universe, SOVAST 26, 132 (1982).
  * Krauss  _et al._ [1984] L. M. Krauss, S. L. Glashow, and D. N. Schramm, Anti-neutrinos Astronomy and Geophysics, Nature 310, 191 (1984).
  * Beacom [2010] J. F. Beacom, The Diffuse Supernova Neutrino Background, Ann. Rev. Nucl. Part. Sci. 60, 439 (2010), arXiv:1004.3311 [astro-ph.HE] .
  * Lunardini [2017] C. Lunardini, Diffuse Neutrino Flux from Supernovae, in  _Handbook of Supernovae_, edited by A. W. Alsabti and P. Murdin (2017) p. 1637.
  * De Gouvêa  _et al._ [2020] A. De Gouvêa, I. Martinez-Soler, Y. F. Perez-Gonzalez, and M. Sen, Fundamental physics with the diffuse supernova background neutrinos, Phys. Rev. D 102, 123012 (2020), arXiv:2007.13748 [hep-ph] .
  * Beacom and Vagins [2004] J. F. Beacom and M. R. Vagins, GADZOOKS! Anti-neutrino spectroscopy with large water Cherenkov detectors, Phys. Rev. Lett. 93, 171101 (2004), arXiv:hep-ph/0309300 .
  * Zhang  _et al._ [2015] H. Zhang _et al._ (Super-Kamiokande), Supernova Relic Neutrino Search with Neutron Tagging at Super-Kamiokande-IV, Astropart. Phys. 60, 41 (2015), arXiv:1311.3738 [hep-ex] .
  * An  _et al._ [2016] F. An _et al._ (JUNO), Neutrino Physics with JUNO, J. Phys. G 43, 030401 (2016), arXiv:1507.05613 [physics.ins-det] .
  * Abe  _et al._ [2018] K. Abe _et al._ (Hyper-Kamiokande), Hyper-Kamiokande Design Report, (2018), arXiv:1805.04163 [physics.ins-det] .
  * Askins  _et al._ [2020] M. Askins _et al._ (Theia), THEIA: an advanced optical neutrino detector, Eur. Phys. J. C 80, 416 (2020), arXiv:1911.03501 [physics.ins-det] .
  * Abi  _et al._ [2020] B. Abi, R. Acciarri, M. A. Acero, G. Adamov, D. Adams, M. Adinolfi, Z. Ahmad, J. Ahmed, T. Alion, S. A. Monsalve, _et al._ , Deep underground neutrino experiment (dune), far detector technical design report, volume ii: Dune physics, arXiv preprint arXiv:2002.03005 (2020).
  * Adams  _et al._ [2013] S. M. Adams, C. S. Kochanek, J. F. Beacom, M. R. Vagins, and K. Z. Stanek, Observing the Next Galactic Supernova, Astrophys. J. 778, 164 (2013), arXiv:1306.0559 [astro-ph.HE] .
  * Pagliaroli  _et al._ [2009] G. Pagliaroli, F. Vissani, E. Coccia, and W. Fulgione, Neutrinos from supernovae as a trigger for gravitational wave search, Phys. Rev. Lett. 103, 031102 (2009).
  * Abbott  _et al._ [2016] B. P. Abbott _et al._ (LIGO Scientific, Virgo), A First Targeted Search for Gravitational-Wave Bursts from Core-Collapse Supernovae in Data of First-Generation Laser Interferometer Detectors, Phys. Rev. D 94, 102001 (2016), arXiv:1605.01785 [gr-qc] .
  * Abe  _et al._ [2021] K. Abe _et al._ (Super-Kamiokande), Search for neutrinos in coincidence with gravitational wave events from the LIGO-Virgo O3a Observing Run with the Super-Kamiokande detector, Astrophys. J. 918, 78 (2021), arXiv:2104.09196 [astro-ph.HE] .
  * Szczepanczyk  _et al._ [2021] M. Szczepanczyk, J. Antelis, M. Benjamin, M. Cavaglia, D. Gondek-Rosinska, T. Hansen, S. Klimenko, M. Morales, C. Moreno, S. Mukherjee, _et al._ , Detecting and reconstructing gravitational waves from the next galactic core-collapse supernova in the advanced detector era, arXiv preprint arXiv:2104.06462 (2021).
  * Burrows and Hayes [1996] A. Burrows and J. Hayes, Pulsar recoil and gravitational radiation due to asymmetrical stellar collapse and explosion, Phys. Rev. Lett. 76, 352 (1996), arXiv:astro-ph/9511106 .
  * Mueller and Janka [1997] E. Mueller and H. T. Janka, Gravitational radiation from convective instabilities in Type II supernova explosions., AAP 317, 140 (1997).
  * Kotake  _et al._ [2009] K. Kotake, W. Iwakami, N. Ohnishi, and S. Yamada, Ray-tracing analysis of anisotropic neutrino radiation for estimating gravitational waves in core-collapse supernovae, The Astrophysical Journal 704, 951 (2009).
  * Muller  _et al._ [2012] E. Muller, H. T. Janka, and A. Wongwathanarat, Parametrized 3D models of neutrino-driven supernova explosions: Neutrino emission asymmetries and gravitational-wave signals, Astron. Astrophys. 537, A63 (2012), arXiv:1106.6301 [astro-ph.SR] .
  * Li  _et al._ [2018] J.-T. Li, G. M. Fuller, and C. T. Kishimoto, Neutrino Burst-Generated Gravitational Radiation From Collapsing Supermassive Stars, Phys. Rev. D 98, 023002 (2018), arXiv:1708.05292 [astro-ph.HE] .
  * Vartanyan and Burrows [2020] D. Vartanyan and A. Burrows, Gravitational Waves from Neutrino Emission Asymmetries in Core-collapse Supernovae, Astrophys. J. 901, 108 (2020), arXiv:2007.07261 [astro-ph.HE] .
  * Richardson  _et al._ [2021] C. Richardson, M. Zanolin, H. Andresen, M. J. Szczepańczyk, K. Gill, and A. Wongwathanarat, Modeling Core-Collapse Supernovae Gravitational-Wave Memory in Laser Interferometric Data, (2021), arXiv:2109.01582 [astro-ph.HE] .
  * Seto  _et al._ [2001] N. Seto, S. Kawamura, and T. Nakamura, Possibility of direct measurement of the acceleration of the universe using 0.1-Hz band laser interferometer gravitational wave antenna in space, Phys. Rev. Lett. 87, 221103 (2001), arXiv:astro-ph/0108011 .
  * Yagi and Seto [2011] K. Yagi and N. Seto, Detector configuration of decigo/bbo and identification of cosmological neutron-star binaries, Phys. Rev. D 83, 044011 (2011).
  * Sato  _et al._ [2017] S. Sato _et al._ , The status of DECIGO, J. Phys. Conf. Ser. 840, 012010 (2017).
  * Kawamura  _et al._ [2021] S. Kawamura, M. Ando, N. Seto, and e. Sato, Shuichi, Current status of space gravitational wave antenna DECIGO and B-DECIGO, Progress of Theoretical and Experimental Physics 2021, 05A105 (2021), arXiv:2006.13545 [gr-qc] .
  * Mukhopadhyay  _et al._ [2021] M. Mukhopadhyay, C. Cardona, and C. Lunardini, The neutrino gravitational memory from a core collapse supernova: phenomenology and physics potential, JCAP 07, 055, arXiv:2105.05862 [astro-ph.HE] .
  * Epstein [1978] R. Epstein, The Generation of Gravitational Radiation by Escaping Supernova Neutrinos, Astrophys. J. 223, 1037 (1978).
  * Turner [1978] M. S. Turner, Gravitational Radiation from Supernova Neutrino Bursts, Nature 274, 565 (1978).
  * Sathyaprakash and Schutz [2009] B. S. Sathyaprakash and B. F. Schutz, Physics, Astrophysics and Cosmology with Gravitational Waves, Living Rev. Rel. 12, 2 (2009), arXiv:0903.0338 [gr-qc] .
  * Moore  _et al._ [2015] C. J. Moore, R. H. Cole, and C. P. L. Berry, Gravitational-wave sensitivity curves, Class. Quant. Grav. 32, 015014 (2015), arXiv:1408.0740 [gr-qc] .
  * Jaranowski and Krolak [2000] P. Jaranowski and A. Krolak, Data analysis of gravitational wave signals from spinning neutron stars. 3. Detection statistics and computational requirements, Phys. Rev. D 61, 062001 (2000), arXiv:gr-qc/9901013 .
  * Abbott  _et al._ [2021] R. Abbott _et al._ (KAGRA, VIRGO, LIGO Scientific), All-sky search for long-duration gravitational-wave bursts in the third Advanced LIGO and Advanced Virgo run, Phys. Rev. D 104, 102001 (2021), arXiv:2107.13796 [gr-qc] .
  * Madau  _et al._ [1998] P. Madau, M. Della Valle, and N. Panagia, On the evolution of the cosmic supernova rates, Mon. Not. Roy. Astron. Soc. 297, 17 (1998), arXiv:astro-ph/9803284 .
  * Ando and Sato [2004] S. Ando and K. Sato, Relic neutrino background from cosmological supernovae, New J. Phys. 6, 170 (2004), arXiv:astro-ph/0410061 .
  * Daigne  _et al._ [2004] F. Daigne, K. A. Olive, E. Vangioni-Flam, J. Silk, and J. Audouze, Cosmic star formation, reionization, and constraints on global chemical evolution, Astrophys. J. 617, 693 (2004), arXiv:astro-ph/0405355 .
  * Phinney [2001] E. S. Phinney, A Practical Theorem on Gravitational Wave Backgrounds, arXiv e-prints , astro-ph/0108028 (2001), arXiv:astro-ph/0108028 [astro-ph] .
  * Buonanno  _et al._ [2005] A. Buonanno, G. Sigl, G. G. Raffelt, H.-T. Janka, and E. Muller, Stochastic gravitational wave background from cosmological supernovae, Phys. Rev. D 72, 084001 (2005), arXiv:astro-ph/0412277 .
  * Crocker  _et al._ [2017] K. Crocker, T. Prestegard, V. Mandic, T. Regimbau, K. Olive, and E. Vangioni, Systematic study of the stochastic gravitational-wave background due to stellar core collapse, Phys. Rev. D 95, 063015 (2017), arXiv:1701.02638 [astro-ph.CO] .
  * Finkel  _et al._ [2022] B. Finkel, H. Andresen, and V. Mandic, Stochastic gravitational-wave background from stellar core-collapse events, Phys. Rev. D 105, 063022 (2022), arXiv:2110.01478 [gr-qc] .
  * Keil  _et al._ [2003] M. T. Keil, G. G. Raffelt, and H.-T. Janka, Monte Carlo study of supernova neutrino spectra formation, Astrophys. J. 590, 971 (2003), arXiv:astro-ph/0208035 .
  * Sukhbold  _et al._ [2016] T. Sukhbold, T. Ertl, S. E. Woosley, J. M. Brown, and H. T. Janka, Core-Collapse Supernovae from 9 to 120 Solar Masses Based on Neutrino-powered Explosions, Astrophys. J. 821, 38 (2016), arXiv:1510.04643 [astro-ph.HE] .
  * Ertl  _et al._ [2016] T. Ertl, H. T. Janka, S. E. Woosley, T. Sukhbold, and M. Ugliano, A two-parameter criterion for classifying the explodability of massive stars by the neutrino-driven mechanism, Astrophys. J. 818, 124 (2016), arXiv:1503.07522 [astro-ph.SR] .
  * Kresse  _et al._ [2021] D. Kresse, T. Ertl, and H.-T. Janka, Stellar Collapse Diversity and the Diffuse Supernova Neutrino Background, Astrophys. J. 909, 169 (2021), arXiv:2010.04728 [astro-ph.HE] .
  * Hirata  _et al._ [1988] K. S. Hirata, T. Kajita, M. Koshiba, M. Nakahata, Y. Oyama, N. Sato, A. Suzuki, M. Takita, Y. Totsuka, T. Kifune, T. Suda, K. Takahashi, T. Tanimori, K. Miyano, M. Yamada, E. W. Beier, L. R. Feldscher, W. Frati, S. B. Kim, A. K. Mann, F. M. Newcomer, R. Van Berg, W. Zhang, and B. G. Cortez, Observation in the kamiokande-ii detector of the neutrino burst from supernova sn1987a, Phys. Rev. D 38, 448 (1988).
  * Abe  _et al._ [2011] K. Abe _et al._ , Letter of Intent: The Hyper-Kamiokande Experiment — Detector Design and Physics Potential —, (2011), arXiv:1109.3262 [hep-ex] .
  * Kyutoku and Kashiyama [2018] K. Kyutoku and K. Kashiyama, Detectability of thermal neutrinos from binary neutron-star mergers and implications for neutrino physics, Phys. Rev. D 97, 103001 (2018).
  * Strumia and Vissani [2003] A. Strumia and F. Vissani, Precise quasielastic neutrino/nucleon cross-section, Phys. Lett. B 564, 42 (2003), arXiv:astro-ph/0302055 .
  * Kunxian [2015] H. Kunxian, Ph.D. thesis, Kyoto University, available at http://www-he.scphys.kyoto-u.ac.jp/theses/doctor/khuang_dt.pdf (2015).
  * Nakamura  _et al._ [2016] K. Nakamura, S. Horiuchi, M. Tanaka, K. Hayama, T. Takiwaki, and K. Kotake, Multimessenger signals of long-term core-collapse supernova simulations: synergetic observation strategies, Mon. Not. Roy. Astron. Soc. 461, 3296 (2016), arXiv:1602.03028 [astro-ph.HE] .
  * Burrows  _et al._ [1992] A. Burrows, D. Klein, and R. Gandhi, The Future of supernova neutrino detection, Phys. Rev. D 45, 3361 (1992).
  * Kochanek  _et al._ [2017] C. S. Kochanek _et al._ , The All-Sky Automated Survey for Supernovae (ASAS-SN) Light Curve Server v1.0, Publ. Astron. Soc. Pac. 129, 104502 (2017), arXiv:1706.07060 [astro-ph.SR] .
  * Hiramatsu  _et al._ [2021] D. Hiramatsu _et al._ , The electron-capture origin of supernova 2018zd, Nature Astron. 5, 903 (2021), arXiv:2011.02176 [astro-ph.HE] .
  * Valenti  _et al._ [2016] S. Valenti _et al._ , The diversity of Type II supernova versus the similarity in their progenitors, Mon. Not. Roy. Astron. Soc. 459, 3939 (2016), arXiv:1603.08953 [astro-ph.SR] .
  * Spiro  _et al._ [2014] S. Spiro, A. Pastorello, M. Pumo, L. Zampieri, M. Turatto, S. Smartt, S. Benetti, E. Cappellaro, S. Valenti, I. Agnoletto, _et al._ , Low luminosity type ii supernovae–ii. pointing towards moderate mass precursors, Monthly Notices of the Royal Astronomical Society 439, 2873 (2014).
  * Hosseinzadeh  _et al._ [2018] G. Hosseinzadeh _et al._ , Short-Lived Circumstellar Interaction in the Low-Luminosity Type IIP SN 2016bkv, Astrophys. J. 861, 63 (2018), arXiv:1801.00015 [astro-ph.HE] .
  * Woosley  _et al._ [2002] S. E. Woosley, A. Heger, and T. A. Weaver, The evolution and explosion of massive stars, Rev. Mod. Phys. 74, 1015 (2002).
  * Sumiyoshi  _et al._ [2006] K. Sumiyoshi, S. Yamada, H. Suzuki, and S. Chiba, Neutrino signals from the formation of black hole: A probe of equation of state of dense matter, Phys. Rev. Lett. 97, 091101 (2006), arXiv:astro-ph/0608509 .
  * O’Connor and Ott [2011] E. O’Connor and C. D. Ott, Black Hole Formation in Failing Core-Collapse Supernovae, Astrophys. J. 730, 70 (2011), arXiv:1010.5550 [astro-ph.HE] .
  * LSST Science Collaboration [2009] LSST Science Collaboration, LSST Science Book, Version 2.0, arXiv e-prints , arXiv:0912.0201 (2009), arXiv:0912.0201 [astro-ph.IM] .


# Resetting Metadynamics

Ofir Blumer  [  Shlomi Reuveni  [  Barak Hirshberg<EMAIL_ADDRESS>[

###### Abstract

Metadynamics is a powerful method to accelerate molecular dynamics
simulations, but its efficiency critically depends on the identification of
collective variables that capture the slow modes of the process.
Unfortunately, collective variables are usually not known a priori and finding
them can be very challenging. We recently presented a collective variables-
free approach to enhanced sampling using stochastic resetting. Here, we
combine the two methods for the first time, showing that it can lead to
greater acceleration than either of them separately. We also demonstrate that
resetting Metadynamics simulations performed with suboptimal collective
variables can lead to speedups comparable with those obtained with optimal
collective variables. Therefore, the application of stochastic resetting can
be an alternative to the challenging task of improving suboptimal collective
variables, at almost no additional computational cost. Finally, we propose a
method to extract unbiased mean first-passage times from Metadynamics
simulations with resetting, resulting in improved tradeoff between speedup and
accuracy. This work opens the way for combining stochastic resetting with
other enhanced sampling methods to accelerate a broad range of molecular
simulations.

###### keywords:

American Chemical Society, LaTeX

TAU] School of Chemistry, Tel Aviv University, Tel Aviv 6997801, Israel. TAU]
School of Chemistry, Tel Aviv University, Tel Aviv 6997801, Israel. [Sackler]
The Center for Computational Molecular and Materials Science, Tel Aviv
University, Tel Aviv 6997801, Israel. [BioSoft] The Center for Physics and
Chemistry of Living Systems, Tel Aviv University, Tel Aviv 6997801, Israel.
TAU] School of Chemistry, Tel Aviv University, Tel Aviv 6997801, Israel.
[Sackler] The Center for Computational Molecular and Materials Science, Tel
Aviv University, Tel Aviv 6997801, Israel. IR,NMR,UV

##  1 Introduction

Molecular dynamics (MD) simulations provide valuable insights into the
dynamics of complex chemical and physical systems. They are a powerful tool
but, due to their atomic spatial and temporal resolution, they cannot be
applied to processes that are longer than a few microseconds, such as protein
folding and crystal nucleation1, 2. Different methods have been developed in
order to overcome this timescale problem, such as umbrella sampling3, 4,
replica-exchange5, free energy dynamics6, 7, Metadynamics (MetaD)8, 9, 10, 11,
12, On-the-fly probability enhanced sampling (OPES)13, 14, 15, 16, and many
others. In this paper, we will focus on MetaD, which relies on identifying
efficient collective variables (CVs), capturing the slow modes of the process,
and introducing an external bias potential to enhance the sampling of phase
space along them. The ability of Metadynamics to accelerate simulations
crucially depends on the quality of the CVs17, 12. An optimal CV is capable of
distinguishing between metastable states of interest as well as describing
their interconversion dynamics. 18, 19 A suboptimal CV can lead to hysteresis,
and poor inference of the unbiased free energy surface or kinetics8, 20, 17,
1, 21.

Very recently, we developed a CV-free approach for enhanced sampling, based on
stochastic resetting (SR)22. Resetting is the procedure of stopping stochastic
processes, at random or fixed time intervals, and restarting them using
independent and identically distributed initial conditions. It has received
much attention recently23, 24, since it is able to expedite various processes
ranging from randomized computer algorithms25, 26, 27 and service in queuing
systems28, to first-passage and search processes29, 30, 31, 32, 33, 34, 35,
36, 37, 38. We demonstrated the power of SR in enhanced sampling of MD
simulations, showing it can lead to speedups of up to an order of magnitude in
simple model systems and a molecular system22. Moreover, we developed a method
to infer the unbiased kinetics, in the absence of SR, from simulations with
SR.

Resetting is an appealing method due to its extreme simplicity: it can be
trivially implemented in MD codes, and no CVs are required. Since finding good
CVs in complex condensed phase systems is a difficult challenge, it is a
potentially significant advantage. Furthermore, unlike other methods, which
continuously add energy to the system, SR does not change the dynamics between
resetting events. However, acceleration by SR is not guaranteed. A sufficient
condition is that the distribution of transition times (also called first-
passage times, FPTs) of the corresponding process is wide enough, i.e., it has
a standard deviation that is greater than its mean39. Many systems, including
the models and molecular example discussed in our previous work22, fulfill
this criterion but there are also counter examples.

The complementary advantages and limitations of MetaD and SR raise an
important question: can we combine them to obtain the best of both worlds?
Since SR can be applied to any random process, we can restart MetaD
simulations as previously done for unbiased simulations. This observation
opens the door for using SR as a complementary tool to existing enhanced
sampling procedures.

In this paper, we combine SR with MetaD for the first time. In one model
system, we show that this approach leads to greater acceleration than any of
them independently, even in comparison to using the optimal CV in MetaD
simulations. In another model system, we restart MetaD simulations, performed
with suboptimal CVs, and obtain accelerations similar to those achieved using
the optimal CV. This result suggests that a straightforward application of SR
can be an alternative to the challenging task of improving a suboptimal CV. We
then demonstrate this for transitions between meta-stable states of alanine
tetrapeptide. Lastly, we develop a procedure to infer the unbiased kinetics
from simulations combining SR and MetaD, showing an improved tradeoff between
speedup and accuracy in comparison to MetaD simulations alone.

##  2 Methods

The simulations of the model systems were performed using the Large-scale
Atomic/Molecular Massively Parallel Simulator40 (LAMMPS), with MetaD
introduced using PLUMED 2.7.141, 42, 43. Alanine tetrapeptide was simulated
using GROMACS 2019.644 and the same PLUMED version. All simulations were
carried out in the NVT ensemble at a temperature of $300\,K$, and $10^{4}$
trajectories were sampled for each case. Full simulation details are given in
the SI.

Initial velocities were sampled from the Maxwell-Boltzmann distribution, while
initial positions were fixed (equivalent to sampling from a delta function
positions distribution). We defined the FPT as the earliest instance in which
a certain criterion was met, as specified below for each system. The COMMITTOR
command in PLUMED was used for testing this criterion and stopping the
simulations when it was fulfilled. For most simulations with SR, we sampled
the time intervals between resetting events from an exponential distribution
with a fixed resetting rate (Poisson resetting) using Python. Simulations
designated for kinetics inference used constant time intervals between
resetting events (Sharp resetting). If a first-passage event did not occur
prior to the next resetting time, the simulation was restarted. We stress that
we continue tallying the overall time until a first-passage occurred,
regardless of the number of resetting events in between. Finally, for
simulations combining SR with MetaD, we emphasize that the MetaD bias
potential was zeroed after each resetting event.

##  3 Results and discussion

###  3.1 A two wells model

We begin by combining SR with MetaD simulations for the first time. A simple
model system is considered first, where the optimal CV is well defined. For
this model, combining SR with MetaD leads to greater speedups than MetaD
independently, even when biasing with the optimal CV. The model is shown in
Figure 1(a). It is a two dimensional harmonic trap, divided into two states
centered at $(x=\pm 3,y=0)\,\AA$ by a barrier of $5\,k_{B}T$ for all y values.

The harmonic trap is soft, such that a particle in one of the wells can easily
travel about $50\,\AA$ away from the central barrier. The exact parameters of
the potential are given in the SI. They were chosen such that the unbiased
mean FPT (MFPT) between the wells is long ($7.5\,ns$), but can still be
sampled in unbiased simulations. The optimal CV is the x-coordinate. We follow
the trajectories of a particle that was initialized at the right minimum and
define the FPT criterion as arriving to the left minimum ($x\leq-3\,\AA$).

For comparison, we first performed SR on standard MD simulations using
different resetting rates. The obtained speedups are shown in Figure 1(b). The
speedup increases with the resetting rate, reaches a maximum value of $\sim
4$, which is obtained at a rate of $50\,ns^{-1}$, and decreases for higher
rates. This non-monotonic trend can be understood since in the limit of high
resetting rates all trajectories are restarted before a transition can occur,
and the speedup drops to zero. The resetting rate which leads to maximum
speedup will be referred to as the optimal resetting rate in this paper.

Next, we performed MetaD simulations without SR using the x-coordinate as a CV
and varying the bias deposition rates. Other bias parameters are given in the
SI. The results are shown as green squares in Figure 1(c). The speedup
increases with the bias rate, with a value of $\sim 20$ attained for a bias
rate of $10^{4}\,ns^{-1}$ (every 100 simulation steps). It is evident that
MetaD leads to larger acceleration than SR for this system, giving speedups
that are greater by a factor of $\sim 5$.

We then combined SR with MetaD and found that even greater speedups are
obtained (Figure 1(c), orange triangles). How the combination of resetting and
MetaD is done in practice is shown in Figure 1(d), using the highest bias
deposition rate ($10^{4}\,ns^{-1}$) as an example. The green square in Figure
1(d) shows the speedup obtained with MetaD and no resetting. Then, for that
given bias deposition rate, we add SR at increasing rates and evaluate the
resulting speedup. We again stress that the MetaD bias is zeroed at every
resetting event. We observe the same qualitative behavior seen in Figure 1(b),
with the speedup increasing until some optimal resetting rate, highlighted
with an orange triangle. We repeat this procedure for all bias rates, and
present the optimal speedup by orange triangles in Figure 1(c). Combining SR
with MetaD gave additional acceleration for all bias deposition rates, with a
maximal speedup of $\sim 130$ at a bias rate of $10^{4}\,ns^{-1}$. The
corresponding optimal resetting rate was found to be $125\,ns^{-1}$, which is
significantly slower than the bias deposition. The fact that SR can further
accelerate MetaD simulations, even when performed with optimal CVs, is the
first key result of this paper.

Figure 1: (a) The two wells model. The star marks the initial position. (b)
Speedups obtained for SR simulations at different resetting rates. (c)
Speedups obtained for MetaD simulations at different bias deposition rates
(green squares), and for combined MetaD + SR simulations at an optimal
resetting rate (orange triangles). (d) Speedups obtained for MetaD + SR
simulations with a bias deposition rate of $10^{4}\,ns^{-1}$ and different
resetting rates.

MetaD practitioners might wonder: 1) How to tell whether SR will accelerate my
simulations?, and 2) How to identify the optimal resetting rate and what would
be the resulting speedup? Next, we show that both questions can be answered at
almost no additional cost, assuming some MetaD trajectories are already
available. To answer the first question, we showed in a recent paper 22 that a
sufficient condition for acceleration of MD simulations by SR is that the
ratio of the standard deviation to the MFPT (the coefficient of variation,
COV) is greater than one 39. Introducing a small resetting rate is then
guaranteed to lead to speedup. We stress that this condition holds also for
resetting MetaD simulations, with the added benefit that enhanced sampling
generates more transitions, and thus gives a much more reliable estimation of
the COV compared to unbiased MD simulations. Moreover, if SR does not
accelerate the unbiased simulations significantly, it does not mean that it
will not do so for MetaD simulations, since biasing alters the FPT
distribution significantly and, consequently, may also change the COV.

Figure 2(a) shows the COV of MetaD simulations (without resetting) as a
function of the bias deposition rate, for the two wells model (Figure 1(a)).
The COV shows non-monotonic behavior with the bias deposition rate. It starts
from a value of $1.05$ without resetting, drops to $0.80$ at an intermediate
bias deposition rate and increases up to a value of $1.44$ at the highest
biasing rate. This shows that MetaD can increase the COV significantly
(leading to much higher speedups).

As for the second question, estimating the optimal resetting rate and the
resulting speedup is also straightforward for MetaD simulations. The MFPT
under a resetting rate $r$ can be estimated using a simple algebraic equation
45,

$\langle\tau\rangle_{r}=\frac{1-\tilde{f}(r)}{r\tilde{f}(r)}.$ (1)

In Eq. 1, $\tilde{f}(r)$ is the Laplace transform of the FPT distribution for
the MetaD simulations, and $\langle\tau\rangle_{r}$ is the MFPT for MetaD
simulations with SR at resetting rate $r$. The Laplace transform is evaluated
as

$\tilde{f}(r)=\langle
e^{-r\tau}\rangle\simeq\frac{1}{N}\sum_{j=1}^{N}e^{-r\tau_{j}},$ (2)

where $N$ is the number of MetaD trajectories, and $\tau_{j}$ is the FPT
obtained from trajectory $j$. Figure 2(b) shows the additional speedups, over
MetaD without resetting, estimated using Eq. 1 (dotted lines). They are
plotted as a function of the resetting rate for the two bias deposition rates
highlighted with colored circles in Figure 2(a). It is evident that the
estimations match results obtained from simulations (full circles). While the
full FPT distribution is required for an exact description of the behavior
under SR, we previously showed22 that as few as a hundred samples are
sufficient for estimating the optimal resetting rate.

Figure 2: (a) COV of MetaD simulations with no SR. (b) Additional speedups
against resetting rate for MetaD simulations with bias deposition rate of
$10^{4}\,ns^{-1}$ (blue) or $20\,ns^{-1}$ (green), for the two wells model.
Full circles present results obtained from simulations, while dotted lines
present evaluations based on the FPT distribution with MetaD and no SR and
using Eq. 1. The dashed gray line indicates no additional speedup.

Surprisingly, for intermediate bias deposition rates, some additional speedup
is gained even though the COV without resetting is smaller than 1. This can be
seen in Figure 2(b), which shows an optimal speedup of $\sim 2$ for a bias
deposition rate of $20\,ns^{-1}$, even though its COV without SR is only
$0.80$. How can it be? A COV $<1$ indicates that introducing a small resetting
rate will increase the MFPT. In other words, the COV only indicates the
initial slope of the speedup curve as a function of the resetting rate.
Interestingly, non-trivial cases where small resetting rates decelerate the
process but larger ones accelerate it are also possible, as can be seen for
the green curve in the inset of Figure 2(b). For comparison, we also give the
results for a bias deposition rate of $10^{4}\,ns^{-1}$, for which the COV
without SR is $>1$, and the initial slope of the speedup with respect to the
resetting rate is positive. The results show that, whether the value of the
COV is greater or smaller than one, it is worthwhile to estimate the speedup
using Eq. 1.

To conclude this example, we combined SR with MetaD for the first time,
leading to faster acceleration than either approach separately. This is
demonstrated even for a system for which the optimal CV is known. Since MetaD
simulations already enhance the sampling of the underlying process, it is
significantly easier to evaluate their COV than for unbiased simulations. If
the COV is larger than 1, then MetaD simulations are guaranteed to be further
accelerated by SR, and Eq. 1 can be used to easily estimate by how much. If
not, they may still be accelerated, which can be easily checked using Eq. 1.

###  3.2 The modified Faradjian-Elber potential

As a second example, we consider a modified version of the two dimensional
potential introduced by Faradjian and Elber when developing the milestoning
enhanced sampling method 46. The potential is shown in Figure 3(a), and full
details are given in the SI. It is also composed of two symmetric wells, with
minima at $(x=\pm 3,y=0)\,\AA$ that are separated by a Gaussian barrier at
$x=0\,\AA$. The barrier is higher than the first example, $12\,k_{B}T$ for
most $y$ values, but has a narrow saddle, only $3\,k_{B}T$ high, around
$y=0\,\AA$. Figure 3(b) shows slices along the x-axis at $y=0,25\,\AA$, as
well as the effective potential integrated over the entire y-axis.

We follow the trajectories of a particle that was initialized at the right
minimum and define the FPT criterion as crossing the barrier and reaching
$x<-1\,\AA$. For this model, we find the same MFPT as in the two wells model.
Employing SR on unbiased simulations gave an optimal speedup of $\sim 15$ at a
resetting rate of $200\,ns^{-1}$. As in the two wells model, MetaD simulations
gave higher speedups than SR, with a speedup of $\sim 212$ when using the
optimal CV, the x-coordinate, at a bias deposition rate of $10^{4}\,ns^{-1}$.
Using this optimal CV and rate, combining SR with MetaD did not lead to
further acceleration of the simulations.

However, in most real systems, the optimal CV is not known, and suboptimal CVs
are almost always used 17. To test the efficiency of SR in such cases, we
gradually reduce the quality of the CV by rotating it. The green squares in
Figure 3(c) show the speedup obtained as a function of the sine of the angle
$\theta$ between the CV and the x-axis, which serves as a measure of the
deviation from the optimal CV. The degradation in the quality of the CV leads
to a decrease in the MetaD speedup with almost no acceleration at an angle of
$24^{\circ}$. However, combining SR with MetaD recovers almost all of the
speedup of the optimal CV, despite the use of suboptimal CVs . This is shown
by the orange triangles in Figure 3(c). Optimizing CVs for condensed phase
systems remains a difficult challenge 9, 47. Our results suggest that SR may
serve as an alternative, or complementary method, to improving CVs. Instead of
using sophisticated algorithms to find better CVs 47, 48, 49, 50, 51, 52, one
can use SR to obtain a similar speedup at a much lower cost. This is the
second key result of this paper.

Figure 3: (a) The modified Faradjian-Elber potential. The star marks the
initial position. (b) The integrated projection of the potential on the x-axis
(blue), and the potential cross section at $y=0,25\,\AA$ (dotted and dashed
black lines, respectively). (c) Speedups obtained for a bias deposition rate
of $10^{4}\,ns^{-1}$ using suboptimal, rotated CVs, in the MetaD simulations
with no SR (green squares) and with optimal SR (orange triangles). The angle
between the CV and the x-axis is denoted by $\theta$.

###  3.3 Alanine tetrapeptide

As a final example, we demonstrate the capabilities of combining MetaD with SR
on a molecular system, alanine tetrapeptide. We focus on two of its
conformers: a “folded” and “unfolded” states, shown in Figure 4(a). Six
dihedral angles serve as important degrees of freedoms, with $\phi_{3}$ being
the slowest one13. Figure 4(b) shows the free-energy surface along $\phi_{3}$,
which has two minima separated by an energy barrier of $\sim 15\,k_{B}T$.
Transitions in unbiased simulations from the unfolded state (upper
configuration in panel (a), left basin in panel (b)) to the folded one (lower
configuration in panel (a), right basin in panel (b)) have an estimated MFPT
of $\sim 5.6\,\mu s$ (see the SI for more details).

To improve the sampling, we performed MetaD simulations using three different
CVs: the angle $\phi_{3}$ serves as the optimal CV, and two adjacent angles,
$\phi_{2}$ and $\psi_{3}$, serve as suboptimal ones. The two dimensional free
energy surfaces as a function of all CVs are presented in panels (c) and (d)
of Figure 4. They show that $\psi_{3}$ has some overlap and does not separate
the two states as well as $\phi_{3}$, while there is almost no separation of
the states in $\phi_{2}$. The simulations were initialized from a fixed,
unfolded configuration, marked with stars in panels (b) to (d). The first-
passage criterion ($0.5<\phi_{3}<1.5\,rad$) is also marked in these panels,
with vertical dashed lines.

Figure 4(e) shows the speedup of MetaD simulations using different protocols,
without SR (green) and with optimal SR (orange). COV values for simulations
with no SR are given in panel (f). As expected, using $\phi_{3}$ as a CV gives
the greatest speedup, and a COV of $\sim 0.3$ for which the optimal resetting
rate is $r^{*}=0$. Thus, there is no benefit from SR in this case. Suboptimal
CVs show similar behavior to that observed for the Faradjian-Elber potential,
but for a realistic system: The speedups obtained for MetaD without SR
decrease when using bad CVs and the COV values increase above one. Namely,
while MetaD simulations using $\phi_{3}$ as the CV gives more than four orders
of magnitude speedup, simulations using $\psi_{3}$ and $\phi_{2}$ lead to
accelerations by factors of only $\sim 580$ and $\sim 4$, respectively.
Concurrently, the COV of $\psi_{3}$ is $\sim 1.24$ while for the worst CV,
$\phi_{2}$, it is $\sim 3$. As a result, SR becomes more effective the poorer
the CV is, giving an additional speedup of $\sim 133$ over MetaD when using
$\phi_{2}$ as a CV.

Figure 4: (a) Two conformers of alanine tetrapeptide. The white, gray, blue,
and red balls represent hydrogen, carbon, nitrogen, and oxygen atoms,
respectively. Free-energy surfaces of alanine tetrapeptide along (b)
$\phi_{3}$, (c) $\phi_{3},\phi_{2}$ and (d) $\phi_{3},\psi_{3}$. (e) Speedup
of MetaD simulations, without SR (green) and with optimal SR (orange). (f) COV
of MetaD simulations without SR for different CVs.

###  3.4 Kinetics inference

To conclude the paper, we demonstrate that SR can improve the inference of
unbiased kinetics from MetaD simulations, using the case of alanine
tetrapeptide as an example. The unbiased FPT distribution can be extracted
from MetaD simulations with no SR through a procedure known as infrequent
MetaD (iMetaD).2 In this method, the MFPT is obtained by rescaling the FPT of
each MetaD trajectory by an acceleration factor that depends exponentially on
the external bias (see Eq. S5 in the SI). iMetaD assumes that no bias is
deposited near the transition state, and that none of the basins are over-
filled. When this assumption is valid, the distribution of the rescaled FPTs
matches the unbiased distribution. However, the assumption does not hold for
suboptimal CVs or high bias deposition rates 1, which result in over-
deposition. Due to the exponential dependence of the acceleration factor on
the bias, trajectories exhibiting over-deposition result in extremely large
acceleration factors. They contribute unrealisticly long FPTs to the rescaled
distribution, shifting the obtained MFPT from the true value by orders of
magnitude.

The inference can be improved, even with suboptimal CVs, by decreasing the
bias deposition rate, resulting in a tradeoff between speedup and accuracy.
This tradeoff is demonstrated in Figure 5. It shows (green squares) the error
in the estimation of the unbiased MFPT as a function of speedup, for iMetaD
simulations of alanine tetrapeptide biasing the $\psi_{3}$ angle, which is a
suboptimal CV. The prediction error is defined as
$\left|\langle\tau\rangle_{true}-\langle\tau\rangle_{est}\right|/\langle\tau\rangle_{true}$
where $\langle\tau\rangle_{true}$ is the true unbiased MFPT and
$\langle\tau\rangle_{est}$ is the estimated MFPT.

Next, we demonstrate that resetting can give a better tradeoff, reducing the
error for all speedups, as shown by the orange triangles in Figure 5. These
results were obtained by adding SR at different resetting rates to the iMetaD
simulations at the highest bias deposition rate (highlighted by a circle).

Full details explaining how to infer the unbiased MFPT from combined SR and
iMetaD simulations are given in the SI. Here, we briefly provide only the key
ingredients and underlying intuition. We note that between resetting events,
the trajectories are standard iMetaD simulations. Moreover, due to SR, the
short trajectories between restarts, can be treated as independent from one
another (recall that restart also zeros previous bias). As a result, we can
use the standard iMetaD rescaling procedure on each short trajectory, and then
evaluate the unbiased survival probability at short times. For short enough
times, even with suboptimal CVs, we will avoid over-deposition and get a good
estimate of the unbiased survival. Finally, we assume that the survival
probability decays exponentially, as commonly done in iMetaD 1, and obtain an
estimate of the unbiased MFPT from its slope. The quality of the linear fit
can be used to assess the reliability of the predicted MFPT, similar to the
Kolmogorov–Smirnov test used by Salvalaglio et. al in standard iMetaD
simulations.1, 53, 54

Figure 5: The error in the prediction of the unbiased MFPT as a function of
speedup for iMetaD simulations at different bias deposition rates (green
squares) and iMetaD simulations with SR at different resetting rates (orange
triangles).

Our results show that, for alanine tetrapetide, applying SR to iMetaD
simulations and using the proposed inference procedure gives a better tradeoff
than decreasing the bias deposition rate. In choosing between the two,
practitioners of iMetaD might consider a simple question: given a fixed
simulation time, would the error obtained be lower in a single long iMetaD
trajectory or on a series of shorter iMetaD trajectories with SR? For
suboptimal CVs, the acceleration factor becomes increasingly unreliable with
time, due to over-deposition of the external bias. In that case, short
trajectories with minimal bias are preferred, and SR may be more favorable.

###  3.5 Conclusions

We combine SR with MetaD simulations for the first time. We show that
resetting can further accelerate MetaD simulations, even when the optimal CV
is used. In practice, the optimal CV is almost never known and suboptimal CVs
are employed. We provide examples in which adding SR to MetaD simulations
performed with poor CVs, leads to speedups comparable to using the optimal CV.
This suggests that resetting may serve as an alternative to the challenging
task of improving suboptimal CVs using sophisticated algorithms.

Testing whether SR can accelerate simulations is very easy. Given a small
number ($\sim 100$) of short MetaD trajectories, showing one transition each,
we can estimate the COV and find whether SR would further accelerate the
simulations, and by how much, using Eq. 1. Resetting can be of benefit even
when it does not provide additional acceleration on top of the one attained by
MetaD. We demonstrate that SR can improve the inference of the unbiased
kinetics from iMetaD simulations performed with suboptimal CVs, giving a
better tradeoff between speedup and accuracy for alanine tetrapeptide.

Finally, we conjecture that benefits coming from combining MetaD and SR are
not limited to the examples presented herein, and are much more general. The
reason is that MetaD, and similar methods, flatten the free energy surface.
Previous work has shown that SR is particularly efficient for flat landscapes
33, 55, 56, with the extreme case being free diffusion 57. Thus, starting from
an arbitrary free energy surface, MetaD changes it to one that is more
amenable to acceleration by SR. Future method development would likely harness
the power of this important observation.

Barak Hirshberg acknowledges support by the Israel Science Foundation (grants
No. 1037/22 and 1312/22) and the Pazy Foundation of the IAEC-UPBC (grant No.
415-2023). Shlomi Reuveni acknowledges support from the Israel Science
Foundation (grant No. 394/19). This project has received funding from the
European Research Council (ERC) under the European Union’s Horizon 2020
research and innovation program (grant agreement No. 947731).

See the SI for full details of the computational setup and of the procedure to
infer the unbiased kinetics from Metadynamics simulations with SR. Example
input files are available on the GitHub repository 58.

## References

  * Salvalaglio et al. 2014 Salvalaglio, M.; Tiwary, P.; Parrinello, M. Assessing the Reliability of the Dynamics Reconstructed from Metadynamics. _Journal of Chemical Theory and Computation_ 2014, _10_ , 1420–1425
  * Tiwary and Parrinello 2013 Tiwary, P.; Parrinello, M. From Metadynamics to Dynamics. _Physical Review Letters_ 2013, _111_ , 230602
  * Torrie and Valleau 1977 Torrie, G.; Valleau, J. Nonphysical Sampling Distributions in Monte Carlo Free-Energy Estimation: Umbrella Sampling. _J. Comput. Phys._ 1977, _23_ , 187–199
  * Kästner 2011 Kästner, J. Umbrella Sampling. _Wiley Interdiscip. Rev. Comput. Mol. Sci._ 2011, _1_ , 932–942
  * Sugita and Okamoto 1999 Sugita, Y.; Okamoto, Y. Replica-exchange molecular dynamics method for protein folding. _Chemical Physics Letters_ 1999, _314_ , 141–151
  * Rosso and Tuckerman 2002 Rosso, L.; Tuckerman, M. E. An Adiabatic Molecular Dynamics Method for the Calculation of Free Energy Profiles. _Mol. Simul._ 2002, _28_ , 91–112
  * Rosso et al. 2002 Rosso, L.; Mináry, P.; Zhu, Z.; Tuckerman, M. E. On the Use of the Adiabatic Molecular Dynamics Technique in the Calculation of Free Energy Profiles. _J. Chem. Phys._ 2002, _116_ , 4389–4402
  * Barducci et al. 2011 Barducci, A.; Bonomi, M.; Parrinello, M. Metadynamics. _Wiley Interdiscip. Rev. Comput. Mol. Sci._ 2011, _1_ , 826–843
  * Valsson et al. 2016 Valsson, O.; Tiwary, P.; Parrinello, M. Enhancing Important Fluctuations: Rare Events and Metadynamics from a Conceptual Viewpoint. _Annual Review of Physical Chemistry_ 2016, _67_ , 159–184
  * Sutto et al. 2012 Sutto, L.; Marsili, S.; Gervasio, F. L. New Advances in Metadynamics. _Wiley Interdiscip. Rev. Comput. Mol. Sci._ 2012, _2_ , 771–779
  * Barducci et al. 2008 Barducci, A.; Bussi, G.; Parrinello, M. Well-Tempered Metadynamics: A Smoothly Converging and Tunable Free-Energy Method. _Physical Review Letters_ 2008, _100_ , 020603
  * Bussi and Laio 2020 Bussi, G.; Laio, A. Using metadynamics to explore complex free-energy landscapes. _Nature Reviews Physics_ 2020, _2_ , 200–212
  * Invernizzi and Parrinello 2020 Invernizzi, M.; Parrinello, M. Rethinking Metadynamics: From Bias Potentials to Probability Distributions. _J. Phys. Chem. Lett._ 2020, _11_ , 2731–2736
  * Invernizzi et al. 2020 Invernizzi, M.; Piaggi, P. M.; Parrinello, M. Unified Approach to Enhanced Sampling. _Phys. Rev. X_ 2020, _10_ , 041034
  * Invernizzi 2021 Invernizzi, M. OPES: On-the-fly Probability Enhanced Sampling method. _Nuovo Cimento C_ 2021, _44_ , 1–4
  * Invernizzi and Parrinello 2022 Invernizzi, M.; Parrinello, M. Exploration vs Convergence Speed in Adaptive-Bias Enhanced Sampling. _Journal of Chemical Theory and Computation_ 2022, _18_ , 3988–3996
  * Invernizzi and Parrinello 2019 Invernizzi, M.; Parrinello, M. Making the Best of a Bad Situation: A Multiscale Approach to Free Energy Calculation. _J. Chem. Theory. Comput._ 2019, _15_ , 2187–2194
  * Demuynck et al. 2018 Demuynck, R.; Wieme, J.; Rogge, S. M. J.; Dedecker, K. D.; Vanduyfhuys, L.; Waroquier, M.; Van Speybroeck, V. Protocol for Identifying Accurate Collective Variables in Enhanced Molecular Dynamics Simulations for the Description of Structural Transformations in Flexible Metal–Organic Frameworks. _Journal of Chemical Theory and Computation_ 2018, _14_ , 5511–5526
  * Peters 2017 Peters, B. In _Reaction Rate Theory and Rare Events Simulations_ ; Peters, B., Ed.; Elsevier: Amsterdam, 2017; pp 539–571
  * Bešker and Gervasio 2012 Bešker, N.; Gervasio, F. L. In _Computational Drug Discovery and Design_ ; Baron, R., Ed.; Methods in Molecular Biology; Springer: New York, NY, 2012; pp 501–513
  * Ray et al. 2022 Ray, D.; Ansari, N.; Rizzi, V.; Invernizzi, M.; Parrinello, M. Rare Event Kinetics from Adaptive Bias Enhanced Sampling. _Journal of Chemical Theory and Computation_ 2022, _18_ , 6500–6509
  * Blumer et al. 2022 Blumer, O.; Reuveni, S.; Hirshberg, B. Stochastic Resetting for Enhanced Sampling. _The Journal of Physical Chemistry Letters_ 2022, _13_ , 11230–11236
  * Evans et al. 2020 Evans, M. R.; Majumdar, S. N.; Schehr, G. Stochastic Resetting and Applications. _J. Phys. A : Math. Theor._ 2020, _53_ , 193001
  * Kundu and Reuveni 2022 Kundu, A.; Reuveni, S. Stochastic Resetting: Theory and Applications [special issue]. _J. Phys. A : Math. Theor._ 2022,
  * Luby et al. 1993 Luby, M.; Sinclair, A.; Zuckerman, D. Optimal Speedup of Las Vegas Algorithms. _Inf. Process. Lett._ 1993, _47_ , 173–180
  * Gomes et al. 1998 Gomes, C.; Selman, B.; Kautz, H. _15th AAAI_ ; The AAAI Press: Madison, WI, 1998; p 431–437
  * Montanari and Zecchina 2002 Montanari, A.; Zecchina, R. Optimizing Searches via Rare Events. _Phys. Rev. Lett._ 2002, _88_ , 178701
  * Bonomo et al. 2022 Bonomo, O. L.; Pal, A.; Reuveni, S. Mitigating Long Queues and Waiting Times with Service Resetting. _PNAS Nexus_ 2022, _1_ , pgac070
  * Bressloff 2020 Bressloff, P. C. Queueing Theory of Search Processes with Stochastic Resetting. _Phys. Rev. E_ 2020, _102_ , 032109
  * Kuśmierz and Gudowska-Nowak 2015 Kuśmierz, L.; Gudowska-Nowak, E. Optimal First-Arrival Times in Lévy Flights with Resetting. _Phys. Rev. E_ 2015, _92_ , 052127
  * Bhat et al. 2016 Bhat, U.; Bacco, C. D.; Redner, S. Stochastic Search with Poisson and Deterministic Resetting. _J. Stat. Mech._ 2016, _8_ , 083401
  * Chechkin and Sokolov 2018 Chechkin, A.; Sokolov, I. Random Search with Resetting: A Unified Renewal Approach. _Phys. Rev. Lett._ 2018, _121_ , 050601
  * Ray et al. 2019 Ray, S.; Mondal, D.; Reuveni, S. Péclet Number Governs Transition to Acceleratory Restart in Drift-Diffusion. 2019; http://arxiv.org/abs/1811.08239, Accessed: November 10, 2022
  * Robin et al. 2019 Robin, T.; Hadany, L.; Urbakh, M. Random Search with Resetting as a Strategy for Optimal Pollination. _Phys. Rev. E_ 2019, _99_ , 052119
  * Evans and Majumdar 2018 Evans, M. R.; Majumdar, S. N. Run and Tumble Particle under Resetting: a Renewal Approach. _J. Phys. A : Math. Theor._ 2018, _51_ , 475003
  * Pal et al. 2020 Pal, A.; Kuśmierz, L.; Reuveni, S. Search with Home Returns Provides Advantage under High Uncertainty. _Phys. Rev. Res._ 2020, _2_ , 043174
  * Bodrova and Sokolov 2020 Bodrova, A. S.; Sokolov, I. M. Resetting Processes with Noninstantaneous Return. _Phys. Rev. E_ 2020, _101_ , 052130
  * Luo et al. 2022 Luo, Y.; Zeng, C.; Huang, T.; Ai, B.-Q. Anomalous Transport Tuned through Stochastic Resetting in the Rugged Energy Landscape of a Chaotic System with Roughness. _Phys. Rev. E_ 2022, _106_ , 034208
  * Pal et al. 2022 Pal, A.; Kostinski, S.; Reuveni, S. The Inspection Paradox in Stochastic Resetting. _J. Phys. A : Math. Theor._ 2022, _55_ , 021001
  * Thompson et al. 2022 Thompson, A. P.; Aktulga, H. M.; Berger, R.; Bolintineanu, D. S.; Brown, W. M.; Crozier, P. S.; in ’t Veld, P. J.; Kohlmeyer, A.; Moore, S. G.; Nguyen, T. D.; Shan, R.; Stevens, M. J.; Tranchida, J.; Trott, C.; Plimpton, S. J. LAMMPS - A Flexible Simulation Tool for Particle-Based Materials Modeling at the Atomic, Meso, and Continuum Scales. _Comp. Phys. Comm._ 2022, _271_ , 108171
  * Bonomi et al. 2009 Bonomi, M.; Branduardi, D.; Bussi, G.; Camilloni, C.; Provasi, D.; Raiteri, P.; Donadio, D.; Marinelli, F.; Pietrucci, F.; Broglia, R. A.; Parrinello, M. PLUMED: A portable plugin for free-energy calculations with molecular dynamics. _Computer Physics Communications_ 2009, _180_ , 1961–1972
  * Bonomi et al. 2019 Bonomi, M. et al. Promoting transparency and reproducibility in enhanced molecular simulations. _Nature Methods_ 2019, _16_ , 670–673
  * Tribello et al. 2014 Tribello, G. A.; Bonomi, M.; Branduardi, D.; Camilloni, C.; Bussi, G. PLUMED 2: New feathers for an old bird. _Computer Physics Communications_ 2014, _185_ , 604–613
  * Abraham et al. 2015 Abraham, M. J.; Murtola, T.; Schulz, R.; Páll, S.; Smith, J. C.; Hess, B.; Lindahl, E. GROMACS: High performance molecular simulations through multi-level parallelism from laptops to supercomputers. _SoftwareX_ 2015, _1-2_ , 19–25
  * Reuveni 2016 Reuveni, S. Optimal Stochastic Restart Renders Fluctuations in First Passage Times Universal. _Physical Review Letters_ 2016, _116_ , 170601
  * Faradjian and Elber 2004 Faradjian, A. K.; Elber, R. Computing time scales from reaction coordinates by milestoning. _The Journal of Chemical Physics_ 2004, _120_ , 10880–10889
  * Bonati et al. 2021 Bonati, L.; Piccini, G.; Parrinello, M. Deep learning the slow modes for rare events sampling. _Proceedings of the National Academy of Sciences_ 2021, _118_ , e2113533118
  * Mendels et al. 2018 Mendels, D.; Piccini, G.; Parrinello, M. Collective Variables from Local Fluctuations. _The Journal of Physical Chemistry Letters_ 2018, _9_ , 2776–2781
  * Piccini et al. 2018 Piccini, G.; Mendels, D.; Parrinello, M. Metadynamics with Discriminants: A Tool for Understanding Chemistry. _Journal of Chemical Theory and Computation_ 2018, _14_ , 5040–5044
  * Sultan and Pande 2018 Sultan, M. M.; Pande, V. S. Automated design of collective variables using supervised machine learning. _The Journal of Chemical Physics_ 2018, _149_ , 094106
  * Sidky et al. 2020 Sidky, H.; Chen, W.; Ferguson, A. L. Machine learning for collective variable discovery and enhanced sampling in biomolecular simulation. _Molecular Physics_ 2020, _118_ , e1737742
  * Karmakar et al. 2021 Karmakar, T.; Invernizzi, M.; Rizzi, V.; Parrinello, M. Collective variables for the study of crystallisation. _Molecular Physics_ 2021, _119_ , e1893848
  * Massey 1951 Massey, F. J. The Kolmogorov-Smirnov Test for Goodness of Fit. _Journal of the American Statistical Association_ 1951, _46_ , 68–78
  * 54 Miller, L.H. (1956), Table of Percentage Points of Kolmogorov Statistics. Journal of the American Statistical Association. Vol. 51, 111-121. http://www.sciepub.com/reference/142693
  * Ray and Reuveni 2020 Ray, S.; Reuveni, S. Diffusion with resetting in a logarithmic potential. _The Journal of Chemical Physics_ 2020, _152_ , 234110
  * Ray and Reuveni 2021 Ray, S.; Reuveni, S. Resetting transition is governed by an interplay between thermal and potential energy. _The Journal of Chemical Physics_ 2021, _154_ , 171103
  * Evans and Majumdar 2011 Evans, M. R.; Majumdar, S. N. Diffusion with Stochastic Resetting. _Phys. Rev. Lett._ 2011, _106_ , 160601
  * 58 Blumer, O.; Reuveni, S.; Hirshberg, B. Input files for ’resettingMetadynamics’. https://github.com/OfirBlumer/resettingMetadynamics