# Physics-Informed Machine Learning of Dynamical Systems for Efficient
Bayesian Inference

Somayajulu L. N. Dhulipala
Idaho National Laboratory
<EMAIL_ADDRESS>
&Yifeng Che
Idaho National Laboratory
<EMAIL_ADDRESS>
Michael D. Shields
Johns Hopkins University
<EMAIL_ADDRESS>

###### Abstract

Although the no-u-turn sampler (NUTS) is a widely adopted method for
performing Bayesian inference, it requires numerous posterior gradients which
can be expensive to compute in practice. Recently, there has been a
significant interest in physics-based machine learning of dynamical (or
Hamiltonian) systems and Hamiltonian neural networks (HNNs) is a noteworthy
architecture. But these types of architectures have not been applied to solve
Bayesian inference problems efficiently. We propose the use of HNNs for
performing Bayesian inference efficiently without requiring numerous posterior
gradients. We introduce latent variable outputs to HNNs (L-HNNs) for improved
expressivity and reduced integration errors. We integrate L-HNNs in NUTS and
further propose an online error monitoring scheme to prevent sampling
degeneracy in regions where L-HNNs may have little training data. We
demonstrate L-HNNs in NUTS with online error monitoring considering several
complex high-dimensional posterior densities and compare its performance to
NUTS.

##  1 Introduction

Hamiltonian Monte Carlo (HMC) is a popular sampling algorithm for scalable
Bayesian inference considering complex and/or high-dimensional posterior
distributions [8]. Given the initial position (or posterior sample) and
momenta, HMC simulates the trajectory of a particle using Hamiltonian dynamics
with gradients of the target posterior and its position at the end time ($T$)
constitutes a new sample from the posterior. Therefore, the choice of this end
time considerably effects the posterior inference quality, and for complex
and/or high-dimensional posterior it can be difficult to manually specify the
optimal end time [4]. No-u-turn sampling (NUTS) [4] has been introduced to
automatically determine the optimal end time. If
$\\{\boldsymbol{q},\boldsymbol{p}\\}$ denote a position-momenta (i.e.,
$\boldsymbol{q}-\boldsymbol{p}$) pair, NUTS builds a binary tree of
$\\{\boldsymbol{q},\boldsymbol{p}\\}$ states and selects a subset of these
states with the aid of slice sampling [7] so as to not violate detailed
balance. For each sample, NUTS terminates the tree building procedure whenever
a u-turn is encountered. NUTS may also terminate the tree building procedure
when the integration error while simulating the particle trajectory is large.
NUTS has been shown to result in a better performance over HMC in terms of the
effective sample size (ESS) [10] per gradient evaluation by several studies
(e.g., [9]). However, it still requires numerous gradients estimations of the
target posterior which can be expensive in practice, especially, when dealing
with large datasets or computational models [2].

Recently, physics-informed machine learning of dynamical (or Hamiltonian)
systems has been receiving a significant interest owing to the wide variety of
applications in the physical sciences. Greydanus et al. [3] introduced
Hamiltonian Neural Networks (HNNs) to learn Hamiltonian systems in an
unsupervised fashion by virtue of a physics-informed loss function. HNNs
satisfy important properties such as time reversibility and Hamiltonian
conservation and were shown to be superior than data-driven neural networks.
Along these lines, several studies proposed improved architectures for
learning Hamiltonian systems such as physics-informed HNNs [6], symplectic
ODE-Net [11], and symplectic neural nets [5]. However, these efforts were all
focused on learning Hamiltonian systems for classical mechanics problems like
n-pendulum behavior and planetary motion, and not to perform Bayesian
inference efficiently.

We propose latent Hamiltonian neural networks (L-HNNs) in NUTS for efficiently
solving Bayesian inference problems by not requiring numerous posterior
gradient estimations [1]. For this, we first introduce latent variable outputs
in HNNs (i.e., L-HNNs) for improved expressivity and reduced integration
errors while simulating the Hamiltonian dynamics. Next, we integrate L-HNNs
with NUTS and propose an online error monitoring scheme to prevent sampling
degeneracy in regions of the uncertainty space where little training data may
be available for L-HNNs (for e.g., tails of the posterior). In essence, this
novel error monitoring scheme temporarily reverts to using posterior gradients
whenever the L-HNNs integration errors are large. We demonstrate L-HNNs in
NUTS with online error monitoring considering several complex posteriors and
compare its performance with traditional NUTS in terms of effectively
capturing the posterior features and the ESS per posterior gradient (i.e.,
computational efficiency).

##  2 Methods

L-HNNs are fully-connected feed-forward neural networks defined as:

$\displaystyle\boldsymbol{u}_{p}=\phi(\boldsymbol{w}_{p-1}~{}\boldsymbol{u}_{p-1}+\boldsymbol{b}_{p-1}),~{}~{}\\{p\in
1,\dots,P\\}$ (1)
$\displaystyle\boldsymbol{\lambda}=\boldsymbol{w}_{P}~{}\boldsymbol{u}_{P}+\boldsymbol{b}_{P}$

where $P$ is the number of hidden layers indexed by $p$, $\boldsymbol{u}_{p}$
are the outputs of the hidden layer $p$, $\boldsymbol{w}_{p}$ and
$\boldsymbol{b}_{p}$ are the weights and biases, respectively, and $\phi(.)$
is the nonlinear activation function. $\boldsymbol{u}_{0}$ are the inputs
defined by the position-momenta pair
$\boldsymbol{z}=\\{\boldsymbol{q},~{}\boldsymbol{p}\\}$. If $d$ is the
dimensionality of the uncertainty space, L-HNNs predict $d$ latent variables
defined by the vector $\boldsymbol{\lambda}$. The sum of these latent
variables is defined as the scalar Hamiltonian predicted by L-HNNs:

$H_{\boldsymbol{\theta}}=\sum_{i=1}^{d}\lambda_{i}$ (2)

Whereas, HNNs directly predict the scalar value of the Hamiltonian [3], L-HNNs
precit $d$ latent variables whose sum is defined as the scalar Hamiltonian.
Introduction of these latent variables improves expressivity and reduces the
integration errors when simulating the Hamiltonian trajectories [1]. The
training data is provided in terms of sets of
$\\{\boldsymbol{q},~{}\boldsymbol{p}\\}$ pairs. Time derivatives of these
$\\{\boldsymbol{q},~{}\boldsymbol{p}\\}$ pairs are computed. Gradients of the
Hamiltonians predicted by L-HNNs [Equation (2)] are also computed. Then, the
following unsupervised physics-based loss function is minimized:

$\mathcal{L}=\Big{\lVert}\frac{\partial
H_{\boldsymbol{\theta}}}{\partial\boldsymbol{p}}-\frac{\Delta\boldsymbol{q}}{\Delta
t}\Big{\lVert}+\Big{\lVert}-\frac{\partial
H_{\boldsymbol{\theta}}}{\partial\boldsymbol{q}}-\frac{\Delta\boldsymbol{p}}{\Delta
t}\Big{\lVert}$ (3)

The training data time derivatives $\frac{\Delta\boldsymbol{q}}{\Delta t}$ and
$\frac{\Delta\boldsymbol{p}}{\Delta t}$ are, respectively, equal to the
gradients of the exact Hamiltonians $\frac{\partial
H}{\partial\boldsymbol{p}}$ and $-\frac{\partial H}{\partial\boldsymbol{q}}$.
Thus, the loss function accounts for the physics governing Hamiltonian
systems.

For Bayesian inference problems, the Hamiltonian is usually the sum of
negative logarithm of the target posterior (i.e., potential energy) and
negative logarithm of the distribution of momenta (i.e., kinetic energy). In
most cases, the distribution of momenta is considered to be a standard
multivariate Gaussian. In NUTS, it is required to simulate the particle
trajectory given the initial $\boldsymbol{z}(0)$. Since the Hamilton’s
equations are first order ODEs, the particle trajectory is obtained by
numerically simulating the following integral with the aid of symplectic
integrators like the leap frog:

$\boldsymbol{z}(T)=\boldsymbol{z}(0)+\int_{0}^{T}\begin{bmatrix}\boldsymbol{0}_{d\times
d}&\boldsymbol{I}_{d\times d}\\\ -\boldsymbol{I}_{d\times
d}&\boldsymbol{0}_{d\times d}\end{bmatrix}~{}\nabla
H_{\boldsymbol{\theta}}(\boldsymbol{z})~{}dt$ (4)

A noteworthy component of this equation is the term $\nabla
H_{\boldsymbol{\theta}}(\boldsymbol{z})$ which is computed through L-HNNs and
not by using gradients of the target posterior.

For posterior sampling using NUTS, for each new sample, several trajectories
of a fictitious particle need to be simulated in sequence while doubling the
trajectory length each time as part of a tree building procedure. The primary
criterion for terminating the tree building procedure is when any of the sub-
trees makes a u-turn, which works towards reducing the serial correlations
between posterior samples. The other criterion for termination is when the
integration errors while simulating the trajectories are large, as defined by:

$\varepsilon\equiv H\big{(}\boldsymbol{z}\big{)}+\ln{u}>\Delta_{max}$ (5)

where $H\big{(}.\big{)}$ is the Hamiltonian value given a
$\boldsymbol{z}=\\{\boldsymbol{q},\boldsymbol{p}\\}$ pair, $u$ is the slice
value simulated from
$Uniform([0,~{}\exp{\\{-H\big{(}\boldsymbol{z}(0)\big{)}\\}}])$ given the
initial $\boldsymbol{z}(0)$ pair, and $\Delta_{max}$ is the error threshold
usually set to 1000. NUTS mostly employs the standard leap frog integrator
using gradients of the target posterior for simulating the particle
trajectories. While the replacement of the target posterior gradients with
L-HNNs gradients [as discussed in Equation 4] is possible and works well for
simple problems, for complex posterior spaces, this may lead to sampling
degeneracy. The reason for such degeneracy is the termination criterion in
Equation (5). Whenever NUTS with L-HNN gradients enters a region of the
posterior space where there was little to no training data for L-HNNs,
integration errors of the particle trajectories can be large, resulting in a
premature termination of the tree building procedure. As a result, there will
be clusters of samples at close vicinity in regions of the posterior space
where L-HNNs had little to no training data. These clusters of samples
constitute the sampling degeneracy problem.

To mitigate the sampling degeneracy problem, we propose an online error
monitoring scheme. We introduce two error thresholds, $\Delta_{max}^{hnn}$,
when using L-HNN gradients in the leap frog, and $\Delta_{max}^{lf}$, when
using the standard leap frog with target posterior gradients. We set
$\Delta_{max}^{hnn}<<\Delta_{max}^{lf}$; for example, $\Delta_{max}^{hnn}=10$
and $\Delta_{max}^{lf}=1000$. The particle trajectories in NUTS are, by
default, simulated using the L-HNN gradients. Whenever the integration errors
using L-HNNs, as denoted by $\varepsilon^{hnn}$, are greater than
$\Delta_{max}^{hnn}$, traditional leap frog with posterior gradients takes
over for a few samples $N^{lf}$. The value of $N^{lf}$ can be in between 5 and
20 and this essentially aids in NUTS returning to regions of high probability
density after which the particle trajectories can again be simulated using the
L-HNN gradients. Equation 6 summarized this error monitoring scheme while
using L-HNNs in NUTS:

$\displaystyle\varepsilon^{hnn}\leq\Delta_{max}^{hnn}$
$\displaystyle\to\textrm{use leap frog with L-HNN gradients}$ (6)
$\displaystyle\varepsilon^{hnn}>\Delta_{max}^{hnn},~{}\varepsilon^{lf}\leq\Delta_{max}^{lf}$
$\displaystyle\to\textrm{use leap frog with posterior gradients for $N^{lf}$
samples}$
$\displaystyle\varepsilon^{hnn}>\Delta_{max}^{hnn},~{}\varepsilon^{lf}>\Delta_{max}^{lf}$
$\displaystyle\to\textrm{terminate tree building; move to next sample}$

where $\varepsilon^{hnn}$ and $\varepsilon^{lf}$ are the integration errors
[equation (5)] when using leap frog with L-HNN and posterior gradients,
respectively.

##  3 Results

We demonstrate Bayesian inference using L-HNNs in NUTS with online error
monitoring considering several case studies. The performance is compared with
traditional NUTS which relies on the posterior gradients. First, we consider a
2-D eight Gaussian mixture density. Each Gaussian has the same covariance
which is an identity matrix but different mean vectors. Figure 1(a) presents
the Hamiltonian evolution with time for different initial
$\\{\boldsymbol{q},\boldsymbol{p}\\}$ consider LHNN-NUTS and NUTS. It is
noticed that LHNN-NUTS satisfactorily conserves the Hamiltonian values in a
similar fashion to that of NUTS. Figures 1(b) and 1(c) present the scatter
plot of $100,000$ simulated samples (first 5000 considered as burn-in) using
LHNN-NUTS and NUTS, respectively. LHNN-NUTS captures all the posterior modes
satisfactorily in a similar fashion to that of NUTS. Moreover, from a
computational perspective, while NUTS requires more than 15 Million posterior
gradients, LHNN-NUTS requires only 0.4 Million posterior gradients for the
training data and during online error monitoring. The ESS per posterior
gradient of LHNN-NUTS and NUTS is 0.0269 and 0.000707, respectively, which is
between 1-2 orders of magnitude in improvement.

(a)

(b)

(c)

Figure 1: A 2-D eight Gaussian mixture is considered as the target posterior.
(a) Hamiltonian conservation for random initial
$\\{\boldsymbol{q},\boldsymbol{p}\\}$ values considering L-HNN gradients and
posterior gradients. (b) Samples from the target posterior simulated using
LHNN-NUTS. (c) Samples from the target posterior simulated using traditional
NUTS.

We have also tested the performance of LHNN-NUTS considering several other
complex high-dimensional posteriors. Table 1 summarizes the results in terms
of the number of gradients and the ESS per posterior gradient and also
presents these values corresponding to NUTS for comparison. Overall, it can be
observed that LHNN-NUTS requires about 1-2 orders of magnitude less posterior
gradients than NUTS. Moreover, LHNN-NUTS results in about an order of
magnitude improvement in the ESS per posterior gradient compared to NUTS.
These results contribute to the scalability and efficiency of the proposed
LHNNs in NUTS with online error monitoring scheme relative to NUTS which is
widely adopted for Bayesian inference.

Table 1: Summary of the performance comparison between LHNN-NUTS and NUTS (adapted from Dhulipala et al. [1]). Posterior density | QOI | LHNN-NUTS | NUTS
---|---|---|---
2-D eight Gaussian mixture |  [c]# gradients |  |
ESS/gradient |  [c]0.4 Million |  |
0.0269 |  [c]15 Million |  |
0.000707 |  |  |
10-D Rosenbrock density |  [c]# gradients |  |
ESS/gradient |  [c]0.42 Million |  |
0.0368 |  [c]7 Million |  |
0.00219 |  |  |
24-D Bayesian logistic regression |  [c]# gradients |  |
ESS/gradient |  [c]0.4 Million |  |
0.243 |  [c]1.2 Million |  |
0.0777 |  |  |
100-D rough well |  [c]# gradients |  |
ESS/gradient |  [c]0.42 Million |  |
0.0065 |  [c]1.28 Million |  |
0.00373 |  |  |

##  4 Summary and conclusions

Physics-based neural networks that learn Hamiltonian systems (e.g., HNNs) have
been receiving a great deal of interest. We proposed using HNNs to efficiently
solve Bayesian inference problems by not requiring numerous posterior
gradients. To increase the expressivity of HNNs and reduce the integration
error while simulating Hamiltonian trajectories, we proposed HNNs with latent
outputs (L-HNNs). We further proposed to use L-HNNs in NUTS with an online
error monitoring scheme that reverts to using the posterior gradients for a
few samples whenever the L-HNNs integration errors are high, which may be in
regions of the uncertainty space where L-HNNs had little training data.
Considering several complex high-dimensional posteriors, we demonstrated LHNN-
NUTS with online error monitoring and compared its performance with NUTS.
Overall, LHNN-NUTS with online error monitoring required 1-2 orders of
magnitude lesser posterior gradients and resulted in about 1 order of
magnitude better ESS per gradient compared to NUTS.

## Acknowledgments

This research is supported through INL’s Laboratory Directed Research and
Development (LDRD) Program under U.S. Department of Energy (DOE) Idaho
Operations Office contract no. DE-AC07-05ID14517.

## Broader impact

Bayesian inference is a principal method for parameter calibration and
uncertainty quantification in many scientific and engineering fields. Since
closed-form solutions usually do not exist for practical Bayesian inference
problems, Markov chain Monte Carlo (MCMC) algorithms are employed to sample
from the target probability density. While random-walk MCMC algorithms are
popular, they have poor scalability with the number of dimensions, and can
feature large serial correlations between the samples. Therefore, algorithms
under the HMC framework (including NUTS) are often employed. However, the HMC-
based algorithms require numerous costly gradient estimations of the target
posterior and can be expensive in practice. The proposed LHNN-NUTS with online
error monitoring removes this computational hurdle and contributes to solving
practical Bayesian inference problems across science and engineering while
guaranteeing robustness.

## References

  * [1] S. L. N. Dhulipala, Y. Che, and M. D. Shields.  Bayesian inference with latent Hamiltonian neural networks, Aug. 2022\.  arXiv:2208.06120.
  * [2] S. L. N. Dhulipala, M. D. Shields, B. W. Spencer, C. Bolisetti, A. E. Slaughter, V. M. Laboure, and P. Chakroborty.  Active learning with multifidelity modeling for efficient rare event simulation.  Journal of Computational Physics, 468:111506, 2022.
  * [3] S. Greydanus, M. Dzamba, and J. Yosinski.  Hamiltonian neural networks, Sept. 2019.  arXiv:1906.01563.
  * [4] M. D. Hoffman and A. Gelman.  The No-U-Turn sampler: adaptively setting path lengths in Hamiltonian Monte Carlo.  Journal of Machine Learning Research, 15(1):1593–1623, 2014.
  * [5] P. Jin, Z. Zhang, A. Zhu, Y. Tang, and G. E. Karniadakis.  SympNets: Intrinsic structure-preserving symplectic networks for identifying Hamiltonian systems.  Neural Networks, 132:166–179, 2020.
  * [6] M. Mattheakis, D. Sondak, A. S. Dogra, and P. Protopapas.  Hamiltonian neural networks for solving equations of motion.  Physical Review E, 105:065305, 2022.
  * [7] R. M. Neal.  Slice sampling.  The Annals of Statistics, 31(3):705–767, 2003.
  * [8] R. M. Neal.  MCMC using Hamiltonian dynamics, June 2012.  arXiv:1206.1901.
  * [9] M. Nishio and A. Arakawa.  Performance of Hamiltonian Monte Carlo and No-U-Turn Sampler for estimating genetic parameters and breeding values.  Genetics Selection Evolution, 51:1–12, 2019.
  * [10] A. Vehtari, A. Gelman, D. Simpson, B. Carpenter, and P.-C. Burkner.  Rank-normalization, folding, and localization: An improved R-hat for assessing convergence of MCMC.  Bayesian Analysis, 16(2):667–718, 2021.
  * [11] Y. D. Zhong, B. Dey, and A. Chakraborty.  Symplectic ode-net: Learning Hamiltonian dynamics with control, Apr. 2019.  arXiv:1909.12077.

## Appendix A: Algorithms for L-HNNs in NUTS with online error monitoring

Algorithm 1 L-HNNs in NUTS with online error monitoring (main loop) [1]

1:Hamiltonian: $H=U(\boldsymbol{q})+K(\boldsymbol{p})$, Samples: M; Starting
sample: $\\{\boldsymbol{q}^{0},~{}\boldsymbol{p}^{0}\\}$; Step size: $\Delta
t$; Threshold for leapfrog: $\Delta_{max}^{lf}$; Threshold for L-HNNs:
$\Delta_{max}^{hnn}$; Number of leapfrog samples: $N_{lf}$

2:Initialize $\boldsymbol{1}_{lf}=0,~{}n_{lf}=0$

3:for $i=1:M$ do

4: $\boldsymbol{p}(0)\sim\mathcal{N}(\boldsymbol{0},~{}\boldsymbol{I}_{d})$

5: $\boldsymbol{q}(0)=\boldsymbol{q}^{i-1}$

6: $u\sim
Uniform\Big{(}\Big{[}0,~{}\exp{\\{-H\big{(}\boldsymbol{q}(0),~{}\boldsymbol{p}(0)\big{)}\\}}\Big{]}\Big{)}$

7: Initialize
$\boldsymbol{q}^{-}=\boldsymbol{q}(0),\boldsymbol{q}^{+}=\boldsymbol{q}(0),\boldsymbol{p}^{-}=\boldsymbol{p}(0),\boldsymbol{p}^{+}=\boldsymbol{p}(0),j=0,\boldsymbol{q}^{*}=\boldsymbol{q}^{i-1},n=1,s=1$

8: if $\boldsymbol{1}_{lf}=1$ then

9: $n_{lf}\leftarrow n_{lf}+1$

10: end if

11: if $n_{lf}=N_{lf}$ then

12: $\boldsymbol{1}_{lf}=0,~{}n_{lf}=0$

13: end if

14: while $s=1$ do

15: Choose direction $\nu_{j}\sim Uniform(\\{-1,~{}1\\})$

16: if $j=-1$ then

