# Quantum dynamics using path integral coarse-graining

Félix Musil  Department of Physics, Freie Universität Berlin, Arnimallee 12,
14195 Berlin, Germany  Iryna Zaporozhets  Department of Physics, Freie
Universität Berlin, Arnimallee 12, 14195 Berlin, Germany  Department of
Chemistry, Rice University, Houston, Texas 77005, United States  Center for
Theoretical Biological Physics, Rice University, Houston, Texas 77005, United
States  Frank Noé  Department of Physics, Freie Universität Berlin, Arnimallee
12, 14195 Berlin, Germany  Center for Theoretical Biological Physics, Rice
University, Houston, Texas 77005, United States  Department of Chemistry, Rice
University, Houston, Texas 77005, United States  Microsoft Research,
Cambridge, UK  Cecilia Clementi<EMAIL_ADDRESS>Department of
Physics, Freie Universität Berlin, Arnimallee 12, 14195 Berlin, Germany
Center for Theoretical Biological Physics, Rice University, Houston, Texas
77005, United States  Department of Chemistry, Rice University, Houston, Texas
77005, United States  Venkat Kapil<EMAIL_ADDRESS>Yusuf Hamied Department of
Chemistry, University of Cambridge, Lensfield Road, Cambridge, CB2 1EW,UK

###### Abstract

Vibrational spectra of condensed and gas-phase systems containing light nuclei
are influenced by their quantum-mechanical behaviour. The quantum dynamics of
light nuclei can be approximated by the imaginary time path integral (PI)
formulation, but still at a large computational cost that increases sharply
with decreasing temperature. By leveraging advances in machine-learned coarse-
graining, we develop a PI method with the reduced computational cost of a
classical simulation. We also propose a simple temperature elevation scheme to
significantly attenuate the artefacts of standard PI approaches and also
eliminate the unfavourable temperature scaling of the computational cost. We
illustrate the approach, by calculating vibrational spectra using standard
models of water molecules and bulk water, demonstrating significant
computational savings and dramatically improved accuracy compared to more
expensive reference approaches. We believe that our simple, efficient and
accurate method could enable routine calculations of vibrational spectra
including nuclear quantum effects for a wide range of molecular systems.

Predictive simulations of thermodynamic and time-dependent properties of
condensed and gas-phase systems lay the foundations of computational materials
design and discovery 1. Accurate modeling of many systems, such as those
containing light nuclei like H, C, N, and O, must account for their quantum-
mechanical behavior to include the zero-point motion of collective modes 2 and
the tunneling of the system across classically inaccessible barriers 3.
Phenomena emerging from the quantum dynamics of light nuclei are ubiquitous in
chemistry and material science of molecular systems, for instance, the
relative diffusion of $\mathrm{H}_{2}$ in clathrate hydrates 4 and kinetic
isotope effects in porous organic crystals 5, proton-transfer rates in
molecular switches 6, the red-shift in the IR spectra of O–H stretch mode in
ice 7 and the characterization of (bio-)molecular systems using vibrational
spectroscopy 8.

Unfortunately, the exact description of quantum dynamics – requiring the
solution of Schrödinger’s equation – is possible only for the smallest of
systems, such as molecules containing a few atoms 9; 10. Extension to larger
systems requires “local" approximations or truncation of the interaction
potential and/or an approximate solution of the many-body Schrödinger equation
11; 12; 13, akin to the electronic structure problem. Alternatives that render
an “approximate but full" quantum-mechanical treatment of all degrees of
freedom are based on the (semi-)classical dynamics of the system 14. In this
context, imaginary time path integral (PI) simulations, although originally
formulated for incorporating the quantum statistical effects 15, are becoming
increasingly popular for studying the approximate quantum dynamics of
distinguishable particles 14. State-of-the-art PI approaches 16; 17; 18; 19
neglect real time quantum coherence but include effects arising from the
quantum statistical distribution of the nuclei 20. These methods give a good
description of the dynamical response of condensed phase systems, for
instance, vibrational spectra of bulk water 21 and ice 19, where quantum
coherence effects last only over short times. On the other hand, predicting
the vibrational response of molecules and clusters is still a challenge, the
reasons being twofold. Firstly, the artifacts of PI methods (arising from the
neglect of real-time coherence and approximations to non-centroid Matsubara
fluctuations 22), such as spurious broadening, frequency shifts,  and an
incorrect temperature-dependence of the relative intensities of quantal modes
23; 24 become worse with a reduction in temperature 25. Secondly, the
computational cost of PI methods increases steeply with inverse temperature 26
making them expensive to obtain a direct comparison with experiments. In the
last few decades, many new approaches have aimed at individually improving the
accuracy 27; 28; 29; 30 or efficiency 31; 19; 32; 33 of PI dynamical methods.
These studies highlight the urgency for an accurate, efficient, and generally-
applicable method that can treat the quantum dynamics of molecules, clusters,
and bulk systems at the same footing, and aid in modeling their vibrational
response.

Figure 1: Temperature dependent IR spectrum of an O–H bond described by a
Morse oscillator. The IR spectrum of a 2D Morse oscillator mimicking an O–H at
600 K, 300 K, 150 K calculated by (a) solving the Schrödingers equation
numerically (QM) 30, and using (b) centroid molecular dynamics (CMD), and (c)
$\textrm{T}_{\mathrm{e}}$ PIGS using $\textrm{T}_{\mathrm{e}}{}$ = 600 K.
Panel (d) shows the dependence of the $\textrm{T}_{\mathrm{e}}$ PIGS IR
spectrum on the temperature of the centroid free energy surface,
($\textrm{T}_{\mathrm{e}}{}$). Insets show the temperature-dependent radial
probability densities calculated using the respective methods and the dashed
grey vertical line is the mode of quantum-mechanical distribution. Note that
the range of $y$ axes is not kept the same to aid clarity.

In this work, we combine PI and coarse-graining methods via machine learning
to render the calculation of quantum vibrational spectra accurate and
computationally affordable and demonstrate its capabilities on paradigmatic
aqueous systems. Our approach builds upon the state-of-the-art centroid
molecular dynamics 34; 16 (CMD) approach, which time-evolves the system
classically on the free energy surface (FES) of the centroid of the imaginary
time path – a modified PES that includes nuclear quantum effects. Its two key
features, a temperature elevation ($\textrm{T}_{\mathrm{e}}{}$) ansatz and
path integral coarse-graining simulations (PIGS), ensure accuracy and
computational efficiency. On one hand, the $\textrm{T}_{\mathrm{e}}{}$ ansatz
alleviate the spurious redshift in CMD, leading to an improvement in the
accuracy of vibrational spectra over state-of-the-art methods. On the other
hand, with PIGS we machine learn the centroid FES in a general manner, i.e.,
without making prior assumptions about the functional form of the system’s
PES, and use it to evolve the system classically. This method, referred to as
$\textrm{T}_{\mathrm{e}}$ PIGS, enables the calculation of the IR spectra of
aqueous systems including nuclear quantum effects, in excellent agreement with
numerically exact or more expensive reference methods. Furthermore, we
demonstrate that our approach is transferable across phases and temperatures
allowing modeling of vibrational spectra at cryogenic temperatures at orders
of magnitude lower computational cost using classical MD.

We first discuss CMD and introduce $\textrm{T}_{\mathrm{e}}$ PIGS in the
context of a simple yet realistic anharmonic system that highlights the
deficiencies of PI methods at low temperatures: a 2D radial Morse oscillator
mimicking an O–H bond, described by the Hamiltonian,

$\displaystyle\hat{H}=(2\mu)^{-1}\left[\hat{p}_{x}^{2}+\hat{p}_{y}^{2}\right]+D\left[1-e^{-\alpha\left(\sqrt{\hat{q}_{x}^{2}+\hat{q}_{y}^{2}}-r_{0}\right)}\right]^{2},$
(1)

where the parameters $\mu$, $D$, $\alpha$, and $r_{0}$ are defined in section
I.A of the supporting information (SI).  For simplicity, we ignore coupling
between the rotations and vibrations of the system, making this a simple but
physically relevant one-dimensional problem that highlights limitations of PI
based dynamical methods. The O–H bond has a large zero point energy
$E_{0}\approx 1843$ $\text{cm}^{-1}$ $\approx 2652$ K and a $0\to 1$
transition energy $E_{1}-E_{0}\approx 3568$ $\text{cm}^{-1}$ $\approx 5134$ K
(see section I.A of the SI for more details), meaning that even at a
temperature of 600 K the system resides almost exclusively in its ground
state, yet probes anharmonic regions of the PES due to zero-point motion. As
shown in Fig. 1(a), this results in a temperature-independent line position of
the IR spectrum, at least up to 600 K, but also a red-shift of around 200
$\text{cm}^{-1}$ with respect to the classical spectrum due to quantum nuclear
motion.

Figure 2: The curvature problem in centroid molecular dynamics Contour plots
of a 2D Morse radial potential up to half the dissociation energy (black to
grey solid lines), the probability density of the system (red to white
indicating high to low), and a snapshot of the ring polymer (red) and its
centroid (blue) at (a) 150 K and (b) 600 K.

The CMD approach 16 is based on the imaginary time path integral isomorphism
35 between the thermodynamics of a quantum system at inverse temperature
$\beta$ and a classical ring polymer made of $P$ replicas of the system at
$\beta_{P}=\beta/P$,

$\displaystyle
Z=\text{Tr}[e^{-\beta\hat{H}}]\propto\lim_{P\to\infty}\int\textrm{d}{}\mathbf{q}~{}e^{-\beta_{P}\left[\sum_{j}U(\mathbf{q}^{(j)})+U^{\textrm{spr}}(\mathbf{q})\right]},$
(2)

where $\mathbf{q}\equiv\\{\mathbf{q}^{(1)},\dots,\mathbf{q}^{(P)}\\}$ is a
shorthand for positions of the $P$ replicas of the system with
$\mathbf{q}^{(j+P)}\equiv\mathbf{q}^{(j)}$ implied, $U(\mathbf{q})$ defines
the classical PES, and $U^{\textrm{spr}}(\mathbf{q})$ is a temperature
dependent spring term 35 that connects consecutive replicas of the system.
Within CMD, the system is time evolved classically on
$U^{\mathrm{CMD}}(\mathbf{q}_{\textrm{c}};\beta)$, defined as the free energy
surface of the centroid of the imaginary time path (modulo a constant) at
$\beta$,

$\displaystyle
U^{\mathrm{CMD}}(\mathbf{q}_{\textrm{c}};\beta)=-\beta^{-1}\log\left\langle~{}\delta\Big{(}\frac{1}{P}\sum_{j=1}^{P}\mathbf{q}^{(j)}-\mathbf{q}_{c}\Big{)}\right\rangle_{\beta},$
(3)

where $\sum_{j=1}^{P}\mathbf{q}^{(j)}/P$ is the centroid of the ring polymer
and $\langle\cdot\rangle$ an average over the path integral Hamiltonian in Eq.
2. The thermodynamic force acting on the centroid can be calculated on the fly
from a constrained PI simulation at each CMD step 16

$\displaystyle\mathbf{f}^{\mathrm{CMD}}(\mathbf{q}_{\textrm{c}};\beta)=-\nabla
U^{\mathrm{CMD}}(\mathbf{q}_{\textrm{c}};\beta)$
$\displaystyle=\left\langle\mathbf{f}_{\textrm{c}}~{}\delta\Big{(}\frac{1}{P}\sum_{j=1}^{P}\mathbf{q}^{(j)}-\mathbf{q}_{\textrm{c}}\Big{)}\right\rangle_{\beta},$
(4)

where, $\mathbf{f}_{\textrm{c}}=\sum_{j=1}^{P}\mathbf{f}^{(j)}/P$, and
$\mathbf{f}^{(j)}$ is the physical force acting on the $j$-th replica.
Alternatively, the centroid can be evolved on its FES within a PI simulation
in a partially adiabatic manner 36 by decoupling the centroid from the rest of
the system. Dynamical properties, such as the IR spectrum, can be easily
calculated via classical time correlation functions (TCFs) on the basis of the
centroid trajectory, akin to classical MD.

As shown in Fig. 1(b), the IR spectrum of the O–H bond computed with CMD is in
excellent agreement with the numerically exact result at 600 K. Unfortunately,
at lower temperatures, the system experiences a spurious red shift that gets
worse as the temperature is reduced. This well-known artifact is referred to
as the “curvature problem" 25 and arises in cases where the ring polymer has a
shape such that its centroid lies outside the ring (see Fig. 2) and therefore
does not represent the quantum-mechanical probability density of the system
(as obtained by the replicas). Interestingly, the curvature problem is purely
a structural artifact, and, as shown in the inset of Fig. 1(b), it can be
diagnosed from PI trajectories: the misalignment of the distribution of the
centroid at 300 K and 150 K to the (physical) probability density obtained
from the replicas. As mentioned above, another issue with CMD, and more
generally with PI methods, is that the required number of replicas scales
inversely with temperature and the maximum physical frequency of the system,
$P\sim\beta\hbar\omega_{\mathrm{max}}$ 37. Therefore its computational cost
increases steeply as the temperature is reduced. These challenges prevent
investigations of systems at cryogenic temperatures where most experimental
data is available.

Recently, Trenins  _et al._ 30 have proposed a quasi centroid molecular
dynamics (QCMD) scheme, that evolves the system on the FES of an ad hoc
curvilinear function of replica positions – a quasi centroid that does not
“fall out" of the hull of the path integral. A careful selection of this
function doesn’t hamper quantum Boltzmann statistics and results in a compact
PI ring polymer that alleviates the curvature problem. This results in an
excellent agreement of the vibrational spectrum of a molecule of water, liquid
water, and proton disordered hexagonal ice 38 with (numerically exact)
reference methods. More recently, Fletcher  _et al._ 33 have proposed an
efficient fast QCMD scheme that avoids the costly on the fly quasi-centroid
forces by precomputing an analytic FES on the basis of a PI trajectory. This
reduces the cost of predicting an accurate vibrational spectrum to that of
classical MD, resulting in an improvement over standard PI methods for the
calculation of vibrational spectra of exemplary molecular systems. 33
Unfortunately, the computational cost of these approaches still grows
unfavorably with temperature due to the need for a low-temperature PI
trajectory for fitting the FES. Moreover, an extension to general systems
requires a general/universal procedure to fit the FES of the centroid and
careful knowledge of appropriate curvilinear coordinates that do not suppress
sampling of the physical regions of the configurational space 39. Nonetheless,
the advancements brought forth by (fast) QCMD 30; 33 suggest that further
improvements in CMD provide a promising route towards the accurate and
efficient calculation of vibrational spectra.

Taking inspiration from fast implementations of CMD 40; 41 we propose path
integral coarse-grained simulations (PIGS) that perform MD on a modified PES –
including quantum nuclear effects. In particular, we leverage the recent
developments in the definition of coarse-grained machine learning potentials
42; 43; 44 and of high-order correlation functions 45; 46 to obtain this
modified potential by coarse-graining the imaginary time path integral to that
of a classical system. In this work, we use PIGS to accelerate CMD in a
general manner, i.e. it is applicable to systems exhibiting a wide range of
interparticle interactions. To obtain the centroid FES at
$\beta_{\mathrm{e}}$, we use the force matching method 47; 48 that is
typically used to build thermodynamically consistent bottom-up coarse-grained
models for macromolecular systems, i.e. the coarse-grained model reproduces
the thermodynamic properties of the all-atom system projected onto the coarse-
grained coordinates. It has been shown 48 that Eq. 3 can be recast as a
variational principle: $U^{\mathrm{CMD}}(\mathbf{q}_{c};\beta)$ corresponds to
the minimum of the force matching functional:

$U^{\mathrm{CMD}}(\mathbf{q}_{c};\beta)=\operatorname*{arg\,min}_{G\in
C(\mathbb{R}^{n})}\left\langle\Big{\|}\mathbf{f}_{c}+\nabla
G(\mathbf{q}_{c})\Big{\|}^{2}\right\rangle_{\beta},$ (5)

where the average is performed at inverse temperature $\beta$ with the PI
centroid constrained at $\mathbf{q}_{\textrm{c}}\in\mathbb{R}^{n}$, and the
minimization is performed over the space $C(\mathbb{R}^{n})$ of all real
continuous functions $G:\mathbb{R}^{n}\rightarrow\mathbb{R}$. In this work we
optimize a machine learning model as a surrogate for the potential of mean
force to reproduce the thermodynamics of the centroid without having to
explicitly simulate $P$ replicas of the system. In practice, we only learn the
difference between the centroid FES and the classical PES, as it is done in
Ref. 40 This keeps the amount of training data to a minimum and provides an
appropriate prior in the absence of data or for collective modes that do not
exhibit quantum nuclear effects.

The centroid potential of mean force associated with an atomic configuration
is expressed as a sum of the classical PES associated with $U$, and atom-
centered contributions

$\tilde{U}^{\mathrm{CMD}}(\mathbf{q}_{c};\bm{\theta})=U(\mathbf{q}_{c})+\sum_{i}A_{\textrm{a}_{i}}(\mathbf{q}_{c}-\mathbf{q}_{c~{}i};\bm{\theta}),$
(6)

where $\bm{\theta}$ is a set of model parameters,
$\mathbf{q}_{c}\equiv\\{\mathbf{q}_{c~{}1},\dots,\mathbf{q}_{c~{}N}\\}$ is a
shorthand for set of atomic positions of a structure, $\mathbf{q}_{c~{}i}$ is
the position of the $i^{\textrm{th}}$ atom of species $\textrm{a}_{i}$, and
the $A_{\textrm{a}_{i}}$ functions are parameterized using a machine learning
model that captures the multi-body interactions emerging from the coarse-
graining procedure 49. We model the atomic potentials of mean force,
$A_{\textrm{a}_{i}}$, by representing the atomic environment with the
normalized SOAP powerspectrum 50; 51; 52 and pass these 3-body features to a
multi-layer perceptron with $3$ hidden layers of width $[400,200,200]$ and the
hyperbolic tangent activation function. The parameters of the model are
obtained by minimizing the force matching loss

$\left\|\mathbf{f}_{c}+\nabla\tilde{U}^{\mathrm{CMD}}(\mathbf{q}_{c};\bm{\theta})\right\|^{2},$
(7)

over a set of centroid sample forces $\mathbf{f}_{c}$ and configurations
$\mathbf{q}_{c}$ obtained from a PIMD simulation at $\beta$ performed with the
i-PI code 53 (see section I.C of the SI for more details). The model and its
training have been implemented in pytorch 54 and the codes are available upon
request.

Figure 3: Quantum dynamics of a water molecule. Comparison of the vibrational
spectrum of a water molecule described by the Partridge-Schwenke model55 at
(a) 600 K, (b) 300 K, (c) 150 K, using discrete value representation (yielding
numerically exact results) in black, the proposed $\textrm{T}_{\mathrm{e}}{}$
PIGS approach using the centroid free energy surface calculated at 600 K in
red, CMD in blue, and TRPMD in green. The insets show the probability
distribution of the O–H bond length calculated from the CMD and
$\textrm{T}_{\mathrm{e}}$ PIGS and the black dashed line indicates the mode of
the distribution obtained from PIMD – constant across temperatures. Panel (d)
displays the temperature-dependent IR spectrum of a water molecule from 600 K
down to 0.5 K. The gray dashed line indicates the reference $0\rightarrow 1$
transition frequency while the black line is shifted by 40 $\text{cm}^{-1}$ to
reflect the results obtained with Matsubara dynamics 22.

We propose a simple and physically motivated temperature elevation
($\textrm{T}_{\mathrm{e}}{}$) ansatz that alleviates the curvature problem of
CMD and also eliminates the unfavorable temperature dependence associated with
the computational cost of vibrational spectra. We note that for a system in
its ground state, i.e. $\hbar\omega>>\beta^{-1}$, the IR spectrum – related by
the dipole correlation function – is only trivially dependent on the
temperature (see section II.A of the SI). It is easy to show that after
rescaling with the inverse temperature $\beta$, the Kubo-transformed time-
correlation function is a constant if the system is in its ground state, i.e.,
$\beta\tilde{C}_{\hat{\bm{\mu}}\hat{\bm{\mu}}}(\omega;\beta)=\mathrm{constant}$.
As shown in Fig. 1(b), CMD does not follow this behavior because of the
curvature problem at low temperatures.  To avoid this artifact, we propose a
$\textrm{T}_{\mathrm{e}}{}$ ansatz, i.e. we rewrite a time-correlation
function in terms of a CMD time-correlation function computed at an elevated
temperature $\textrm{T}_{\mathrm{e}}{}$,