17:
$\boldsymbol{q}^{-},\boldsymbol{p}^{-},\\_,\\_,\boldsymbol{q}^{\prime},\boldsymbol{p}^{\prime},n^{\prime},s^{\prime},\boldsymbol{1}_{lf}=\textrm{{BuildTree}}(\boldsymbol{q}^{-},\boldsymbol{p}^{-},u,\nu_{j},j,\Delta
t,\boldsymbol{1}_{lf})$

18: else

19:
$\\_,\\_,\boldsymbol{q}^{+},\boldsymbol{p}^{+},\boldsymbol{q}^{\prime},\boldsymbol{p}^{\prime},n^{\prime},s^{\prime},\boldsymbol{1}_{lf}=\textrm{{BuildTree}}(\boldsymbol{q}^{+},\boldsymbol{p}^{+},u,\nu_{j},j,\Delta
t,\boldsymbol{1}_{lf})$

20: end if

21: if $s^{\prime}=1$ then

22: With probability $\textrm{min}\\{1,~{}\frac{n^{\prime}}{n}\\}$, set
$\\{\boldsymbol{q}^{i},\boldsymbol{p}^{i}\\}\leftarrow\\{\boldsymbol{q}^{\prime},\boldsymbol{p}^{\prime}\\}$

23: end if

24: $n\leftarrow n+n^{\prime}$

25: $s\leftarrow
s^{\prime}\boldsymbol{1}\big{[}(\boldsymbol{q}^{+}-\boldsymbol{q}^{-})\cdot\boldsymbol{p}^{-}\geq
0\big{]}\boldsymbol{1}\big{[}(\boldsymbol{q}^{+}-\boldsymbol{q}^{-})\cdot\boldsymbol{p}^{+}\geq
0\big{]}$

26: $j\leftarrow j+1$

27: end while

28:end for

Algorithm 2 L-HNNs in NUTS with online error monitoring (build tree function)
[1]

1:function $\textrm{{BuildTree}}(\boldsymbol{q},\boldsymbol{p},u,\nu,j,\Delta
t,\boldsymbol{1}_{lf})$

2:if $j=0$ then

3: Base case taking one leapfrog step

4: $\boldsymbol{q}^{\prime},\boldsymbol{p}^{\prime}\leftarrow$ Algorithm 3
with initial conditions:
$\boldsymbol{z}(0)=\\{\boldsymbol{q},~{}\boldsymbol{p}\\}$, Steps: $1$; End
Time: $\Delta t$

5:
$\boldsymbol{1}_{lf}\leftarrow\boldsymbol{1}_{lf}~{}\textrm{{or}}~{}\boldsymbol{1}\big{[}H\big{(}\boldsymbol{q}^{\prime},~{}\boldsymbol{p}^{\prime}\big{)}+\ln{u}>\Delta_{max}^{hnn}\big{]}$

6:
$s^{\prime}\leftarrow\boldsymbol{1}\big{[}H\big{(}\boldsymbol{q}^{\prime},~{}\boldsymbol{p}^{\prime}\big{)}+\ln{u}\leq\Delta_{max}^{hnn}\big{]}$

7: if $\boldsymbol{1}_{lf}=1$ then

8: $\boldsymbol{q}^{\prime},\boldsymbol{p}^{\prime}\leftarrow$ Leapfrog
integration with initial conditions:
$\boldsymbol{z}(0)=\\{\boldsymbol{q},~{}\boldsymbol{p}\\}$, Steps: $1$; End
Time: $\Delta t$

9:
$s^{\prime}\leftarrow\boldsymbol{1}\big{[}H\big{(}\boldsymbol{q}^{\prime},~{}\boldsymbol{p}^{\prime}\big{)}+\ln{u}\leq\Delta_{max}^{lf}\big{]}$

10: end if

11:
$n^{\prime}\leftarrow\boldsymbol{1}\big{[}u\leq\exp{\\{-H\big{(}\boldsymbol{q}^{\prime},~{}\boldsymbol{p}^{\prime}\big{)}\\}}\big{]}$

12: return
$\boldsymbol{q}^{\prime},~{}\boldsymbol{p}^{\prime},\boldsymbol{q}^{\prime},~{}\boldsymbol{p}^{\prime},\boldsymbol{q}^{\prime},\boldsymbol{p}^{\prime},n^{\prime},s^{\prime},\boldsymbol{1}_{lf}$

13:else

14: Recursion to build left and right sub-trees (follows from Algorithm 3 in
[4], with $\boldsymbol{1}_{lf}$ additionally passed to and retrieved from
every BuildTree evaluation)

15: return
$\boldsymbol{q}^{-},~{}\boldsymbol{p}^{-},\boldsymbol{q}^{+},~{}\boldsymbol{p}^{+},\boldsymbol{q}^{\prime},\boldsymbol{p}^{\prime},n^{\prime},s^{\prime},\boldsymbol{1}_{lf}$

16:end if

Algorithm 3 Latent Hamiltonian neural networks evaluation in leapfrog
integration [1]

1:Hamiltonian: $H$; Initial conditions:
$\boldsymbol{z}(0)=\\{\boldsymbol{q}(0),~{}\boldsymbol{p}(0)\\}$; Dimensions:
$d$; Steps: $N$; End time: $T$

2:$\Delta t=\frac{T}{N}$

3:for $j=0:N-1$ do

4: $t=j~{}\Delta t$

5: Compute HNN output gradient $\frac{\partial
H_{\boldsymbol{\theta}}}{\partial\boldsymbol{q}(t)}$

6: for $i=1:d$ do

7: $q_{i}(t+\Delta t)=q_{i}(t)+\frac{\Delta t}{m_{i}}~{}p_{i}(t)-\frac{\Delta
t^{2}}{2m_{i}}~{}\frac{\partial H_{\boldsymbol{\theta}}}{\partial q_{i}(t)}$

8: end for

9: Compute HNN output gradient $\frac{\partial
H_{\boldsymbol{\theta}}}{\partial\boldsymbol{q}(t+\Delta t)}$

10: for $i=1:d$ do

11: $p_{i}(t+\Delta t)=p_{i}(t)-\frac{\Delta t}{2}~{}\bigg{(}\frac{\partial
H_{\boldsymbol{\theta}}}{\partial q_{i}(t)}+\frac{\partial
H_{\boldsymbol{\theta}}}{\partial q_{i}(t+\Delta t)}\bigg{)}$

12: end for

13:end for

## Checklist

  1. 1.

For all authors…

    1. (a)

Do the main claims made in the abstract and introduction accurately reflect
the paper’s contributions and scope? [Yes]

    2. (b)

Did you describe the limitations of your work? [No]

    3. (c)

Did you discuss any potential negative societal impacts of your work? [N/A]

    4. (d)

Have you read the ethics review guidelines and ensured that your paper
conforms to them? [Yes]

  2. 2.

If you are including theoretical results…

    1. (a)

Did you state the full set of assumptions of all theoretical results? [N/A]

    2. (b)

Did you include complete proofs of all theoretical results? [N/A]

  3. 3.

If you ran experiments…

    1. (a)

Did you include the code, data, and instructions needed to reproduce the main
experimental results (either in the supplemental material or as a URL)? [N/A]

    2. (b)

Did you specify all the training details (e.g., data splits, hyperparameters,
how they were chosen)? [N/A]

    3. (c)

Did you report error bars (e.g., with respect to the random seed after running
experiments multiple times)? [N/A]

    4. (d)

Did you include the total amount of compute and the type of resources used
(e.g., type of GPUs, internal cluster, or cloud provider)? [N/A]

  4. 4.

If you are using existing assets (e.g., code, data, models) or
curating/releasing new assets…

    1. (a)

If your work uses existing assets, did you cite the creators? [N/A]

    2. (b)

Did you mention the license of the assets? [N/A]

    3. (c)

Did you include any new assets either in the supplemental material or as a
URL? [N/A]

    4. (d)

Did you discuss whether and how consent was obtained from people whose data
you’re using/curating? [N/A]

    5. (e)

Did you discuss whether the data you are using/curating contains personally
identifiable information or offensive content? [N/A]

  5. 5.

If you used crowdsourcing or conducted research with human subjects…

    1. (a)

Did you include the full text of instructions given to participants and
screenshots, if applicable? [N/A]

    2. (b)

Did you describe any potential participant risks, with links to Institutional
Review Board (IRB) approvals, if applicable? [N/A]

    3. (c)

Did you include the estimated hourly wage paid to participants and the total
amount spent on participant compensation? [N/A]


# Nambu-Goto Dynamics of Field Theory Cosmic String Loops

Jose J. Blanco-Pillado  Daniel Jiménez-Aguilar  Joanes Lizarraga  Asier Lopez-
Eiguren  Ken D. Olum  Ander Urio  and Jon Urrestilla

###### Abstract

We perform a detailed comparison of the dynamics of cosmic string loops
obtained in cosmological field theory simulations with their expected motion
according to the Nambu-Goto action. We demonstrate that these loops follow the
trajectories predicted within the NG effective theory except in regions of
high curvature where energy is emitted from the loop in the form of massive
radiation. This energy loss continues for all the loops studied in this
simulation until they self-intersect or become small enough that they
annihilate and disappear well before they complete a single oscillation. We
comment on the relevance of this investigation to the interpretation of the
results from cosmological field theory simulations as well as their
extrapolation to a cosmological context.

##  1 Introduction

Many non-linear field theories have in their spectrum solitonic solutions that
play a significant role in the dynamics of the theory. These objects have been
extensively studied over the past few decades in connection to many different
fields of research from condensed matter to cosmology. One important aspect of
these studies has been the investigation of the dynamics and interactions of
these solitons.

One can explore these issues by performing numerical experiments in a lattice
simulation. This lattice does not know anything about the solitonic
configuration and simply evolves every point in space as any other one
following the equations of motion. One can then recover the information of the
evolution of the fields in the lattice and interpret the results in terms of a
collection of solitons in motion possibly interacting with each other. This is
a useful way to analyse the microphysics of the solitons since one can probe
the profiles of the fields at scales smaller than the characteristic soliton
size111This is always the case since one of the requirements of our lattice
spacing should be that it is always smaller than the soliton characteristic
scale.. In particular, one could use this technique to find out whether the
solitons remain in their lowest energy state or whether they are typically
excited at some point during the simulation.

On the other hand, this procedure is computationally very expensive if one
wants to simulate a volume that is large compared to the size of the soliton.
This is particularly relevant in simulations that involve many solitons or
whenever one is interested in investigating an effect whose typical scale is
much larger than the soliton’s thickness. In these cases one is forced to look
for some effective theory that captures the degrees of freedom that are
relevant for the problem without having to simulate every point in a lattice.
This drastic reduction in the number of degrees of freedom that one needs to
compute to simulate the soliton’s dynamics suggests the possibility of another
kind of simulation. Such simulations based on effective theories allow for a
much larger dynamic range in the simulation, which in turn will help us obtain
a better understanding of the large scale dynamics of the problem. However,
one needs to make sure that there are no microphysical effects that are
missing in the effective theory that could potentially become relevant for the
large scale dynamics that one wants to faithfully reproduce in the simulation.