$\displaystyle\tilde{C}^{\textrm{T}_{\mathrm{e}}{}}_{\hat{\bm{\mu}}\hat{\bm{\mu}}}(\omega,\beta;\beta_{\mathrm{e}})=\frac{\beta_{\mathrm{e}}}{\beta}\tilde{C}^{~{}\mathrm{CMD}}_{\hat{\bm{\mu}}\hat{\bm{\mu}}}(\omega,\beta_{\mathrm{e}}).$
(8)

This ansatz is valid for any quantum system (of distinguishable particles) and
can be used to improve PI based methods that exhibit diminishing performance
at low temperatures. In order to ease the calculation of quantum dynamical
properties for general (high-dimensional) systems, we approximate this
relation as

$\displaystyle\tilde{C}^{\textrm{T}_{\mathrm{e}}{}}_{\hat{\bm{\mu}}\hat{\bm{\mu}}}(\omega,\beta;\beta_{\mathrm{e}})\approx
Z^{-1}\int\textrm{d}{}\mathbf{q}^{\prime}~{}e^{-\beta
U^{\mathrm{CMD}}(\mathbf{q}^{\prime};\beta_{e})}~{}\bm{\mu}(0)\cdot\bm{\mu}(t),$
(9)

so that it is simply estimated by evolving the centroid at $\beta$ on a FES
calculated at a high temperature $\beta_{\mathrm{e}}$. This approximate form
of Eq.8 bears advantages such as it is exact in the harmonic and the classical
limit, and doesn’t require any posterior rescaling of the TCF. These limits
also suggest that Eq.9 should give a good description of stiff modes, weakly
coupled with the rest of the system, without perturbing the dynamics of the
low frequency (classical) modes.  The use of an elevated temperature also
improves sampling efficiency, akin to the use of a high temperature in the
adiabatic free energy dynamics to improve sampling of a high dimensional
configurational space. The elevated temperature also permits the use of a
small number of replicas to discretize the imaginary time path. In this work,
we perform fully adiabatic CMD by separately computing the FES at a higher
temperature, however, we also plan to implement and study a partially
adiabatic implementation of the $\textrm{T}_{\mathrm{e}}{}$ ansatz with CMD.

As shown in Fig. 1(b), the CMD IR spectrum remains in excellent agreement with
the exact result at 600 K and the radial distribution function of the centroid
of the O–H bond is aligned with the physical distribution. These observations
suggest that the system doesn’t exhibit the curvature problem at 600 K and
thus we test the $\textrm{T}_{\mathrm{e}}{}$ anstaz for $T_{\mathrm{e}}=600$ K
where $\beta_{\mathrm{e}}=(k_{\mathrm{B}}{}T_{\mathrm{e}})^{-1}$. As shown in
Fig. 1(c) and Fig. 1(d), increasing the “elevated temperature" progressively
improves the description of the IR spectrum at 150 K. Moreover, using
$T_{\mathrm{e}}=600$ K leads to an excellent agreement of the temperature
dependent IR spectrum with the exact result. Note that for the Morse potential
there exists a wide window of suitable $T_{\mathrm{e}}$ (see section II.B of
the SI for more details), i.e. large enough to alleviate the curvature problem
but also small enough for the harmonic approximation in Eq. 9 to be valid, and
we observed similar features with the water molecule and bulk water. Fig. 1(c)
shows that the line position of the predicted IR spectra and the radial
distribution of the O–H bond are largely temperature independent using
$T_{\mathrm{e}}=600$ K, as expected for a system in its ground state. Finally,
the number of replicas needed to calculate
$U^{\mathrm{CMD}}(\mathbf{q}^{\prime};\beta_{e})$ is
$P\approx\beta_{e}\hbar\omega_{\mathrm{max}}<\beta\hbar\omega_{\mathrm{max}}$
which is independent of $\beta$. Thus, within the $T_{e}$ ansatz, the cost of
simulating the quantum dynamics of a system in its ground state at $\beta$
doesn’t scale with temperature.

The workflow for computing the quantum vibrational spectrum of a system using
PIGS and the $\textrm{T}_{\mathrm{e}}{}$ ansatz can be summarized as follow.
First, we perform short PIMD simulations exploring a range of temperatures and
select the lowest temperature $T_{\mathrm{e}}$ by checking for the alignment
between the centroid and the physical radial distributions.  In this study, we
use radial distribution functions (RDFs) to check for this alignment, however,
more general systems might require using many-body correlation functions such
as SOAP 50 or ACE 56 (of which RDFs form a subset), combined with state-of-
the-art dimensionality reduction schemes to compare centroid and physical
distributions. Then, we use Eqs. 7 and 6 to machine-learn the centroid FES at
$\beta_{\mathrm{e}}$ as a sum of local atom-centered components. In the final
step, we predict quantum dynamical properties by performing MD at the desired
temperatures. These simple steps enable the development of an effective PES
that includes quantum nuclear effects for dynamics – an FES trained on a
single high-temperature PI trajectory – and that is transferable across
temperature. The computational details of all the simulations are described in
sections I.B, I.D, and I.E of the SI.

Figure 4: Quantum dynamics of bulk water. Comparison of the vibrational
spectrum of (a) liquid water at 300 K and (b) hexagonal ice at 150 K,
described by a q-TIP4P/f 57 water model, calculated using QCMD in black, the
proposed $\textrm{T}_{\mathrm{e}}$ PIGS approach using the centroid free
energy surface calculated at 600 K in red, CMD in blue, and TRPMD in green.
Panel (c) displays the $\textrm{T}_{\mathrm{e}}$ PIGS IR spectrum predictions
for hexagonal ice from 150 K down to 0.5 K.

We demonstrate the capabilities of the $\textrm{T}_{\mathrm{e}}$ PIGS approach
by applying it to the IR spectrum of a water molecule – a challenging system
that exhibits a strong red-shift of the stretching modes due to zero-point
motion and a fine temperature-dependent splitting of the vibrational peaks due
to the coupling of rotations and vibrations. We study the IR spectrum at 150
K, 300 K, and 600 K using the well-known Partridge-Schwenke model 55, which
exhibits spectroscopic accuracy, and compare with experimental and numerically
exact results 58. We also compare with state-of-the-art approaches like CMD
and thermostatted ring polymer molecular dynamics (TRPMD) 18.

The water molecule resides in its vibrational ground state at least up to 600
K, but exhibits a sensitive dependence of the ro-vibrational splitting of the
modes with temperature, as seen from the numerically exact IR spectra 30 in
Fig. 3. The TRPMD approach largely captures the correct line positions of the
(envelopes of) the stretching and the bending bands but is artificially
broadened 18. The CMD approach gives a good description of the full IR
spectrum at 600 K but is artificially red-shifted at 300 K and 150 K due to
the onset of the curvature problem. This is confirmed by observing in Fig. 3
that the centroid distributions at 300 K and 150 K indicate a lower
(unphysical) O–H bond length and misalignment with the physical radial
distribution. Nevertheless, at 600 K we do not observe any “symptoms" of the
curvature problem, and thus we use $T_{\mathrm{e}}=600$ K for the
$\textrm{T}_{\mathrm{e}}$ PIGS approach.

As shown in Fig. 3, $\textrm{T}_{\mathrm{e}}$ PIGS spectra at 600 K are in
excellent agreement with the exact and the independent CMD results, underlying
the accuracy of the learned FES. More importantly, $\textrm{T}_{\mathrm{e}}$
PIGS describes the fundamental frequencies and the rotational splittings of
the stretching and the bending modes at 300 K and 150 K, a substantial
improvement over state-of-the-art PI methods in terms of accuracy.  We note
however that, like CMD, $\textrm{T}_{\mathrm{e}}$ PIGS does not capture the
subtle temperature-dependence of the relative intensities of the bending mode
w.r.t to the stretching mode. This is largely an artefact of PI-based
approximate quantum dynamics methods which lack quantum coherence in their
dynamics 23 arising from the momenta not being drawn from their quantum
Boltzmann distribution 24. On the other hand, $\textrm{T}_{\mathrm{e}}$ PIGS
is computationally efficient since the FES is estimated from a 100 ps PI
trajectory with 8 replicas and a timestep of 0.5 fs, suggesting computational
gains of at least a factor of $500\times$, $1000\times$ and $2000\times$ at
600 K, 300 K and 150 K, respectively, if one knows an appropriate
$\textrm{T}_{\mathrm{e}}{}$, compared to a 1 ns long CMD simulations needed to
converge the TCFs.  The overall cost for identifying
$\textrm{T}_{\mathrm{e}}{}$, namely a set of 100 ps PIMD simulations across
100 K-600 K, was at least two orders of magnitude lower than that of a CMD
simulation needed to calculate the spectrum at 150 K. Moreover, we expect this
cost to be overestimated, due to the over-conservative nature of our
benchmark. These could easily be made an order of magnitude less-expensive
using out-of-the-box accelerated-sampling methods 15 based on generalized
Langevin equation thermostats 59; 60; 61 and high-order splittings of the
Boltzmann operator 62; 63; 64; 32 combined with replica-exchange 65 across 100
K-600 K. Other accelerated techniques like ring polymer contraction and
multiple time stepping-37; 66 promise a classical computational-cost for the
calculation of quantum dynamical properties 21; 31, contingent to finding an
inexpensive surrogate for the high-frequency modes of the system, however,
they bear the inaccuracies of traditional path-integral methods that we
address using the $\textrm{T}_{\mathrm{e}}{}$ ansatz.

To showcase the absence of computational scaling with temperature, we
calculate the IR spectrum of a water molecule down to 0.5 K. Decrease in
temperature expectedly reduces the rotational splitting of the vibrational
modes, and below 10 K, the system falls into its rotational ground state with
well resolved peaks for the three vibrational modes. This temperature is in
excellent agreement with the experimentally known rotational constants of a
water molecule, i.e. 13-35 K 67. Furthermore, the frequencies of the ground
state vibrational bending mode matches the experimental / exact results 58 up
to 8 $\text{cm}^{-1}$ – the resolution of calculated spectra – and those of
the two stretching modes are expectedly blue-shifted by around 40
$\text{cm}^{-1}$ in excellent agreement with the results obtained from
Mastubara dynamics 22 – the true reference for PI based dynamical methods. We
emphasize that the calculation of quantum mechanical spectrum at cryogenic
temperatures by means of classical dynamics at this level of accuracy-to-cost
ratio is unprecedented to our knowledge.

As a final test of our approach, we study the IR spectrum of condensed phase
aqueous systems: bulk water at 300 K and hexagonal ice at 150 K. The presence
of inter-molecule or crystal modes that couple with high-frequency modes
should constitute a challenge for $\textrm{T}_{\mathrm{e}}$ PIGS as it has
been the case with previous PI approaches 28; 27. An additional challenge is
that ice melts below the onset temperature of the curvature problem which
could complicate the calculation of the FES at $\beta_{\mathrm{e}}$. To
circumvent this issue, we exploit the local nature of the FES (see Eq. 6) and
insights from deep inelastic neutron scattering experiments 7 that probe the
quantum nuclear motion of atoms. These experiments suggest that the local
potential felt by the nuclei due to quantum delocalization is short ranged 68
and sensitive/unique to their local environments 69. Given that Eq. 6 makes
the fitted FES size-extensive and local environments of H and O atoms in ice
are present in liquid water 70, we conjecture that the
$\textrm{T}_{\mathrm{e}}{}$ FES of hexagonal ice can be constructed from a
high-temperature simulation of liquid water. We estimate vibrational spectra
using the q-TIP4P/f water potential 57 and a linear dipole moment surface.
Although this model doesn’t exhibit “experimental" or “spectroscopic"
accuracy, it has been extensively used for comparing the performance of
various PI 71; 19; 38; 30 and wavefunction-based11; 12 methods and exhibits
good agreement with the experimental spectra when combined with an appropriate
dipole moment surface 72.

As shown in Fig. 4, we obtain well resolved spectra for both liquid water at
300 K and hexagonal ice at 150 K using a $T_{\mathrm{e}}=600$ K
$k_{\mathrm{B}}{}$ FES fitted on a 10 ps PI simulation of liquid water. Our
results are in good agreement with CMD, TRPMD and QCMD for room-temperature
liquid water, where the curvature problem is small 73. In the case of
hexagonal ice at 150 K, we see a quantitative and qualitative improvement on
the description of the high frequency band with respect to CMD and TRPMD and
an excellent agreement with QCMD 30. In addition to the increased accuracy,
our approach is over three orders of magnitude less expensive (as measured by
the number of force evaluations) than state-of-the-art approaches at 300 K and
150 K. Furthermore, the absence of temperature-dependent scaling behaviour of
the cost, allows us to also compute the IR spectrum of hexagonal ice all the
way down to 0.5 K. As expected, we observe that the system resides in its
ground state and does not exhibit any line shifts of the high frequency modes.

In summary, we propose a new approach, $\textrm{T}_{\mathrm{e}}$ PIGS that
accurately simulates the quantum dynamics of light nuclei at the cost of
classical MD by combining PI quantum mechanics, machine learning and bottom up
coarse-graining. This combination is both timely and significant as
exemplified by other recent work(that we discovered out during the review
process) that implement a fast version of CMD 74; 75 using machine learning
potentials. We use a physically-motivated temperature elevation ansatz that
moves the system on the centroid FES obtained at a high temperature, and
fitted efficiently by using machine-learning-based coarse-graining approaches
42; 43; 44. The $\textrm{T}_{\mathrm{e}}{}$ ansatz is exact in the high-
frequency harmonic limit where the vibrational mode essentially lives in its
ground state, as well as for low-frequency modes for which the ring polymer
distribution collapses on the centroid. Furthermore, our study of condensed
aqueous phases water shows that it also performs well for intermediate
frequencies $\hbar\omega\sim\beta^{-1}$ , for instance, the librational /
rotational modes. These limits suggest that our approach could be useful for
studying a wide range of systems with high-frequency modes that exhibiting
weak coupling to rest of the system. We show that the
$\textrm{T}_{\mathrm{e}}$ PIGS FES is _transferable_ across temperatures and
phases, but is also size extensive, meaning that it allows for further
computational savings by learning the FES on a smaller system. We believe that
our approach constitutes a substantial improvement in accuracy over routinely
used state-of-the-art methods like CMD and TRPMD. In addition, its low
computational cost and the absence of its scaling behavior with temperature
allows for accurate the IR spectra predictions even at cryogenic temperatures
– considered prohibitive with state-of-the-art methods.

The routine use of the $\textrm{T}_{\mathrm{e}}{}$ PIGS approach on materials
and chemical system will require a careful and thorough study on a diverse set
of “difficult" systems which could challenge the accuracy of the
$\textrm{T}_{\mathrm{e}}{}$ ansatz. For instance, one could imagine highly
fluxional molecules, system exhibiting a near-continuum of strongly coupled
high-frequency modes, or dynamical processes dominated by “rare" quantum
tunnelling events, being testing cases probing the regimes in which the
$\textrm{T}_{\mathrm{e}}{}$ ansatz is not exact or expected to be accurate.
Similarly, the extension of $\textrm{T}_{\mathrm{e}}{}$ PIGS to general
systems described by first-principles electronic structure methods, will
require development of a hierarchical framework leveraging recent advances in
active-learning strategies for the development of accurate and reliable
machine-learning PES, enhanced sampling of the quantum Boltzmann distribution
using accelerated path-integral methods, and subsequently the calculation of
the centroid FES at an elevated temperature. Additionally, these systems would
also form an ideal test bed for understanding what range of elevated
temperatures could be used to treat a frozen ground-state vibrational mode.
Could a “universal" range of $\textrm{T}_{\mathrm{e}}{}\in[500,600]$ K. be
used for chemical systems, displaying high-frequency modes in the 3000 - 5000
cm${}^{\text{-1}}$ regime, to completely avoid process of selecting an
elevated temperature for a general system? And finally, one could envisage the
use of the PIGS approach for calculations of equilibrium properties. We plan
to address all these directions in future studies.  Overall, we believe that
the simplicity, low cost and high accuracy of $\textrm{T}_{\mathrm{e}}{}$ PIGS
could open up prospects for routine modeling of quantum-vibrational spectra of
general systems and direct comparisons with experiments, often performed at
low temperatures 76.

## Acknowledgements

We thank Start Althorpe, George Trenins, Christoph Schran, Angelos
Michaelides, and members of the Clementi’s group at FU for insightful
discussions and comments on the manuscript. We also thank George Trenins for
sharing QCMD results. V.K. acknowledges funding from the Swiss National
Science Foundation (SNSF) under Project $\text{P2ELP2}\\_\text{191678}$ and
the Ernest Oppenheimer Fund, allocation of CPU hours by CSCS under Project ID
s1000, and support from Churchill College, University of Cambridge. C.C.
acknowledges funding from the Deutsche Forschungsgemeinschaft DFG (SFB/TRR
186, Project A12; SFB 1114, Projects B03 and A04; SFB 1078, Project C7; and
RTG 2433, Project Q05), the National Science Foundation (CHE-1738990,
CHE-1900374, and PHY-2019745), and the Einstein Foundation Berlin (Project
0420815101). F.N. acknowledges funding from the Deutsche
Forschungsgemeinschaft DFG (SFB 1114, Projects A04 and B08), The Berlin
Mathematics center MATH+ (Projects AA1-6, AA2-8), The Berlin Institute for the
Foundations of Learning Data (BIFOLD) and the European Commission (ERC CoG
772230). F.M. acknowledges support from the SNSF under the Postdoc.Mobility
fellowship P500PT_203124 and from the Physics department at FU Berlin for the
computational time.

## Source data

All the raw data and code required to reproduce the results, and scripts and
raw data needed to make the figures is made available in the online
repository: venkatkapil24/pigs-methodology.