These two types of simulations are therefore complementary. One can use the
lattice simulations to learn the important field theory effects that need to
be accounted for in the effective theory of the solitons. Once this is done,
one should be able to find some common ground where both these simulations can
be compared and where an agreement can be reached on the important dynamics to
study. Once this is achieved, an extrapolation to the interesting scales can
be safely done using the effective theory.

In this paper we would like to take the first step towards showing this
agreement between these two approaches in the context of local cosmic string
networks. In this case, we will consider the Abelian-Higgs model as the field
theory where local cosmic strings occur as solitons [1] and the Nambu-Goto
(NG) action as the effective theory at low energies [2, 3]. Cosmological
simulations of both types, lattice field theory [4, 5, 6, 7, 8, 9, 10, 11] and
Nambu-Goto dynamics [12, 13, 14, 15, 16, 17, 18, 19, 20, 21], have been
extensively studied in the past. However, there seems to be an important
disagreement about the abundance of non-self-intersecting (NSI) loops between
these two numerical techniques. As its name suggests, non-self-intersecting
loops are loops that in their dynamics do not self-intersect. All the loops
found so far in field theory simulations do self intersect, and the loops seem
to decay in a short period of time. Nambu-Goto simulations also produce self-
intersecting loops, but over the course of their evolution a large number of
non-self-intersecting loops appear. These non-self-intersecting loops are
important, since their main energy loss mechanism is via gravitational waves,
whereas loops that continually intersect lose energy also via massive
radiation.

This apparent disagreement has been the source of a debate for a very long
time, since some of the first field theory simulations of string networks [4].
One possible explanation could be that field theory simulations have not been
lucky enough to produce a large enough NSI loop. This is a reasonable
explanation, since field theory simulations have much fewer loops in general,
and Nambu-Goto dynamics shows that NSI conditions in a random loop are rare
compared to self-intersecting ones.

It is important to remember at this point that the NG action is only an
approximation to the actual dynamics of the strings and, as we will describe
in detail in this paper, can certainly break down under some special
circumstances. However, field theory simulations of individual smooth strings
have shown conclusively that these strings follow almost exactly the Nambu-
Goto dynamics [22, 23]. This makes this disagreement with field theory
simulations more puzzling and has prompted some authors to suggest that
another reason for the disagreement may be due to the presence of excitations
on the strings in the network simulations [8, 24].

These field theory excitations have been found in many soliton solutions, in
particular in cosmic strings [25, 26, 27, 28], and have been recently studied
in a series of papers in different models in [29, 30, 31]. The results of
these studies show that many solitonic solutions may store a significant
amount of energy in the form of excitations [29, 32]. In particular, the phase
transition that creates these solitons could lead to an excited initial state
due to the extra energy floating around the bulk as the soliton is formed.
Furthermore, these papers show that the time scale for the decay of these
excitations could be much larger than the typical time scale of the soliton,
the light-crossing time of the width of the soliton [33, 29]. This suggests
that the effect of these excitations could play some role in the dynamics of
strings in field theory simulations. In particular, the presence of these
extra modes localized on the string could modify their equation of state
rendering their dynamics quite different from the one predicted by the Nambu-
Goto action. However, the existence of this extra energy on cosmic strings has
never been conclusively shown in any field theory simulation of the string
network222The closest study of this kind in the literature is the one
presented in [32] where the authors perform cosmological simulations of $2+1$
dimensional global vortices. Their results indicate that the amount of energy
in this case is rather low..

In this paper we will study this problem by carefully analysing the evolution
of loops extracted directly from field theory simulations. Specifically, we
obtain the position and velocity of some of the largest loops found in the
course of a field theory simulation of a network of strings and compare their
evolution with the one predicted by the Nambu-Goto action.

The results indicate that these loops follow the same trajectories as their NG
counterparts, except in localized regions where the curvature of the strings
is large compared with the string core thickness, where the NG approximation
is not good by definition. As we will argue in the main part of the text, this
shows that there is no significant deviation from the Nambu-Goto dynamics due
to a new equation of state for the strings in the parts where the curvature is
not high. In other words, it seems that the strings in our simulations do not
have a large amount of energy stored in them in the form of localized
excitations.

The organization of the paper is the following. In section 2, we comment on
the characteristics of the field theory simulations that were used to generate
the field theory trajectories of the loops. In section 3, we describe the
techniques we use to compare the dynamics of field theory loops and their
Nambu-Goto predictions. In section 4, we show our results with a few snapshots
of the string trajectories comparing both field theory and Nambu-Goto.
Finally, in section 5, we comment on the implications of these results for the
cosmological extrapolation of field theory cosmic string networks.

##  2 Field theory simulations of cosmic string loops

The field theory that we will investigate in this paper is the Abelian-Higgs
model, whose Lagrangian density,

$L=D_{\mu}\phi
D^{\mu}\phi^{*}-\frac{\lambda}{4}\left(|\phi|^{2}-\eta^{2}\right)^{2}-\frac{1}{4e^{2}}F^{\mu\nu}F_{\mu\nu}~{},$
(2.1)

describes the dynamics of a complex scalar field, $\phi(x)$, coupled to a
vector field, $A_{\mu}(x)$, through the covariant derivative,
$D_{\mu}\phi=(\partial_{\mu}-iA_{\mu})\phi$. Furthermore, the usual field
strength for the vector field is given by
$F_{\mu\nu}=\partial_{\mu}A_{\nu}-\partial_{\nu}A_{\mu}$. We will consider the
case $\beta=\lambda/(2e^{2})=1$, which means that the masses of the
excitations in the vacuum for the vector and scalar fields are equal:
$m=m_{s}=m_{v}$.

It is well known that the equations of motion obtained from this theory allow
for the existence of solitonic field theory vortices [1]. In $3+1$ dimensions
the vortices become strings whose energy is concentrated in a core thickness
of the order of $\delta\sim m^{-1}$.

In order to investigate the loops in this theory we follow the prescription
detailed in [24], where network loops are created from some random initial
conditions in the lattice. We point the interested reader to [24] for details,
but we summarize here the necessary basic information. After discretization of
the Hamiltonian that corresponds to the Lagrangian (2.1), we obtain the
equations of motion and solve them in cubic lattices with periodic boundary
conditions.

The initial configuration of the system is chosen to be such that all fields
are set to zero except for the scalar field, which is set to be a stationary
Gaussian random field with a power spectrum given by

$P_{\phi}=Ae^{-kl_{\phi}}\;,$ (2.2)

where the amplitude $A$ is chosen so that
$\left\langle\left|\phi\right|^{2}\right\rangle=\eta^{2}$. The free parameter
which rules the initial randomness is the correlation length $l_{\phi}$, which
can be set initially to different values.

These random initial conditions lead to a considerable excess of energy in the
simulation volume. Therefore, a cooling process is applied using a diffusive
period of evolution, and once a smooth field distribution is obtained, the
network evolves following the true equations of motion in flat space.

It is important to mention that after the diffusive period the network is at
rest, which means that any loop at this initial stage will start from a static
configuration. We will comment on these primordial loops in Appendix A, but
let us for now comment on the fact that they are not the main interest to us
since they are not representative of the typical loops in a network
simulation. In the course of the evolution of the network, loops of strings
will be formed by either self-intersections or intercommutation of long
strings and these are the loops that we are mainly interested in. By this
stage the string is in motion and so the loops formed are not static. Note
that the numerical evolution of the string network was performed in flat
spacetime, in other words on a Minkowski lattice spacetime.

In order to localize the position of the strings and the subsequent loops we
first identify all the plaquettes with a non-trivial winding in our
simulations. In this way the connection between all the centers of these
plaquettes will constitute the string. We follow the evolution of the network
outputting the windings and confirm the formation of loops coming from
intersections by visual inspection.

Since, as mentioned, the simulations are done using periodic boundary
conditions, all the strings in the box can be considered to be closed loops;
but these loops can be “broken” by the periodic boundary conditions. Thus, if
one were to plot the loop directly, some would not appear to be a connected
piece of string. In order to avoid this we reconstruct the loops that are
broken by the periodic boundary conditions applying spatial translations and
assigning new coordinates to the string positions. In this manner, we have a
list of connected positions in space for all loops.

In total we have analyzed 7 loops (and their descendants) from the simulations
in [24]. These loops were obtained using two different correlation lengths,
$l_{\phi}=15$ and $25$ in $\eta^{-1}$ units. All of them were produced using
lattices of $N=1024$ points per dimension with a spatial resolution of $\delta
x=0.125$ and temporal resolution of $\delta t=0.2\delta x$, again in
$\eta^{-1}$ units. We refer the reader to [24] for more specific details on
the preparation of these simulations. Moreover, the loops, which are
reconstructed following the above prescription, are output at each time step
of the evolution so that the field theory information available for the Nambu-
Goto prediction is the most accurate possible. Furthermore, as explained in
[24] by the time the loops we use here get formed, the large scale dynamics of
the string network is consistent with a scaling regime 333Recall, however,
that the simulation here is done in flat space. It would be very interesting
to analyze loops from a lattice simulation on an expanding background. Some
effort in this direction is already underway..

This data, extracted directly from the field theory simulations, is the
starting point of our analysis.

##  3 Comparison of Field Theory data with the Nambu-Goto action

As we said above, the effective action that describes the dynamics of local
strings is expected to be the Nambu-Goto action [2, 3]. This can be justified
by making a judicious choice of coordinate system around the center of the
string and integrating the action along the transverse directions of the
string [34]: the resulting action can be shown to be of the Nambu-Goto form.
This argument rests on several assumptions that we now list in detail.

First, it assumes that the local curvature of the string is small compared to
its thickness. In fact, one can consider the Nambu-Goto action as the lowest
order approximation of an infinite expansion in terms of the ratio of
$\delta/R$, where $\delta$ is the nominal thickness of the string soliton and
$R$ is the radius of curvature of the string in space. In a cosmological
setting, truncating this series keeping only the first term seems quite
reasonable since the separation of scales from the microphysical size of the
strings to any cosmologically relevant scale is phenomenally large. Of course,
this separation of scales is not so large in a field theory simulation of a
string network.

The second assumption, which is somewhat related to the previous one, is that
the string does not lose energy by radiation in the course of its evolution.
This is of course built in the NG action since, as we discuss below, there is
a conservation law for the invariant energy of a loop. However, from the point
of view of field theory, one might imagine a situation where solitonic strings
lose part of their energy into radiation in the form of propagating modes in
the bulk. Of course this cannot happen for a relaxed static string, since by
definition this object is the lowest energy configuration with the particular
boundary conditions, the winding of the scalar field. Boosting this object can
not lead to radiation either. So the only way this string can radiate is due
to acceleration. This is easy to achieve in strings since during their
evolution they can develop regions of curvature that will induce acceleration.
The question is then a quantitative one. How much energy is radiated from the
typical acceleration present in the evolution of strings? In order to answer
this question we should remind ourselves that all the propagating modes in
this model of local strings are massive. This suggests that one should wiggle
the string with a frequency at least of the order of this mass, $m$ in our
case, if one wants to produce any radiation. Below this frequency, the source
for the radiation does not have a large enough frequency to produce
propagating modes. This argument has been extensively used in the past and
demonstrated explicitly in numerical simulations of strings in [23] and more
recently in the analogous situation in domain wall strings in [30]444In fact,
the situation is a little more complicated since the thickness of the source
(the oscillating string) is typically also of the order of the inverse of the
mass of the radiated particle, namely, $\delta\sim m^{-1}$. This means that
for higher frequencies the radiation is also cut-off due to interference
effects [30]..

There are of course moments where the string can release part of its energy.
The simplest way to visualize this is in the lower dimensional process of
vortex-antivortex annihilation. In these events, it is clear that the
arguments leading to the conclusions above do not apply since the topological
stability of the solitons disappears. Extending this to $3+1$ dimensions, we
can classify other instances where similar processes occur. A clear example of
these kinds of events is the interaction between two long strings, the so-
called intercommutation process by which string loops can detach themselves
from long strings [35]. Of course, this process cannot be described by the NG
action. Other processes closely related to these are cusp formation and kink-
kink collisions. During the formation of the cusp, part of the string
annihilates with itself releasing energy in the process [22]. Similarly, kink-
kink collisions [36] or in general the appearance of very high curvature
regions [23, 30] lead to a similar energy ejection from the string. All these
non-perturbative processes are by now well understood and need to be accounted
for separately from the NG evolution.

Finally, another important assumption in the use of the NG action to describe
strings is the general expectation that, in their rest frame, the solitonic
strings would be well approximated by the static solution of lowest energy.
The underlying idea for this expectation is the supposition that excitations
on the string will decay in a time scale of the order of the light-crossing
time of the thickness of the string. This is a very small time scale to have
any relevance for the evolution of a loop even in a field theory simulation.
However, this assumption has been recently brought into question [24] due to
the existence of localized excitations of the field theory string that can
have a long lifetime [33, 29, 32]. The presence of this extra energy can
change the equation of state of the strings and so modify the trajectory of
strings. These ideas have been recently explored in several papers in lower
dimensional models with solitons [29, 32, 30]. The results of these studies
seem to indicate that even though some of these solitons could have some extra
energy at the moment of formation, it is difficult to see how they can achieve
the necessary significance to alter the evolution of the solitons.

Here we set out to investigate the relevance of all these possible effects in
the evolution of cosmic string loops from cosmological network simulations. In
particular, we will compare the evolution of the string extracted from field
theory following the procedure we indicated earlier with the one that one
would infer from NG dynamics. In the following section we will explain how to
obtain the prediction of the NG dynamics from the field theory data at any
moment in time.

###  3.1 The Nambu-Goto dynamics for a cosmic string loop

The Nambu-Goto action for a relativistic string is given by

$S_{NG}=-\mu\int{d^{2}\xi\sqrt{-\gamma}}~{},$ (3.1)

where $\mu$ is the energy per unit length of the string and $\gamma$ is the
induced metric on the worldsheet parametrized by the coordinates $\xi^{1,2}$.
The equations of motion for a string propagating in flat spacetime can be
obtained from this action and are given by555See, for example, [37] for an
account of all the details of the NG action.

${\bf\ddot{x}}-{\bf x^{\prime\prime}}=0~{},$ (3.2)

where ${\bf x}(t,\sigma)$ parametrizes the position of the string and dotted
and primed quantities denote their differentiation with respect to the two
worldsheet parameters, $(t,\sigma)$. Moreover, we have also imposed the gauge
conditions

${\bf\dot{x}}\cdot{\bf x^{\prime}}=0~{},$ (3.3)

so the only physical velocity is perpendicular to the string, and

${\bf\dot{x}}^{2}+{\bf x^{\prime}}^{2}=1~{},$ (3.4)

which means that the spacelike parameter, $\sigma$, is proportional to the
energy per unit length along the string. These equations can be integrated, so
the most general solution is of the form

${\bf x}(t,\sigma)=\frac{1}{2}\left({\bf a}(\sigma-t)+{\bf
b}(\sigma+t)\right)~{},$ (3.5)

where the constraints impose the conditions

$\displaystyle|{\bf a^{\prime}}|$ $\displaystyle=$ $\displaystyle 1~{},$ (3.6)
$\displaystyle|{\bf b^{\prime}}|$ $\displaystyle=$ $\displaystyle 1~{}.$ (3.7)

This means that all one needs to do in practice to obtain the evolution of a
string is to find the form of the two functions ${\bf a}(\sigma_{-})$ and
${\bf b}(\sigma_{+})$ with respect to their argument. In the following we will
describe an algorithm to obtain these functions from the position of the
string at two different time steps. This will allow us to apply this procedure
to the data obtained from the field theory simulation.

Note that after obtaining these functions at a particular moment, eq. (3.5)
will allow us to find the position of the string at any subsequent time. This
can in turn be compared with the position in field theory. If the string moves
exactly as the NG predicts, we could obtain the form of these functions at any
moment in time and the result would be the same. In the following we will
explain that it will be convenient to repeat this procedure at several times.

Another important point about the NG dynamics is the fact that one can find
the integral of the parameter $\sigma$ along the string. This is a constant of
motion of any loop of string and so can also be used as a measure of the NG
dynamics.

###  3.2 Obtaining the NG dynamics from field theory data

As we mentioned earlier, we extract the information about individual loops
from the lattice field theory simulation of the network by identifying their
position at any moment in time. Here we will explain in detail how to
transform that into a prediction of the NG evolution.

We start with the position of the string through the lattice as a list of
spatial positions of plaquette centers, ${\bf p}_{n}$, $n=1,2,\ldots$. The
first thing we do is to smooth out these position vectors since otherwise we
will have big jumps related to the discrete nature of the lattice. One
reasonable possibility is to smooth this data by a Gaussian window function
whose width is given by a few lattice spacings, $\text{M}\delta x$. In order
to justify this choice, let us first remember that in order to have a faithful
simulation of the relevant dynamics the thickness of the string is somewhat
larger than the lattice spacing. Therefore, one should not expect to know the
position of the center of the strings with a precision much larger than this
width666We use $\text{M}=2$ for the loops presented in this paper. We have
checked that our results do not change significantly using $\text{M}=5$ ..

Using this smoothed data, we obtain the list of vectors that describe the
normalized tangent vector of the string by computing

${\bf\hat{p}^{\prime}_{n}}=\frac{{\bf p}_{n+1}-{\bf p}_{n}}{|{\bf
p}_{n+1}-{\bf p}_{n}|}~{}.$ (3.8)

Now we need to compute the velocity vectors for each point of the string. In
order to do that, we will assume the NG evolution and consider that the string
moves in the direction perpendicular to its tangent vector. This allows us to
compute this velocity using the following algorithm.

First, we find the point of the string at a later time described by $t+\Delta
T$ that is the intersection of the plane perpendicular to the tangent vector
at the original position of the string (at time $t$) with the string at
$t+\Delta T$. Let’s denote this point by $\tilde{\bf p}_{n}(\Delta T)$. Then
we estimate the velocity vector of this segment of the string in the NG
approximation by taking

${\bf\dot{x}}_{n}=(\tilde{\bf p}_{n}(\Delta T)-{\bf p}_{n})/\Delta T~{}.\\\ $
(3.9)

Here we should comment on an important point. In the previous algorithm we
have not specified what is the relation between the time interval between the
two sets of strings positions we have in the original data and the computing
time interval in the numerical field theory simulation, $\delta t$. In our
case, we have experimented with different values and finally settled on
$\Delta T=8~{}\delta t$. The reason to take these two snapshots of the string
separated in time larger than the minimal possible time separation, $\delta
t$, is to try to smooth out possible errors in the estimate of the
velocity777Using a smaller value of $\Delta T$ induces the presence of points
of super-luminal motion. This is of course an error induced in regions of high
velocity of the string. Note also that regions of self-annihilation would give
rise to these problems. However, it is clear that this is to be expected since
in those regions the string does not behave as NG, so this algorithm should
definitely fail..

Using this velocity and the normalized tangent vectors we can now compute the
tangent vector correctly parametrized according to the usual NG gauge. This
means that we can define the new tangent vectors as

${\bf
x}^{\prime}_{n}=\left(\sqrt{1-|{\bf\dot{x}}_{n}|^{2}}\right){\bf\hat{p}^{\prime}_{n}}~{}.$
(3.10)

Using the velocity and the tangent vectors we can in turn compute the
functions ${\bf a}^{\prime}$ and ${\bf b}^{\prime}$ using

$\displaystyle{\bf a}^{\prime}_{n}={\bf x}^{\prime}_{n}-{\bf\dot{x}}_{n}\,,$
(3.11) $\displaystyle{\bf b}^{\prime}_{n}={\bf
x}^{\prime}_{n}+{\bf\dot{x}}_{n}\,,$ (3.12)

and from here it is easy to find the position of the string at any moment in
time following the prescription of the NG solution in eq. (3.5).

Finally, using this data we can easily compute the local Lorentz factor
associated with each segment, namely,

$\Gamma_{n}=\frac{1}{\sqrt{1-|{\bf\dot{x}}_{n}|^{2}}}~{},$ (3.13)

as well as the amount of $\sigma$ parameter in each segment of the string by
computing

$\Delta\sigma_{n}=\Gamma_{n}|{\bf x}_{n+1}-{\bf x}_{n}|~{}.$ (3.14)

Integrating this over our list of elements of the string we obtain the
invariant energy of the loop, $E_{NG}=\mu\sum_{n}{\Delta\sigma_{n}}$.

As we mentioned earlier, in the course of the reconstruction there are points
that lead to an estimation of the velocity from $|{\bf\dot{x}}_{n}|$ very
close to the speed of light or even above it. In order to suppress the
pathological behaviour that these points could have on the total energy we
have decided to put an artificial cap to the velocity of each individual
point. We replace any estimate of $|{\bf\dot{x}}_{n}|>0.9$ by
$v_{\text{max}}=0.9$. We have tried other regularization procedures and
checked that the total energy of the loop is not significantly affected by the
different procedures.

##  4 Results

Using the techniques we have outlined above, we can compare the evolution of
field theory loops from our simulations with the motion predicted from the NG
dynamics. Using the data from two steps in the field theory simulation
separated by $\Delta T$ since the formation of the loop, we build the ${\bf
a}(\sigma)$ and ${\bf b}(\sigma)$ functions. This allows us to plot the
predicted string NG position for all times. This comparison shows that the
field theory string follows the NG solution very accurately for the large
majority of the string length. This seems to suggest that these loops
extracted directly from the simulation are not endowed with a significant
amount of extra energy as conjectured in [24], or at least, not enough to
change the trajectory of the string perceptibly. A large amount of energy in
bound states would change the equation of state of the string and its local
velocity would be changed with respect to the one obtained in NG. This is not
observed for most parts of the string length.