## References

  * Houk and Liu (2017) K. N. Houk and F. Liu, Accounts of Chemical Research 50, 539 (2017).
  * Pereyaslavets  _et al._ (2018) L. Pereyaslavets, I. Kurnikov, G. Kamath, O. Butin, A. Illarionov, I. Leontyev, M. Olevanov, M. Levitt, R. D. Kornberg, and B. Fain, Proceedings of the National Academy of Sciences 115, 8878 (2018).
  * Benoit  _et al._ (1998) M. Benoit, D. Marx, and M. Parrinello, Nature 392, 258 (1998).
  * Cendagorta  _et al._ (2021) J. R. Cendagorta, H. Shen, Z. Bačić, and M. E. Tuckerman, Advanced Theory and Simulations 4, 2000258 (2021).
  * Liu  _et al._ (2019) M. Liu, L. Zhang, M. A. Little, V. Kapil, M. Ceriotti, S. Yang, L. Ding, D. L. Holden, R. Balderas-Xicohténcatl, D. He, R. Clowes, S. Y. Chong, G. Schütz, L. Chen, M. Hirscher, and A. I. Cooper, Science 366, 613 (2019).
  * Litman  _et al._ (2019) Y. Litman, J. O. Richardson, T. Kumagai, and M. Rossi, Journal of the American Chemical Society 141, 2526 (2019).
  * Burnham  _et al._ (2006) C. J. Burnham, G. F. Reiter, J. Mayers, T. Abdul-Redah, H. Reichert, and H. Dosch, Physical Chemistry Chemical Physics 8, 3966 (2006).
  * Rossi  _et al._ (2015) M. Rossi, W. Fang, and A. Michaelides, The Journal of Physical Chemistry Letters 6, 4233 (2015).
  * Li and Guo (2001) G. Li and H. Guo, Journal of Molecular Spectroscopy 210, 90 (2001).
  * Larsson  _et al._ (2022) H. R. Larsson, M. Schröder, R. Beckmann, F. Brieuc, C. Schran, D. Marx, and O. Vendrell, “State-resolved infrared spectrum of the protonated water dimer: Revisiting the characteristic proton transfer doublet peak,”  (2022), arXiv:2206.12029 [physics, physics:quant-ph].
  * Wang and Bowman (2011) Y. Wang and J. M. Bowman, The Journal of Chemical Physics 134, 154510 (2011).
  * Gruenbaum  _et al._ (2013) S. M. Gruenbaum, C. J. Tainter, L. Shi, Y. Ni, and J. L. Skinner, Journal of Chemical Theory and Computation 9, 3109 (2013).
  * Qu  _et al._ (2021) C. Qu, P. L. Houston, R. Conte, A. Nandi, and J. M. Bowman, The Journal of Physical Chemistry A 125, 5346 (2021).
  * Althorpe (2021) S. C. Althorpe, The European Physical Journal B 94, 155 (2021).
  * Markland and Ceriotti (2018) T. E. Markland and M. Ceriotti, Nature Reviews Chemistry 2, 1 (2018).
  * Cao and Voth (1994) J. Cao and G. A. Voth, The Journal of Chemical Physics 100, 5106 (1994).
  * Craig and Manolopoulos (2004) I. R. Craig and D. E. Manolopoulos, The Journal of Chemical Physics 121, 3368 (2004).
  * Rossi  _et al._ (2014a) M. Rossi, M. Ceriotti, and D. E. Manolopoulos, The Journal of Chemical Physics 140, 234116 (2014a).
  * Kapil  _et al._ (2020) V. Kapil, D. M. Wilkins, J. Lan, and M. Ceriotti, The Journal of Chemical Physics 152, 124104 (2020).
  * Hele  _et al._ (2015) T. J. H. Hele, M. J. Willatt, A. Muolo, and S. C. Althorpe, The Journal of Chemical Physics 142, 134103 (2015).
  * Marsalek and Markland (2017) O. Marsalek and T. E. Markland, The Journal of Physical Chemistry Letters 8, 1545 (2017).
  * Trenins and Althorpe (2018) G. Trenins and S. C. Althorpe, The Journal of Chemical Physics 149, 014102 (2018), publisher: American Institute of Physics.
  * Benson and Althorpe (2021) R. L. Benson and S. C. Althorpe, The Journal of Chemical Physics 155, 104107 (2021).
  * Plé  _et al._ (2021) T. Plé, S. Huppert, F. Finocchi, P. Depondt, and S. Bonella, The Journal of Chemical Physics 155, 104108 (2021).
  * Witt  _et al._ (2009) A. Witt, S. D. Ivanov, M. Shiga, H. Forbert, and D. Marx, The Journal of Chemical Physics 130, 194510 (2009).
  * Uhl  _et al._ (2016) F. Uhl, D. Marx, and M. Ceriotti, The Journal of Chemical Physics 145, 054101 (2016).
  * Rossi  _et al._ (2017) M. Rossi, V. Kapil, and M. Ceriotti, The Journal of Chemical Physics 148, 102301 (2017).
  * Willatt  _et al._ (2018) M. J. Willatt, M. Ceriotti, and S. C. Althorpe, The Journal of Chemical Physics 148, 102336 (2018).
  * Liu (2014) J. Liu, The Journal of Chemical Physics 140, 224107 (2014).
  * Trenins  _et al._ (2019) G. Trenins, M. J. Willatt, and S. C. Althorpe, The Journal of Chemical Physics 151, 054109 (2019).
  * Kapil  _et al._ (2016a) V. Kapil, J. VandeVondele, and M. Ceriotti, The Journal of Chemical Physics 144, 054111 (2016a).
  * Shepherd  _et al._ (2021) S. Shepherd, J. Lan, D. M. Wilkins, and V. Kapil, The Journal of Physical Chemistry Letters , 9108 (2021).
  * Fletcher  _et al._ (2021) T. Fletcher, A. Zhu, J. E. Lawrence, and D. E. Manolopoulos, The Journal of Chemical Physics 155, 231101 (2021).
  * Cao and Voth (1993) J. Cao and G. A. Voth, The Journal of Chemical Physics 99, 10070 (1993).
  * Chandler and Wolynes (1981) D. Chandler and P. G. Wolynes, The Journal of Chemical Physics 74, 4078 (1981).
  * Hone  _et al._ (2006) T. D. Hone, P. J. Rossky, and G. A. Voth, The Journal of Chemical Physics 124, 154103 (2006).
  * Markland and Manolopoulos (2008) T. E. Markland and D. E. Manolopoulos, The Journal of Chemical Physics 129, 024105 (2008).
  * Benson  _et al._ (2019) R. L. Benson, G. Trenins, and S. C. Althorpe, Faraday Discussions 221, 350 (2019).
  * Haggard  _et al._ (2021) C. Haggard, V. G. Sadhasivam, G. Trenins, and S. C. Althorpe, The Journal of Chemical Physics 155, 174120 (2021).
  * Hone  _et al._ (2005) T. D. Hone, S. Izvekov, and G. A. Voth, Journal of Chemical Physics 122, 054105 (2005).
  * Paesani  _et al._ (2006) F. Paesani, W. Zhang, D. A. Case, T. E. Cheatham, and G. A. Voth, The Journal of Chemical Physics 125, 184507 (2006).
  * Wang  _et al._ (2019) J. Wang, S. Olsson, C. Wehmeyer, A. Pérez, N. E. Charron, G. de Fabritiis, F. Noé, and C. Clementi, ACS Central Science 5, 755 (2019).
  * Husic  _et al._ (2020) B. E. Husic, N. E. Charron, D. Lemm, J. Wang, A. Pérez, M. Majewski, A. Krämer, Y. Chen, S. Olsson, G. de Fabritiis, F. Noé, and C. Clementi, J. Chem. Phys. 153, 194101 (2020).
  * Wang  _et al._ (2020) J. Wang, S. Chmiela, K.-R. Müller, F. Noé, and C. Clementi, J. Chem. Phys. 152, 194106 (2020).
  * Musil  _et al._ (2021a) F. Musil, A. Grisafi, A. P. Bartók, C. Ortner, G. Csányi, and M. Ceriotti, Chemical Reviews 121, 9759 (2021a), arXiv:2101.04673 .
  * Unke  _et al._ (2021) O. T. Unke, S. Chmiela, H. E. Sauceda, M. Gastegger, I. Poltavsky, K. T. Schütt, A. Tkatchenko, and K. R. Müller, “Machine Learning Force Fields,”  (2021), arXiv:2010.07067 .
  * Izvekov and Voth (2005) S. Izvekov and G. A. Voth, Journal of Chemical Physics 123, 134105 (2005).
  * Noid  _et al._ (2008) W. G. Noid, J. W. Chu, G. S. Ayton, V. Krishna, S. Izvekov, G. A. Voth, A. Das, and H. C. Andersen, Journal of Chemical Physics 128, 244114 (2008).
  * Wang  _et al._ (2021) J. Wang, N. Charron, B. Husic, S. Olsson, F. Noé, and C. Clementi, Journal of Chemical Physics 154, 164113 (2021).
  * Bartók  _et al._ (2013) A. P. Bartók, R. Kondor, and G. Csányi, Physical Review B 87, 184115 (2013), arXiv:1209.3140v2 .
  * Darby  _et al._ (2021) J. P. Darby, J. R. Kermode, and G. Csányi (2021) arXiv:2112.13055 .
  * Musil  _et al._ (2021b) F. Musil, M. Veit, A. Goscinski, G. Fraux, M. J. Willatt, M. Stricker, T. Junge, and M. Ceriotti, Journal of Chemical Physics 154, 114109 (2021b), arXiv:2101.08814 .
  * Kapil  _et al._ (2019a) V. Kapil, M. Rossi, O. Marsalek, R. Petraglia, Y. Litman, T. Spura, B. Cheng, A. Cuzzocrea, R. H. Meißner, D. M. Wilkins, B. A. Helfrecht, P. Juda, S. P. Bienvenue, W. Fang, J. Kessler, I. Poltavsky, S. Vandenbrande, J. Wieme, C. Corminboeuf, T. D. Kühne, D. E. Manolopoulos, T. E. Markland, J. O. Richardson, A. Tkatchenko, G. A. Tribello, V. Van Speybroeck, and M. Ceriotti, Computer Physics Communications 236, 214 (2019a).
  * Paszke  _et al._ (2019) A. Paszke, S. Gross, F. Massa, A. Lerer, J. Bradbury, G. Chanan, T. Killeen, Z. Lin, N. Gimelshein, L. Antiga, A. Desmaison, A. Kopf, E. Yang, Z. DeVito, M. Raison, A. Tejani, S. Chilamkurthy, B. Steiner, L. Fang, J. Bai, and S. Chintala, in  _Advances in Neural Information Processing Systems 32_, edited by H. Wallach, H. Larochelle, A. Beygelzimer, F. d'Alché-Buc, E. Fox, and R. Garnett (Curran Associates, Inc., 2019) pp. 8024–8035.
  * Partridge and Schwenke (1997) H. Partridge and D. W. Schwenke, The Journal of Chemical Physics 106, 4618 (1997).
  * Drautz (2019) R. Drautz, Physical Review B 99, 014104 (2019).
  * Habershon  _et al._ (2009) S. Habershon, T. E. Markland, and D. E. Manolopoulos, The Journal of Chemical Physics 131, 024501 (2009).
  * Gao  _et al._ (2018) Z. Gao, N. Giovambattista, and O. Sahin, Scientific Reports 8, 6228 (2018).
  * Ceriotti  _et al._ (2011) M. Ceriotti, D. E. Manolopoulos, and M. Parrinello, The Journal of Chemical Physics 134, 084104 (2011).
  * Ceriotti and Manolopoulos (2012) M. Ceriotti and D. E. Manolopoulos, Physical Review Letters 109, 100604 (2012).
  * Mauger  _et al._ (2021) N. Mauger, T. Plé, L. Lagardère, S. Bonella, É. Mangaud, J.-P. Piquemal, and S. Huppert, The Journal of Physical Chemistry Letters 12, 8285 (2021).
  * Pérez and Tuckerman (2011) A. Pérez and M. E. Tuckerman, The Journal of Chemical Physics 135, 064104 (2011).
  * Kapil  _et al._ (2016b) V. Kapil, J. Behler, and M. Ceriotti, The Journal of Chemical Physics 145, 234103 (2016b).
  * Kapil  _et al._ (2019b) V. Kapil, J. Wieme, S. Vandenbrande, A. Lamaire, V. Van Speybroeck, and M. Ceriotti, Journal of Chemical Theory and Computation 15, 3237 (2019b).
  * Sugita and Okamoto (1999) Y. Sugita and Y. Okamoto, Chemical Physics Letters 314, 141 (1999).
  * Tuckerman  _et al._ (1992) M. Tuckerman, B. J. Berne, and G. J. Martyna, The Journal of Chemical Physics 97, 1990 (1992).
  * Hall and Dowling (1967) R. T. Hall and J. M. Dowling, The Journal of Chemical Physics 47, 2454 (1967).
  * Lin  _et al._ (2010) L. Lin, J. A. Morrone, R. Car, and M. Parrinello, Physical Review Letters 105, 110602 (2010).
  * Andreani  _et al._ (2016) C. Andreani, G. Romanelli, and R. Senesi, The Journal of Physical Chemistry Letters 7, 2216 (2016).
  * Gasparotto  _et al._ (2018) P. Gasparotto, R. H. Meißner, and M. Ceriotti, Journal of Chemical Theory and Computation 14, 486 (2018).
  * Rossi  _et al._ (2014b) M. Rossi, H. Liu, F. Paesani, J. Bowman, and M. Ceriotti, The Journal of Chemical Physics 141, 181101 (2014b).
  * Liu  _et al._ (2016) H. Liu, Y. Wang, and J. M. Bowman, The Journal of Physical Chemistry B 120, 1735 (2016).
  * Paesani and Voth (2010) F. Paesani and G. A. Voth, The Journal of Chemical Physics 132, 014105 (2010).
  * Loose  _et al._ (2022) T. D. Loose, P. G. Sahrmann, and G. A. Voth, “Centroid Molecular Dynamics Can Be Greatly Accelerated Through Neural Network Learned Centroid Forces Derived from Path Integral Molecular Dynamics,”  (2022), arXiv:2208.07973 [physics].
  * Wu  _et al._ (2022) C. Wu, R. Li, and K. Yu, Frontiers in Molecular Biosciences 9 (2022).
  * Verma  _et al._ (2019) D. Verma, R. M. P. Tanyag, S. M. O. O’Connell, and A. F. Vilesov, Advances in Physics: X 4, 1553569 (2019).