There are, however, regions where we see a departure of the field theory
string position when compared to NG. Many of the places where we see this
departure are regions where the NG dynamics predicts the existence of a high
curvature section on the string. As we mentioned earlier it is therefore not a
surprise that the field theory string does not follow the NG prediction in
those regions.

An example of such local departure from NG dynamics is shown in figure 1.

Figure 1: Several snapshots of the evolution of one of the field theory loops.
We show in blue the position of the string obtained directly from the lattice
simulation. In red is the predicted position obtained from the reconstruction
of the NG data at the initial conditions and evolved using the NG dynamics
until the time shown. The agreement between these two descriptions is very
good for most of the loop’s evolution. We have zoomed in on a region of the
string at a particular moment where there is a visible departure between them.

We present several snapshots of this loop’s evolution in field theory (in
blue) and in the NG dynamics (in red) obtained using the algorithm described
in the previous section. In order to represent the field theory string, we
give it a width of the order of $\delta$ (the thickness of the solitonic
object). The predicted position of the NG action is hidden inside of the blue
tube describing the position of the field theory string for most of the
string. The two curves only deviate from one another in a small section of the
whole string. In that region, the NG string curves itself at the scale of the
order of the string thickness, but the field theory string does not do that
and finds a shortcut.

These episodes of high curvature act as a source of energy loss from the
string. Some of these events resemble the cusp annihilation simulated several
years ago in [22]. Others just correspond to the interaction of wiggles on the
string that produce high curvature regions. In some cases, these interactions
lead to the formation of tiny daughter loops that immediately annihilate in
the field theory side.

The subsequent evolution of the field theory string does not follow the NG
prediction after those episodes. The reason for this is also clear. The NG
action conserves energy, and therefore it does not account for this energy
loss mechanism. This means that the evolution from field theory would start
being different in the region where energy is radiated. As time passes, this
departure from FT and NG spreads over the rest of the string. If one waits
long enough, the difference becomes quite visible, and if one were to continue
the comparison forward, the shapes of the loops would grow more different.
However, this is not a real measure of the different local dynamics.

In order to do a better job of identifying the reason for the different
evolution we follow a procedure also used in the past in [22]. After we have
identified one of these high curvature events on the string, we reconstruct
the NG data again. This yields different ${\bf a}$ and ${\bf b}$ functions
that should be valid for the subsequent evolution. The interesting point is
that the field theory after these events have passed is again accurately
described by the new NG data. We show in figure 2

Figure 2: The evolution of the field theory string seems to follow the NG
prediction obtained from the reconstruction of the string after the high
curvature event. There is no visible departure between the NG and the FT
descriptions.

a comparison of the position of the string obtained from field theory to the
NG reconstruction after the first episode of high curvature radiation. We
notice that there is no visible departure of the field theory evolution from
the prediction of NG.

This behaviour continues for a while until a new episode occurs. This is shown
in figure 3,

Figure 3: Another high curvature event on the same loop. The NG description
(in red) is the one obtained after the first event in figure 1.

where we clearly see another region where the NG prediction deviates from the
field theory result in a localized region. It is clear that the string in
field theory does not want to curve itself so much as the NG predicts and
takes a shortcut. This process radiates again some portion of the energy of
the string.

We have seen a similar behaviour in all our loops. Some of the examples are
clear, but some other ones are harder to visualize since more than one of
these high curvature events happen to have some non-trivial overlap in time
888However, we would like to emphasize that the results we present here, with
this particular loop, are indeed a good representation of what we obtain in
our analysis of all the other $6$ loops from [24]. . In fact, this already
happens in our example loop. We show in figure 4

Figure 4: Another high curvature event on the same loop. The NG description
(in red) is the one obtained after the first event in figure 1. The evolution
between these episodes is well represented by the NG dynamics.

a third event situated quite far away in space from the previous one but that
overlaps in time with the event represented in figure 3.

We also look at the energy computed from the local reconstruction of the NG
string obtained from the field theory data at any moment in time. In a purely
NG dynamics, this quantity (the total amount of $\sigma$ of the loop) should
be a constant of motion. We plot in figure 5 this energy for the same loop
that we discussed before. We observe that the energy overall tends to go down.
There are some episodic events where the decrease in energy is sharper, and
some of those can be linked with the high-curvature events. We mark in figure
5

Figure 5:  Total amount of invariant energy (total amount of $\sigma$) for the
NG reconstruction of the field theory data of the loop analyzed in other
figures. The purple regions correspond to the high curvature events presented
in figures 1,3 and 4. The grey region corresponds to figure 2, where the
energy remains roughly constant. The white region (between $t\sim 400$ and
$550$) corresponds to figure 6, where there is a small deviation from the NG
dynamics.

the three different episodes that we have been discussing above by shading in
light purple the ranges of times displayed in figures 1, 3, and 4. Looking at
the energy, there seems to be a connection between these events and the
regions where the energy starts to decrease.

There are other instances (between $t\sim 200$ and $t\sim 400$) where the
energy seems to be constant, represented in grey in figure 5. These correspond
to the times where the NG reconstruction was made after the first episode, as
shown in figure 2, where there is a very good agreement with the FT data.

Figure 6: The evolution of the field theory string in times where the energy
decreases. One can not identify any single event of high curvature associated
with this decline in energy. Instead we see several regions of high curvature
where the deviation from NG is clearly visible but not as striking as in the
previous examples.

There are, however, regions where the energy slowly decreases that are not so
obviously associated with any of these individual events (see figure 6, which
corresponds to times 450, 500 and 550). The reasons for these are not so
clear. The curvature does not seem to be very high, and could also be
understood as small scale structure in the loop that leads to energy loss by
radiation. In other words, this could be due to events similar to the ones
already mentioned that are not so clearly visible in our procedure. Looking at
figure 6 one indeed sees several regions of slight departure between NG and FT
dynamics, especially at $t=550$. Note that identifying small regions of high
curvature from the NG reconstruction is sometimes quite difficult due to
numerical error. In the future we will design new methods to quantify this
effect and understand the reason for this energy loss in these regions as
well. It is remarkable, and somewhat puzzling, that though the energy goes
down roughly by the same percentage as in a high-curvature event, the visual
inspection of the loop dynamics does show a rather small deviation from the NG
trajectory. This could just be due to the combined effect of several smaller
regions where the deviation is small instead of a single large event like in
the other cases.

##  5 Conclusions

In this paper we have compared the evolution of field theory loops obtained in
the course of a cosmological lattice simulation with their expected dynamics
in the NG approximation. Understanding the discrepancy between these two
approaches is of paramount importance in order to make an accurate prediction
of the observational signatures of strings. In particular, it is crucial in
the estimate of the gravitational wave signature from strings in current and
future gravitational wave observatories (see [38] and references therein).

Our investigations show that loops in field theory seem to behave according to
the NG action in regions where the curvature is not high. The visual
comparison with the NG motion does not support the need for any departure in
the equation of state of field theory loops, also, in regions where the
curvature is not high. The strings move with the local trajectories dictated
by NG. However, we have found that the strings lose part of their energy in
the course of their evolution. Some of these energy loss regions correspond to
high curvature sections where some portion of the string annihilates with an
adjacent part. These high curvature events are indeed predicted by the NG
evolution obtained from the original reconstruction of the field theory
string. Of course, the NG dynamics cannot account accurately for the
subsequent evolution of the field theory string since part of the energy of
the string is lost in these events. This can be bypassed by reconstructing the
NG data again after one of these incidents. The result we obtain by following
this prescription seems to show that the evolution of the strings is again
described by NG with this new data. There are other instances where the string
loses energy, which cannot be so clearly pinpointed as regions of high
curvature, at least following our visual inspection. Remarkably, though, the
trajectory of the string does not seem to be altered perceptibly in these
events. They deserve further investigation.

The picture that emerges from our detailed comparison of both descriptions of
the string motion is the following. Most of the time, loops behave as NG
predicts, but there are instances where the NG action breaks down and one
needs to interrupt this comparison for a while until the NG behaviour resumes
again. The study of the conserved NG energy backs up that there are instances
where the energy drops that correspond to high curvature events.

This localized energy loss mechanism makes the loops shrink and sometimes
self-intersect before they have a chance to oscillate for a full period. So at
the end the resultant loops are too small to expect them to behave as NG and
they finally disappear. This could explain why we do not get at the end any
non-self-intersecting loop from our simulation even though the dynamics of
loops is well explained by NG for most of their evolution.

The question arises then: if the loops behave almost everywhere like NG, then,
would one expect to get NSI loops also in FT, and thus a big chunk of the
energy of the network be released as GW? The direct obvious way of answering
this question might be to keep simulating loops of this kind until a NSI loop
is found in field theory. This is not a good strategy, because, as indicated
previously, NSI conditions are not so easy to come by. Many large NSI loops
are found in NG simulations [20], because one simulates a much, much larger
volume with many more loops [20]. Unfortunately, we do not have the dynamic
range in field theory to do such simulations. Of course, we may be lucky and
find one such loop in our simulations after a large number of them.

Another idea would be to start with a different set of loops. For example, we
could get loops from field theory simulations in the radiation or matter era.
These loops should be smoother and have a greater chance to become non-self-
intersecting.

One would also be tempted to look for larger initial loops. However, even
though it is important to have large loops, so their size is large compared to
their thickness, what we have seen in these simulations is that this is not
the most important fact. One can have a very large loop with wiggles that lead
to high curvature regions, which would thus lose energy by this mechanism.

The best scenario would be to start with a large enough loop that radiates
most of this energy in the high curvature regions in its first few moments
leaving behind a smoother loop that now should behave mostly as NG (except
maybe for the presence of cusps). This expectation is based on the results
obtained in field theory simulations of the collisions of wiggles that lead to
high curvature regions [23, 30]. That is the analogous situation to what we
are seeing here in these loops, but in long, infinite strings. The results
there ndicate that the radiation from these events decreases quite fast after
their first encounter. The wiggles become milder, and their subsequent
interaction is not so violent. This argues for a period of smoothing of the
loops of the order of one oscillation time after which the loops would become
quite close to NG in all their evolution.

If this is the correct view, it means that loops created in a real
cosmological network have a transient period where they emit massive radiation
from these highly curved regions. After this initial stage (of the order of
the period of the loop) this effect should smooth out the loop and loops would
behave as the NG action predicts. Furthermore, as the universe evolves, if the
above picture is right, the sizes of structures on loops would increase
proportionally to the horizon distance, while the string core size remains
fixed. Thus, over cosmological time, the curvature radii seen on loops would
become many orders of magnitude larger than the string thickness, the
radiative processes we see here would disappear, and the loop motion would be
accurately given by NG.