# Standard candles and sirens rescue $H_{0}$

Aniket Agrawal11affiliationmark:  , Teppei Okumura11affiliationmark:
22affiliationmark:  , and Toshifumi Futamase33affiliationmark:  1 Institute of
Astronomy and Astrophysics, Academia Sinica, 11F of AS/NTU Astronomy-
Mathematics Building, No.1, Sec. 4, Roosevelt Road, Taipei 10617, Taiwan,
R.O.C<EMAIL_ADDRESS>
2 Kavli Institute for the Physics and Mathematics of the Universe (WPI),
UTIAS, The University of Tokyo, Kashiwa, Chiba 277-8583, Japan
3 Department of Astrophysics and Meteorology, Kyoto Sangyo University, Kita-
ku, Kyoto 603-8555, Japan

###### Abstract

We show that future observations of binary neutron star systems with
electromagnetic counterparts together with the traditional probes of low- and
high-redshift Type Ia supernovae (SNe Ia) can help resolve the Hubble tension.
The luminosity distance inferred from these probes and its scatter depend on
the underlying cosmology. By using the gravitational lensing of light or
gravitational waves emitted by, and peculiar motion of, these systems we
derive constraints on the sum of neutrino masses, the equation of state of
dark energy parametrized in the form $w_{0}+w_{a}(1-a)$, along with the Hubble
constant and cold dark matter density in the universe. We show that even after
marginalizing over poorly constrained physical quantities, such as the sum of
neutrino masses and the nature of dark energy, low-redshift gravitational-wave
observations, in combination with SNe Ia, have the potential to rule out new
physics as the underlying cause of the Hubble tension at $\gtrsim 5.5\sigma$.

###### Subject headings:

Cosmological parameters (339), Dark energy (351), Neutrino masses (1102),
Gravitational waves (678), Type Ia supernovae (1728)

††slugcomment: Accepted to The Astrophysical Journal: 10/09/2020

##  1\. Introduction

The standard model of cosmology is under stress. There is at least a $4\sigma$
discrepancy between the values of the expansion rate of the universe today,
the Hubble constant $h\equiv H_{0}/(100$ km s-1Mpc-1), measured using early
and late universe probes (Verde et al., 2019). A solution to this discrepancy
will either revolutionize our understanding of the physical world or help us
discover previously unknown systematics. Therefore, significant theoretical
and observational effort is currently underway to uncover the cause of this
‘Hubble tension’. On the theory side, several models have been proposed to
explain the origin of this discrepancy, from early dark energy (Poulin et al.,
2019), self-interacting neutrinos (Kreisch et al., 2020), or scalar fields
that inject energy locally around the matter-radiation equality (Agrawal et
al., 2019b), to effects of a local inhomogeneity (Kasai & Futamase, 2019).
Since there is still disagreement as to the theoretical underpinnings of any
of these models, another approach is to make better measurements of $h$ using
as many probes as possible. Different systematics affecting different probes
allow a check on systematics as an origin of the Hubble tension. It is shown
in Verde et al. (2019) that low-redshift probes such as Type Ia supernovae
(SNe Ia), strong lenses, water masers, and surface brightness fluctuations
seem to be consistent with each other and give a value $h\sim 0.73$. High-
redshift probes including the cosmic microwave background (CMB) measurements
from Planck (Planck Collaboration et al., 2020) and those obtained from
combining galaxy clustering with early universe Big Bang nucleosynthesis also
agree with each other and predict $h\sim 0.67$. There is, thus, a disagreement
between the values measured at high and low redshift.

Direct detection of gravitational waves (GWs) has opened up a new window into
the universe (Abbott et al., 2016). By measuring the time variation of GWs we
can measure the luminosity distance to their sources. These distances are
poised to be measured with extremely high precision. Optical follow-up of
sources that also emit light allows us to measure their redshift (Abbott et
al., 2017). Technological limits restrict us to detecting GWs from low-
redshift sources only. Nevertheless, they provide a completely independent way
to measure $h$ at these redshifts. In this paper we show that by combining
future measurements of $h$ using GW sources out to $z\lesssim 0.1$ and SNe Ia
out to $z\lesssim 1.7$, we can potentially resolve the Hubble tension by
providing extremely tight constraints on $h$ in these two redshift ranges. As
we show, the two values are expected to disagree by more than $5.5\sigma$
allowing us to conclude that at least some of the Hubble tension originates
from systematic uncertainties and not from any physical effects. This
conclusion is robust to any assumptions about the nature of dark energy or the
sum of neutrino masses that we marginalize over.

##  2\. The Magnitude-Redshift relation

Cosmological information from both standard candles (SNe Ia), and standard
sirens (GW sources), is contained in the luminosity distance,
$d_{\text{L}}(z)=(1+z)\chi(z)\,,$ where $\chi(z)$ is the comoving distance at
the same redshift,

$\displaystyle\chi(z)=\frac{1}{H_{0}}\int_{0}^{z}dz^{\prime}\frac{1}{E(z^{\prime})}\,,$
(1)

with

$\displaystyle E^{2}(z)=$
$\displaystyle\Omega_{r}(1+z)^{4}+\Omega_{M}(1+z)^{3}+\Omega_{K}(1+z)^{2}$
$\displaystyle+\Omega_{\Lambda}(1+z)^{3(1+w_{0}+w_{a})}e^{-3w_{a}z/(1+z)}\,,$
(2)

where $\Omega_{r},\,\Omega_{M},\,\Omega_{K}$, and $\Omega_{\Lambda}$ are the
energy density fractions of radiation, matter, curvature, and dark energy,
respectively, and $w\equiv w_{0}+w_{a}(1-a)$ is the time-varying equation of
state for dark energy, parameterized by $w_{0}$, which characterizes the
constant part, and $w_{a}$, which represents the amplitude of time variation
(Chevallier & Polarski, 2001; Linder, 2003). Conventionally, the observed
quantity is the apparent magnitude, $m(z)$, which is related to $d_{L}(z)$ as
$m(z)=5\,\text{log}_{10}d_{\text{L}}(z)+M\,,$ where $M$ is the absolute
magnitude of a source, which can be determined. For SNe Ia $M$ is calibrated
using their observed peak luminosity and its subsequent decay, while for GW
sources it is calculated using the spectral and temporal variation of the
emitted GWs. Thus, by measuring $m(z)$, and thus $d_{L}(z)$, one can constrain
the cosmological parameters.

These equations for luminosity distance hold for sources that have no peculiar
velocity and are observed in a homogeneous universe, such that the emitted
electromagnetic or GWs are not gravitationally lensed. In an inhomogeneous
universe, matter along the line of sight (l.o.s.) affects the propagation of
these waves via lensing and Sachs–Wolfe and integrated Sachs–Wolfe effects
(Sachs & Wolfe, 1967). In the matter domination era, lensing is the dominant
effect among these. Lensing changes the observed brightness of a given source,
making it appear either closer or farther than it actually is, thus changing
the observed luminosity distance.

Peculiar motion of the source affects the observed redshift. As a result, it
affects the predicted luminosity distance to the source. The total peculiar
motion is a sum of two components—one of cosmological origin, which is sourced
by the large-scale structure of the universe, and the other of astrophysical
origin, which is sourced by the small-scale dynamics of the host galaxy.
Cosmological information can only be gleaned from the first component, which
can be measured separately by measuring the peculiar velocity of the galaxy as
a whole. For example, if the host galaxy is part of a cluster we need to use
the cluster redshift (Léget et al., 2018). Thus, in the rest of the paper, we
consider only cosmological peculiar velocities.

As shown in Hui & Greene (2006) lensing and peculiar motion respectively alter
the inferred luminosity distance by

$\displaystyle\delta d_{\text{L, lens}}(z_{s},\hat{\textbf{n}})=$
$\displaystyle-\frac{3H^{2}_{0}\Omega_{m0}}{2}\int_{0}^{\chi_{s}}d\chi\frac{\chi(\chi_{s}-\chi)}{\chi_{s}}$
$\displaystyle\times(1+z)\delta_{m}(z,\hat{\textbf{n}})\,,$ (3)
$\displaystyle\delta d_{\text{L, vel}}(z_{s},\hat{\textbf{n}})=$
$\displaystyle\Bigg{[}\frac{\textbf{v}_{o}\cdot\hat{\textbf{n}}-\textbf{v}_{s}\cdot\hat{\textbf{n}}}{a_{s}H_{s}\chi_{s}}\Bigg{]}+\textbf{v}_{s}\cdot\hat{\textbf{n}}\,.$
(4)

Here $\hat{\textbf{n}}$ is the unit vector in the _observed_ l.o.s. direction,
$\chi_{s}$ is the comoving distance at _observed_ redshift $z_{s}$ of the
source, $a_{s}$ is the scale factor corresponding to $z_{s}$, $H_{s}$ is the
Hubble rate at redshift $z_{s}$, $\delta_{m}(z,\hat{\textbf{n}})$ is the
matter density fluctuation at redshift $z(\chi)$ in the direction
$\hat{\textbf{n}}$, and $\textbf{v}_{o}$ and $\textbf{v}_{s}$ are the peculiar
velocities of the observer and the host galaxy, respectively. These equations
follow from linearizing the Einstein equation for small metric perturbations,
so that terms of higher order can be dropped. Despite the linear approximation
it can still be used to account for some nonlinearity in density and velocity
perturbations because they are second and first derivatives of the metric
perturbation, respectively (assuming a linear relation between density and
velocity). It is also noteworthy that propagating GWs obey the _same_
equations as photons, as long as their amplitude is small and their wavelength
is long enough so that the ray optics limit is realized (Misner et al., 1973).