Nevertheless, the study reported here is quite preliminary. We have analyzed
only a few loops, and we do not understand their dynamics completely. So it is
possible that more is going on than the simple description above. In that case
many alternative scenarios [4, 7, 8, 24] may be possible.

In summary, we believe we have made a major step forward in our understanding
of the dynamics of loops from field theory simulations. It is clear that loops
appear to move as NG for most of their evolution. One clear situation in which
this is not happening is, not surprisingly, in events of high-curvature. There
are other instances in which the energy of the loops does not seem to follow a
NG prediction, and yet the behaviour of the loop does not present big
departures from the NG trajectory. The details of this picture need to be
confirmed by further study. In the future we will perform new numerical
experiments and use different field theory simulations that could corroborate
this picture. Some of this work is already underway.

## Acknowledgments

We are grateful to Mark Hindmarsh for discussion and for letting us use the
data from his work. We also thank Jose R.C.C.C. Correia, Carlos J. A. P.
Martins, Jose M. Queiruga, Tanmay Vachaspati and Alex Vilenkin for useful
discussions and collaboration on related work. This study is supported in part
by the PID2021-123703NB-C21 grant funded by MCIN/ AEI /10.13039/501100011033/
and by ERDF; ”A way of making Europe”, the Basque Government grant
(IT-1628-22), the Basque Foundation for Science (IKERBASQUE) and the National
Science Foundation under grant number 2207267\. Part of this work has been
performed under the Project HPC- EUROPA3 (INFRAIA-2016-1-730897), with the
support of the EC Research Innovation Action under the H2020 Programme; in
particular, AU gratefully acknowledges the support of the Department of
Physics of the University of Helsinki and the computer resources and technical
support provided by CSC Finland, as well as support from the University of the
Basque Country grant (PIF20/151).

##  Appendix A Primordial loops

As we explained in the main part of the text, we create our string network
with an initial period of diffusion. One can of course look at the initial
evolution of some of the loops created at this time as well. However, we
should note that they are quite different from the ones we have analyzed in
the rest of the paper. First of all, they are much smoother due to the period
of diffusion and furthermore they are all created at rest, meaning all the
segments of the string start their evolution without any initial velocity.
These properties make these loops rather special from the point of view of
their Nambu-Goto evolution. It is easy to show that an initially static loop
will overlap with itself along the entirety of its physical length in a half
of its period if it moves according to the NG description [39]. It is
therefore clear that we cannot use these loops to illustrate the typical
behaviour of a loop in a realistic cosmological setting. However, as we will
describe in the following, we can use these primordial loops to check the
validity of our results and our conclusions.

We show in figure 7

Figure 7: Snapshots of the evolution of a primordial loop. The loop starts
from rest. We notice how the evolution of the field theory (in blue) is very
close to the NG (red) except towards the end of evolution where the loop has
shrunk by a large fraction and the NG predicts a complete overlap of the
extent of the string on itself.

the comparison of the evolution of these primordial loops with their predicted
NG dynamics. The pictures demonstrate that the evolution is pretty much
identical in both cases all the way until moments before the predicted overlap
of the loop. There are no departures from the NG dynamics due to high
curvature regions for much of its evolution. This is easy to understand since
the loop is indeed much smoother due to the diffusion period. However, as the
loop shrinks, we start seeing some deviations from the NG behaviour, although
not so dramatic as in the non-primordial loops. As the loop comes close to its
overlap, the difference between both field theory and NG becomes more
apparent. This is also to be expected since the interaction of different
regions of the string in this pathological self-intersection is of course not
handled by the NG dynamics. Nevertheless, the fact that up to this point both
descriptions agree with one another can be seen as a validation of both the
field theory and NG reconstruction codes.

We also show in figure 8 the energy of this loop using the NG reconstruction
at each moment in time. We notice that the energy is pretty much constant
until $t\sim 40$, which is also the moment where there is the first signal of
deviation from NG dynamics in the loop’s evolution (see figs. 7). The
deviation from the NG prediction happens in several places in the loop and
even though they are associated with high curvature regions, they are not so
obvious, at least not visually, as the ones presented earlier in the non-
primordial loop. This is somewhat similar to what happens in figure 6. The
energy of the loop is decreasing, and yet the dynamics of the loop seems to
follow quite closely that of NG.

Figure 8: Total amount of invariant energy (total amount of $\sigma$) for the
NG reconstruction of the field theory data of the primordial loop.

## References

  * [1] H. Nielsen and P. Olesen, _Vortex-line models for dual strings_ , _Nuclear Physics B_ 61 (1973) 45.
  * [2] Y. Nambu _Lectures at the Copenhagen Summer Symposium_ (1970) .
  * [3] T. Goto, _Relativistic quantum mechanics of one-dimensional mechanical continuum and subsidiary condition of dual resonance model_ , _Prog. Theor. Phys._ 46 (1971) 1560.
  * [4] G. Vincent, N.D. Antunes and M. Hindmarsh, _Numerical simulations of string networks in the Abelian Higgs model_ , _Phys. Rev. Lett._ 80 (1998) 2277 [hep-ph/9708427].
  * [5] J.N. Moore, E.P.S. Shellard and C.J.A.P. Martins, _On the evolution of Abelian-Higgs string networks_ , _Phys. Rev. D_ 65 (2002) 023503 [hep-ph/0107171].
  * [6] N. Bevis, M. Hindmarsh, M. Kunz and J. Urrestilla, _CMB power spectrum contribution from cosmic strings using field-evolution simulations of the Abelian Higgs model_ , _Phys. Rev. D_ 75 (2007) 065015 [astro-ph/0605018].
  * [7] M. Hindmarsh, S. Stuckey and N. Bevis, _Abelian Higgs Cosmic Strings: Small Scale Structure and Loops_ , _Phys. Rev. D_ 79 (2009) 123504 [0812.1929].
  * [8] M. Hindmarsh, J. Lizarraga, J. Urrestilla, D. Daverio and M. Kunz, _Scaling from gauge and scalar radiation in Abelian Higgs string networks_ , _Phys. Rev. D_ 96 (2017) 023525 [1703.06696].
  * [9] T. Helfer, J.C. Aurrekoetxea and E.A. Lim, _Cosmic String Loop Collapse in Full General Relativity_ , _Phys. Rev. D_ 99 (2019) 104028 [1808.06678].
  * [10] J.R.C.C.C. Correia and C.J.A.P. Martins, _Abelian–Higgs cosmic string evolution with multiple GPUs_ , _Astron. Comput._ 34 (2021) 100438 [2005.14454].
  * [11] J.R.C.C.C. Correia and C.J.A.P. Martins, _High resolution calibration of the cosmic strings velocity dependent one-scale model_ , _Phys. Rev. D_ 104 (2021) 063511 [2108.07513].
  * [12] A. Albrecht and N. Turok, _Evolution of Cosmic Strings_ , _Phys. Rev. Lett._ 54 (1985) 1868.
  * [13] D.P. Bennett and F.R. Bouchet, _Cosmic string evolution_ , _Phys. Rev. Lett._ 63 (1989) 2776.
  * [14] B. Allen and E.P.S. Shellard, _Cosmic string evolution: a numerical simulation_ , _Phys. Rev. Lett._ 64 (1990) 119.
  * [15] V. Vanchurin, K.D. Olum and A. Vilenkin, _Scaling of cosmic string loops_ , _Phys. Rev. D_ 74 (2006) 063527 [gr-qc/0511159].
  * [16] V. Vanchurin, K. Olum and A. Vilenkin, _Cosmic string scaling in flat space_ , _Phys. Rev. D_ 72 (2005) 063514 [gr-qc/0501040].
  * [17] C.J.A.P. Martins and E.P.S. Shellard, _Fractal properties and small-scale structure of cosmic string networks_ , _Phys. Rev. D_ 73 (2006) 043515 [astro-ph/0511792].
  * [18] C. Ringeval, M. Sakellariadou and F. Bouchet, _Cosmological evolution of cosmic string loops_ , _JCAP_ 02 (2007) 023 [astro-ph/0511646].
  * [19] K.D. Olum and V. Vanchurin, _Cosmic string loops in the expanding Universe_ , _Phys. Rev. D_ 75 (2007) 063521 [astro-ph/0610419].
  * [20] J.J. Blanco-Pillado, K.D. Olum and B. Shlaer, _Large parallel cosmic string simulations: New results on loop production_ , _Phys. Rev. D_ 83 (2011) 083514 [1101.5173].
  * [21] J.J. Blanco-Pillado, K.D. Olum and B. Shlaer, _The number of cosmic string loops_ , _Phys. Rev. D_ 89 (2014) 023512 [1309.6637].
  * [22] K.D. Olum and J.J. Blanco-Pillado, _Field theory simulation of Abelian Higgs cosmic string cusps_ , _Phys. Rev. D_ 60 (1999) 023503 [gr-qc/9812040].
  * [23] K.D. Olum and J.J. Blanco-Pillado, _Radiation from cosmic string standing waves_ , _Phys. Rev. Lett._ 84 (2000) 4288 [astro-ph/9910354].
  * [24] M. Hindmarsh, J. Lizarraga, A. Urio and J. Urrestilla, _Loop decay in Abelian-Higgs string networks_ , _Phys. Rev. D_ 104 (2021) 043519 [2103.16248].
  * [25] H. Arodz, _Bound states of the vector field with a vortex in the Abelian Higgs model_ , _Acta Phys. Polon._ B22 (1991) 511.
  * [26] M. Goodband and M. Hindmarsh, _Bound states and instabilities of vortices_ , _Phys. Rev._ D52 (1995) 4621.
  * [27] A. Alonso-Izquierdo, W. Garcia Fuertes and J. Mateos Guilarte, _A note on BPS vortex bound states_ , _Phys. Lett._ B753 (2016) 29 [1509.06632].
  * [28] T. Kojo, H. Suganuma and K. Tsumura, _Peristaltic modes of single vortex in the Abelian Higgs model_ , _Phys. Rev._ D75 (2007) 105015 [hep-ph/0702014].
  * [29] J.J. Blanco-Pillado, D. Jiménez-Aguilar and J. Urrestilla, _Exciting the domain wall soliton_ , _JCAP_ 01 (2021) 027 [2006.13255].
  * [30] J.J. Blanco-Pillado, D. Jiménez-Aguilar, J.M. Queiruga and J. Urrestilla, _The dynamics of Domain Wall Strings_ , 2209.12945.
  * [31] J.J. Blanco-Pillado, D. Jiménez-Aguilar, J.M. Queiruga and J. Urrestilla, _Parametric Resonances in Axionic Cosmic Strings_ , 2212.06194.
  * [32] J.J. Blanco-Pillado, D. Jiménez-Aguilar, J.M. Queiruga and J. Urrestilla, _Internal excitations of global vortices_ , _JCAP_ 10 (2021) 047 [2107.02215].
  * [33] N.S. Manton and H. Merabet, _Phi**4 kinks: Gradient flow and dynamics_ , _Nonlinearity_ 10 (1997) 3 [hep-th/9605038].
  * [34] D. Forster, _Dynamics of Relativistic Vortex Lines and their Relation to Dual Theory_ , _Nucl. Phys. B_ 81 (1974) 84.
  * [35] R.A. Matzner, _Interaction of U(1) cosmic strings: Numerical intercommutation_ , _Comput. Phys._ 2 (1988) 51.
  * [36] D. Matsunami, L. Pogosian, A. Saurabh and T. Vachaspati, _Decay of Cosmic String Loops Due to Particle Radiation_ , _Phys. Rev. Lett._ 122 (2019) 201301 [1903.05102].
  * [37] A. Vilenkin and E.S. Shellard, _Cosmic Strings and Other Topological Defects_ , Cambridge University Press (7, 2000).
  * [38] J.J. Blanco-Pillado and K.D. Olum, _Stochastic gravitational wave background from smoothed cosmic string loops_ , _Phys. Rev. D_ 96 (2017) 104046 [1709.02693].
  * [39] T.W.B. Kibble and N. Turok, _Selfintersection of Cosmic Strings_ , _Phys. Lett. B_ 116 (1982) 141.


Fokker Planck equation

The research of Hassani and Shen is supported by NSF Grants 1837253, 1943064, AFOSR Grant FA9550-20-1-0111, DCIST-CRA, and the AI Institute for Learning-Enabled Optimization at Scale (TILOS).
Zhenfu Wang is supported by the National Key R&D Program of China, Project Number 2021YFA1002800, NSFC grant No.12171009, Young Elite Scientist Sponsorship Program by China Association for Science and Technology (CAST) No. YESS20200028 and the start-up fund from Peking University.
Amin Karbasi acknowledges funding in direct support of this work from NSF (IIS-1845032), ONR (N00014-19-1-2406),  NSF (2112665), and the AI Institute for Learning-Enabled Optimization at Scale (TILOS).

[Ambrosio et al., 2005]
Luigi Ambrosio, Nicola Gigli, and Giuseppe Savaré.
Gradient flows: in metric spaces and in the space of
probability measures.
Springer Science & Business Media, 2005.

[Bailo et al., 2018]
Rafael Bailo, Jose A Carrillo, and Jingwei Hu.
Fully discrete positivity-preserving and energy-dissipating schemes
for aggregation-diffusion equations with a gradient flow structure.
arXiv preprint arXiv:1811.11502, 2018.

[Bernton, 2018]
Espen Bernton.
Langevin monte carlo and jko splitting.
In Conference On Learning Theory, pages 1777–1798. PMLR,

[Blechschmidt and Ernst, 2021]
Jan Blechschmidt and Oliver G Ernst.
Three ways to solve partial differential equations with neural
networks—a review.
GAMM-Mitteilungen, 440 (2):0 e202100006, 2021.

[Cardaliaguet and Porretta, 2020]
Pierre Cardaliaguet and Alessio Porretta.
An introduction to mean field game theory.
In Mean Field Games, pages 1–158. Springer, 2020.

[Carrillo et al., 2015]
José A Carrillo, Alina Chertock, and Yanghong Huang.
A finite-volume method for nonlinear nonlocal equations with a
gradient flow structure.
Communications in Computational Physics, 170
(1):0 233–258, 2015.

[Carrillo et al., 2021]
Jose A Carrillo, Katy Craig, Li Wang, and Chaozhen Wei.
Primal dual methods for wasserstein gradient flows.
Foundations of Computational Mathematics, pages 1–55, 2021.

[Chen et al., 2018]
Ricky TQ Chen, Yulia Rubanova, Jesse Bettencourt, and David K Duvenaud.
Neural ordinary differential equations.
Advances in neural information processing systems, 31, 2018.

[Chizat and Bach, 2018]
Lenaic Chizat and Francis Bach.
On the global convergence of gradient descent for over-parameterized
models using optimal transport.
Advances in neural information processing systems, 31, 2018.

[Choromanski et al., 2020]
Krzysztof M Choromanski, Jared Quincy Davis, Valerii Likhosherstov, Xingyou
Song, Jean-Jacques Slotine, Jacob Varley, Honglak Lee, Adrian Weller, and
Vikas Sindhwani.
Ode to an ode.
Advances in Neural Information Processing Systems,
33:0 3338–3350, 2020.

[Chourasia et al., 2021]
Rishav Chourasia, Jiayuan Ye, and Reza Shokri.
Differential privacy dynamics of langevin diffusion and noisy
gradient descent.
Advances in Neural Information Processing Systems, 34, 2021.

[Dalalyan, 2017]
Arnak S Dalalyan.
Theoretical guarantees for approximate sampling from smooth and
log-concave densities.
Journal of the Royal Statistical Society: Series B (Statistical
Methodology), 790 (3):0 651–676, 2017.

[Degond and Mustieles, 1990]
Pierre Degond and Francisco-José Mustieles.
A deterministic approximation of diffusion equations using particles.
SIAM Journal on Scientific and Statistical Computing,
110 (2):0 293–310, 1990.

[Dupont et al., 2019]
Emilien Dupont, Arnaud Doucet, and Yee Whye Teh.
Augmented neural odes.
Advances in Neural Information Processing Systems, 32, 2019.

[Fang et al., 2021]
Cong Fang, Jason Lee, Pengkun Yang, and Tong Zhang.
Modeling from features: a mean-field framework for over-parameterized
deep neural networks.
In Conference on learning theory, pages 1887–1936. PMLR,

[Gomes et al., 2014]
Diogo A Gomes et al.
Mean field games models—a brief survey.
Dynamic Games and Applications, 40 (2):0
110–154, 2014.

[Han et al., 2018]
Jiequn Han, Arnulf Jentzen, and E Weinan.
Solving high-dimensional partial differential equations using deep
Proceedings of the National Academy of Sciences, 1150
(34):0 8505–8510, 2018.

[Hashimoto et al., 2016]
Tatsunori Hashimoto, David Gifford, and Tommi Jaakkola.
Learning population-level diffusions with generative rnns.
In International Conference on Machine Learning, pages
2417–2426. PMLR, 2016.

[Jabin and Wang, 2016]
Pierre-Emmanuel Jabin and Zhenfu Wang.
Mean field limit and propagation of chaos for vlasov systems with
bounded forces.
Journal of Functional Analysis, 2710 (12):0
3588–3627, 2016.
ISSN 0022-1236.

[Jordan et al., 1998]
Richard Jordan, David Kinderlehrer, and Felix Otto.
The variational formulation of the fokker–planck equation.
SIAM Journal on Mathematical Analysis, 290
(1):0 1–17, 1998.

[Li et al., 2021]
Ruilin Li, Hongyuan Zha, and Molei Tao.
Sqrt (d) dimension dependence of langevin monte carlo.
arXiv preprint arXiv:2109.03839, 2021.

[Li et al., 2019]
Xuechen Li, Yi Wu, Lester Mackey, and Murat A Erdogdu.
Stochastic runge-kutta accelerates langevin monte carlo and beyond.
Advances in neural information processing systems, 32, 2019.

[Liu et al., 2019]
Chang Liu, Jingwei Zhuo, and Jun Zhu.
Understanding mcmc dynamics as flows on the wasserstein space.
In International Conference on Machine Learning, pages
4093–4103. PMLR, 2019.

[Liu et al., 2020]
Shu Liu, Wuchen Li, Hongyuan Zha, and Haomin Zhou.
Neural parametric fokker-planck equations.
arXiv preprint arXiv:2002.11309, 2020.

[Long et al., 2018]
Zichao Long, Yiping Lu, Xianzhong Ma, and Bin Dong.
Pde-net: Learning pdes from data.
In International Conference on Machine Learning, pages
3208–3216. PMLR, 2018.

[Long et al., 2019]
Zichao Long, Yiping Lu, and Bin Dong.
Pde-net 2.0: Learning pdes from data with a numeric-symbolic hybrid
deep network.
Journal of Computational Physics, 399:0 108925, 2019.

[Lucia and Gervino, 2015]
Umberto Lucia and Gianpiero Gervino.
Fokker-planck equation and thermodynamic system analysis.
Entropy, 170 (2):0 763–771, 2015.

[Markowich and Villani, 2000]
Peter A Markowich and Cédric Villani.
On the trend to equilibrium for the fokker-planck equation: an
interplay between physics and functional analysis.
Mat. Contemp, 19:0 1–29, 2000.

[Mokrov et al., 2021]
Petr Mokrov, Alexander Korotin, Lingxiao Li, Aude Genevay, Justin M Solomon,
and Evgeny Burnaev.
Large-scale wasserstein gradient flows.
Advances in Neural Information Processing Systems, 34, 2021.

[Qi and Majda, 2016]
Di Qi and Andrew J Majda.
Low-dimensional reduced-order models for statistical response and
uncertainty quantification: Two-layer baroclinic turbulence.
Journal of the Atmospheric Sciences, 730
(12):0 4609–4639, 2016.

[Raissi et al., 2019]
Maziar Raissi, Paris Perdikaris, and George E Karniadakis.
Physics-informed neural networks: A deep learning framework for
solving forward and inverse problems involving nonlinear partial differential
Journal of Computational physics, 378:0 686–707,

[Sirignano and Spiliopoulos, 2020]
Justin Sirignano and Konstantinos Spiliopoulos.
Mean field analysis of neural networks: A central limit theorem.
Stochastic Processes and their Applications, 1300
(3):0 1820–1852, 2020.

[Sohl-Dickstein et al., 2015]
Jascha Sohl-Dickstein, Eric Weiss, Niru Maheswaranathan, and Surya Ganguli.
Deep unsupervised learning using nonequilibrium thermodynamics.
In International Conference on Machine Learning, pages
2256–2265. PMLR, 2015.

[Solin et al., 2021]
Arno Solin, Ella Tamir, and Prakhar Verma.
Scalable inference in sdes by direct matching of the
fokker–planck–kolmogorov equation.
Advances in Neural Information Processing Systems, 34, 2021.

[Song et al., 2020]
Yang Song, Jascha Sohl-Dickstein, Diederik P Kingma, Abhishek Kumar, Stefano
Ermon, and Ben Poole.
Score-based generative modeling through stochastic differential
In International Conference on Learning Representations, 2020.

[Sonoda and Murata, 2019]
Sho Sonoda and Noboru Murata.
Transport analysis of infinitely deep neural network.
The Journal of Machine Learning Research, 200
(1):0 31–82, 2019.

[Villani, 2009]
Cédric Villani.
Optimal transport: old and new, volume 338.
Springer, 2009.