The observed magnitude changes by (Hada & Futamase, 2014, 2016, 2019)

$\displaystyle\delta m_{\text{obs}}(z,\hat{\textbf{n}})=\,$ $\displaystyle
5\,\text{log}_{10}(1+\delta
d_{\text{L}}(z,\hat{\textbf{n}}))\simeq\frac{5}{\ln\,10}\,\delta
d_{\text{L}}(z,\hat{\textbf{n}})\,,$ (5)

where we have assumed that the fluctuation in luminosity distance is small.
The variance in the luminosity distance can then be written as

$\displaystyle\left\langle\delta
d^{2}_{\text{L}}(z,\hat{\textbf{n}})\right\rangle=\left\langle\delta
d^{2}_{\text{L,lens}}(z,\hat{\textbf{n}})\right\rangle+\left\langle\delta
d^{2}_{\text{L,vel}}(z,\hat{\textbf{n}})\right\rangle\,.$ (6)

The cross-correlation between lensing and peculiar velocities vanishes because
peculiar velocities are integrated along the l.o.s. due to the lensing kernel
and so average out to zero (Hui & Greene, 2006). The cross-correlation between
the l.o.s. peculiar velocities of the host galaxy and the observer is
negligible compared to the variance of the peculiar velocity of the host
galaxy (Hui & Greene, 2006). We choose the redshift bins to be small enough
that the variances do not change appreciably within each bin. The lensing
contribution to the variance is then given as

$\displaystyle\sigma_{\text{lens}}^{2}(z,\hat{\textbf{n}})\equiv$
$\displaystyle\Big{[}\frac{5}{\ln 10}\Big{]}^{2}\left\langle\delta
d^{2}_{\text{L,lens}}(z,\hat{\textbf{n}})\right\rangle$ $\displaystyle=$
$\displaystyle\Big{[}\frac{15H^{2}_{0}\Omega_{m0}}{2\,\ln
10}\Big{]}^{2}\int_{0}^{\chi_{s}}d\chi\Big{[}\frac{\chi(\chi_{s}-\chi)}{\chi_{s}}\Big{]}^{2}$
$\displaystyle\times(1+z)^{2}\int\frac{dk}{2\pi}kP_{\text{nl}}(k,z)\,,$ (7)

where we have used Limber’s approximation (see Appendix D of Hui & Greene
(2006) for more details) and $P_{\text{nl}}(k,z)$ is the nonlinear matter
power spectrum at redshift $z$. The velocity contribution is given by

$\displaystyle\sigma_{\text{vel}}^{2}(z,\hat{\textbf{n}})\equiv$
$\displaystyle\Big{[}\frac{5}{\ln 10}\Big{]}^{2}\left\langle\delta
d^{2}_{\text{L,vel}}(z,\hat{\textbf{n}})\right\rangle$
$\displaystyle=\Big{[}\frac{5}{\ln
10}\Big{]}^{2}\Big{[}1-\frac{1}{a_{s}H_{s}\chi_{s}}\Big{]}^{2}$
$\displaystyle\times\int\frac{dk}{6\pi^{2}}[D^{{}^{\prime}}(k,z)]^{2}P_{\text{nl}}(k,z=0)\,,$
(8)

where $D^{\prime}(k,z)\equiv-H(z)\frac{d\,D(k,z)}{dz}$ with $D(k,z)$ the
linear growth factor for matter fluctuations. Note that it is a function of
wavenumber because of nonzero neutrino mass.

Figure 1.— (Left) Logarithmic derivatives of the predicted apparent magnitude,
$m_{\text{homo}}$ (top), (scaled) velocity variance $\sigma^{2}_{\text{vel}}$
(bottom left), and (scaled) lensing variance $\sigma^{2}_{\text{lens}}$
(bottom right), w.r.t. the five cosmological parameters of interest,
$\omega_{\text{CDM}}$ (red dotted), sum of neutrino masses $m_{\nu}$ (blue
dashed), dark energy equation of state parameters $w_{0}$ (green
dotted–dashed) and $w_{a}$ (yellow dotted–dotted–dashed), and the Hubble
constant $h$ (black solid). Note that for $w_{a}$ the standard derivative that
has been shown as a logarithm of $0$ is not defined. The black dashed line in
the upper panel corresponds to $z=0.1$, which is the highest redshift for our
low-$z$ sample, and below which the velocity variance dominates over lensing.
(Right) Expected distribution of standard sirens and SNe Ia from ZTF (Graham
et al., 2019) (left) and SNe Ia from WFIRST (Spergel et al., 2015) (right).

We use the Fisher matrix formalism as described in Agrawal et al. (2019a) to
make forecasts for cosmological constraints using standard sirens and candles.
The likelihood is taken to be Gaussian, and different events are assumed to be
independent of each other. The resultant Fisher matrix becomes (Agrawal et
al., 2019a)

$\displaystyle\mathcal{F}_{ab,\text{tot}}=\sum_{i=1}^{N}$
$\displaystyle\Bigg{[}\frac{1}{\sigma^{2}_{\text{tot}}}\frac{\partial
m_{\text{homo}}(z_{i})}{\partial\theta_{a}}\frac{\partial
m_{\text{homo}}(z_{i})}{\partial\theta_{b}}$
$\displaystyle+\frac{1}{2\,\sigma^{4}_{\text{tot}}}\frac{\partial\sigma^{2}_{\text{tot}}(z_{i})}{\partial\theta_{a}}\frac{\partial\sigma^{2}_{\text{tot}}(z_{i})}{\partial\theta_{b}}\Bigg{]}\,.$
(9)

where the sum runs over different observed events, $\theta_{a}$ denotes the
cosmological parameter of interest, $m_{\text{homo}}$ is the magnitude
predicted in a homogeneous universe and $\sigma^{2}_{\text{tot}}$ is the total
variance of the observed magnitude, which is given by the quadrature sum of
the lensing, peculiar velocity and intrinsic contributions,
$\sigma^{2}_{\text{tot}}=\sigma^{2}_{\text{lens}}+\sigma^{2}_{\text{vel}}+\sigma^{2}_{\text{in}}\,.$
We use $\sigma_{\text{in}}=0.02$ for GW sources (Holz & Hughes, 2005) and
$\sigma_{\text{in}}=0.12$ for SNe Ia (Spergel et al., 2015).

The top panel of Figure 1 shows the absolute value of derivatives of
$m_{\text{homo}}$, $\sigma^{2}_{\text{lens}}$, and $\sigma^{2}_{\text{vel}}$
that enter in the Fisher matrix, Equation (2), w.r.t. five cosmological
parameters that we wish to constrain—$h$, cold dark matter density at $z=0$,
$\omega_{CDM}\equiv\Omega_{CDM}h^{2}$, $w_{0}$, $w_{a}$, and the sum of
neutrino masses that we denote by $m_{\nu}$. The variances have been scaled by
the total variance (with $\sigma_{\text{in}}=0.12$) to make them equivalent to
the homogeneous magnitude in Equation (2). The first thing to notice is the
relative sizes of these derivatives. At low redshifts the derivative of the
homogeneous magnitude dominates the lensing contribution for all parameters,
whereas at higher redshifts the derivatives of the lensing variance start to
become comparable to that of the homogeneous magnitude. In contrast, for all
parameters other than $h$ the velocity contribution dominates at the lowest
redshifts, but it drops quite sharply as redshift increases. We have verified
that the constraints are in fact dominated by the homogeneous part at low
redshifts. However, including variances helps to break degeneracies among
parameters, and to put much tighter constraints on the dark energy equation of
state and sum of neutrino masses as shown in Agrawal et al. (2019a) and Hada &
Futamase (2019, 2016).

The derivative of the homogeneous magnitude w.r.t. $h$ is almost an order of
magnitude higher than that w.r.t. other parameters, and is $\sim 2$. Thus, at
the lowest redshift, even without including information from the lensing or
velocity variance, we obtain a fractional error on $h$ (with all other
parameters fixed) of $\sim 1$% from just one GW event of the precision
considered. This, while being exceptionally promising, should not be
surprising since a perfect measurement of the magnitude corresponds to a
perfect measurement of the expansion rate, which at $z\sim 0$ is just the
Hubble constant.

Figure 2.— (Left) Joint constraints on $\omega_{\text{CDM}}$, $m_{\nu}$,
$w_{0}$, and $w_{a}$ with all other cosmological parameters fixed. Red dashed
lines denote constraints obtained when using high-$z$ supernovae from WFIRST
and low-$z$ supernovae from ZTF, while blue filled regions show constraints
obtained when supernovae are combined with low-$z$ GW events. (Right) Same as
the left panel, but with $h$ also varying now. As shown, a constraint of $\sim
0.3$% on $h$ can be obtained by combining GW sources with SNe Ia, even when
the sum of neutrino masses and time-varying dark energy equation of state are
simultaneously varied.

##  3\. Results

Future experiments such as the Einstein Telescope (Punturo et al., 2010) will
observe GW events from sources such as binary black hole (BH-BH), black hole-
neutron star (BH-NS) or neutron star-neutron star (NS-NS) mergers. Of these,
electromagnetic radiation is emitted along with gravitational radiation for
BH-NS and NS-NS mergers. By optical follow-up a redshift can be measured for
these events (Abbott et al., 2017). Alternatively, one can cross-correlate
measurements of GW events in luminosity distance space with measurements of
galaxy distribution in the same space, to obtain the most likely standard
siren redshift (Kopparapu et al., 2008). Using these redshifts in combination
with the precise measurement of luminosity distance for these standard sirens
provides extremely tight constraints on fundamental cosmology as we now
demonstrate.

Future experiments will observe $2810$ GW events per $\text{Gpc}^{3}$ (Abbott
et al., 2020). We assume that they are uniformly distributed with this number
density out to $z=0.1$ and that redshifts for these events can be perfectly
measured. We scale this number by $1/(1+z)^{2}$ to account for the decrease in
optical flux with redshift. Then, the number of such events observed in each
redshift bin of $\Delta z=0.005$ out to $z=0.1$ is shown in the left plot of
the right panel of Figure 1.

For the supernova sample, we use ZTF, which will observe low-redshift
supernovae in the range $0.01\lesssim z\lesssim 0.1$ (Graham et al., 2019) and
WFIRST, which will observe high-redshift supernovae in the range $0.2\lesssim
z\lesssim 1.7$ (Spergel et al., 2015). The left plot of the right panel of
Figure 1 shows the expected distribution of supernovae in different redshift
bins for ZTF as the light green bars. The right plot of the right panel of
Figure 1 shows the expected distribution of supernovae in different redshift
bins for WFIRST. Note that while it is dominated by supernovae in the redshift
range $z\leq 0.6$, where lensing effects are subdominant compared to the
intrinsic uncertainty, the overall number of supernovae for $z\leq 0.6$ and
$z>0.6$ is approximately the same so that lensing effects must be taken into
account.

Using the above distributions in Equation (2) we can determine the expected
constraints on cosmological parameters using just the supernovae or combining
the information from both standard candles and low-redshift standard sirens.
Note that for high-redshift events the scatter in luminosity distance due to
lensing dominates over that from peculiar motion, while for low-redshift
events the one from peculiar motion dominates (Agrawal et al., 2019a).
Therefore, we do not include the scatter from lensing at low redshifts and the
scatter from peculiar velocities at high redshifts.

The top panel of Figure 2 shows the constraints obtained using these two
probes in the 4-parameter space considered in Agrawal et al. (2019a). The red
dashed contours represent constraints obtained when considering only
information from supernovae. We include effects only from lensing due to
matter along the l.o.s. for the high-redshift probes as peculiar motion has
negligible contribution to the luminosity distance for these redshifts (Hui &
Greene, 2006). The blue filled regions show constraints obtained when we
combine supernovae with low-redshift GW events. Note that for low-redshift
probes we include the contribution from peculiar motion alone as lensing is
not significant at low redshifts (Hui & Greene, 2006).

Equation (2) shows that cosmological information in the apparent magnitude
neatly splits into contributions coming from the homogeneous universe and
inhomogeneous universe, the first and second terms, respectively. For the
homogeneous contribution, $\omega_{\text{CDM}}$ and $m_{\nu}$ are completely
degenerate, both contributing via the total matter density only. Using
information about the inhomogeneities allows us to break this degeneracy.
Since the homogeneous contribution dominates, the shape and orientation of
contours for these parameters are similar but not exactly the same. We have
verified that if we only consider the homogeneous contribution, while fixing
one of these parameters, then the contours do indeed look the same.

We find that since the luminosity distance is measured to much higher accuracy
for a GW source than for an SN Ia, only about 600 GW events add significant
information about the sum of neutrino masses over the 4000 SNe Ia from ZTF and
WFIRST. This increased constraining power comes from the larger contribution
of peculiar velocities once the intrinsic uncertainty becomes lower. These
tight constraints also point to the possibility of allowing other parameters
to be free. In particular, with an eye to the Hubble tension we consider
letting the Hubble constant, $h$, be free.

The resulting constraints in the 5-parameter space are shown in the bottom
panel of Figure 2. We see a drastic reduction in the constraints on neutrino
mass, with the $2\sigma$ constraint from including GW events almost as good as
the $1\sigma$ constraint from just supernovae. The other significant reduction
is seen in $h$ as expected. We find that one can constrain $h$ to better than
0.3% by combining GW events with the supernova measurements even when the
time-varying dark energy equation of state and the sum of neutrino masses are
simultaneously varied. A key reason for this reduction when using GWs is that
they allow us to measure $h$ with extremely high precision leading to an
almost fixed $h$. Thus the blue contours are not much worse than the same in
the top panel of Figure 2.

Figure 3.— Constraints in $h$-$\omega_{\text{CDM}}$ space after marginalizing
over the sum of neutrino masses and dark energy equation of state parameters
($w_{0}$, $w_{a}$). Red contours show constraints obtained when using only
supernovae from ZTF and WFIRST, while blue contours show constraints obtained
when using only low-$z$ GW events. The clear separation of blue and red
contours presents an opportunity to resolve the Hubble tension.

In order to emphasize the promise of GW observations further, we show the
constraints on $h$ obtained when marginalizing over all other parameters
except $\omega_{\text{CDM}}$ in Figure 3. The blue contours show constraints
obtained with a fiducial value of $h=0.6727$, consistent with CMB measurements
(Ade et al., 2016) and using low-redshift GW events as the probe. The red
contours show constraints from SNe Ia from both low (ZTF) and high (WFIRST)
redshifts with $h=0.733$ as the fiducial value, consistent with local
measurements (Verde et al., 2019). Note that these contours are impossible to
draw considering only the homogeneous contribution as $\omega_{\text{CDM}}$ is
completely degenerate with $m_{\nu}$ and so no constraints can be obtained if
either is marginalized over. It is only with peculiar velocity information
that we can constrain one after marginalizing over the other.

If the Hubble tension arises from unknown systematics, GW and SNe Ia
observations should see different values of $h$. In that case it is crucial
that we can conclusively say that these two probes are inconsistent with each
other, for which the errors on each measurement should be small enough to
completely exclude the other. As Figure 3 demonstrates for the probes
considered here, this is the case. The two fiducial values we adopt for $h$
are more than $5.5\sigma$ separated from each other (after marginalizing over
$\omega_{\text{CDM}}$). Alternatively, if the two contours shown in Figure 3
overlapped, despite different central values, we could not conclude that they
are inconsistent with each other and the Hubble tension could be said to still
originate from either new physics or unknown systematics. By having
nonoverlapping contours we are able to ensure that the difference arises from
systematic errors and thus rule out new physics. If, on the other hand, we
used $h=0.71$ and $0.68$ for low and high redshifts, respectively, we still
obtain a $2.5\sigma$ separation between the two contours, which can be further
improved by observation of a larger number of GW events or by using priors
from other measurements. It is also important to note that these conclusions
follow even after marginalizing over the sum of neutrino masses or dark energy
equation of state parameters, which are poorly constrained by the CMB, or by
baryon acoustic oscillation measurements in the absence of a CMB prior. They
are thus robust to changing dark energy models whose equation of state can be
written in the form $w=w_{0}+w_{a}(1-a)$ or the sum of neutrino masses.

##  4\. Conclusion

We have shown that a combination of future GW and SNe Ia observations can
resolve the Hubble tension if the universe is assumed to be homogeneous and
isotropic on scales spanning the supernovae and GW events we have considered.
The numbers considered here can be realized in the next decade. We find that
only about 600 GW sources with a measurement of their associated redshifts out
to $z\sim 0.1$ are needed for $2\sigma$ constraints on $h\sim 0.3$%. These
constraints are obtained even without assuming any prior on the sum of
neutrino masses and the dark energy equation of state. Peculiar motion of GW
sources is indispensable to break the degeneracy between the sum of neutrino
masses and the cold dark matter density in the universe and to obtain tight
constraints on the, as yet undetermined, sum of neutrino masses. The Hubble
tension presents a fantastic opportunity to test our current understanding of
the universe. Future observations made with GWs from inspiraling binary
neutron star systems and light from SNe Ia will clarify the origin of the
Hubble tension, whether it is coming from new physics or unknown systematics
in a homogeneous universe, or that the universe is inhomogeneous.

T.O. acknowledges support from the Ministry of Science and Technology of
Taiwan under grants No. MOST 106-2119-M-001-031-MY3 and No. MOST
109-2112-M-001-027- and the Career Development Award, Academia Sinica (AS-
CDA-108-M02) for the period of 2019 to 2023. T.F. is supported by Grant-in-
Aids for Scientific Research from JSPS (No. 17K05453 and No. 18H04357).

## References

  * Abbott et al. (2016) Abbott, B., et al. 2016, Phys. Rev. Lett., 116, 061102
  * Abbott et al. (2017) —. 2017, Phys. Rev. Lett., 119, 161101
  * Abbott et al. (2020) —. 2020, Astrophys. J. Lett., 892, L3
  * Ade et al. (2016) Ade, P. A. R., et al. 2016, Astron. Astrophys., 594, A13
  * Agrawal et al. (2019a) Agrawal, A., Okumura, T., & Futamase, T. 2019a, Phys. Rev. D, 100, 063534
  * Agrawal et al. (2019b) Agrawal, P., Cyr-Racine, F.-Y., Pinner, D., & Randall, L. 2019b, arXiv:1904.01016
  * Chevallier & Polarski (2001) Chevallier, M., & Polarski, D. 2001, Int. J. Mod. Phys., D10, 213
  * Graham et al. (2019) Graham, M. J., Kulkarni, S. R., Bellm, E. C., et al. 2019, PASP, 131, 078001
  * Hada & Futamase (2014) Hada, R., & Futamase, T. 2014, JCAP, 1412, 042
  * Hada & Futamase (2016) —. 2016, Astrophys. J., 828, 112
  * Hada & Futamase (2019) Hada, R., & Futamase, T. 2019, J. Cosmology Astropart. Phys, 2019, 033
  * Holz & Hughes (2005) Holz, D. E., & Hughes, S. A. 2005, Astrophys. J., 629, 15
  * Hui & Greene (2006) Hui, L., & Greene, P. B. 2006, Phys. Rev., D73, 123526
  * Kasai & Futamase (2019) Kasai, M., & Futamase, T. 2019, PTEP, 2019, 073E01
  * Kopparapu et al. (2008) Kopparapu, R. K., Hanna, C., Kalogera, V., et al. 2008, Astrophys. J., 675, 1459
  * Kreisch et al. (2020) Kreisch, C. D., Cyr-Racine, F.-Y., & Doré, O. 2020, Phys. Rev. D, 101, 123505
  * Léget et al. (2018) Léget, P. F., et al. 2018, Astron. Astrophys., 615, A162
  * Linder (2003) Linder, E. V. 2003, Physical Review Letters, 90, 091301
  * Misner et al. (1973) Misner, C. W., Thorne, K., & Wheeler, J. 1973, Gravitation (San Francisco: W. H. Freeman)
  * Planck Collaboration et al. (2020) Planck Collaboration, Aghanim, N., Akrami, Y., et al. 2020, A&A, 641, A6
  * Poulin et al. (2019) Poulin, V., Smith, T. L., Karwal, T., & Kamionkowski, M. 2019, Phys. Rev. Lett., 122, 221301
  * Punturo et al. (2010) Punturo, M., et al. 2010, Class. Quant. Grav., 27, 194002
  * Sachs & Wolfe (1967) Sachs, R. K., & Wolfe, A. M. 1967, ApJ, 147, 73
  * Spergel et al. (2015) Spergel, D., et al. 2015, arXiv:1503.03757
  * Verde et al. (2019) Verde, L., Treu, T., & Riess, A. G. 2019, Nature Astronomy, 3, 891