# Description of inclusive $(d,d^{\prime}x)$ reaction with the semiclassical
distorted wave model

Hibiki Nakada<EMAIL_ADDRESS>Research Center for Nuclear Physics
(RCNP), Osaka University, Ibaraki, Osaka 567-0047, Japan  Kazuki Yoshida
Advanced Science Research Center, Japan Atomic Energy Agency, Tokai, Ibaraki
319-1195, Japan  Kazuyuki Ogata  Department of Physics, Kyushu University,
Fukuoka 819-0395, Japan  Research Center for Nuclear Physics (RCNP), Osaka
University, Ibaraki, Osaka 567-0047, Japan

###### Abstract

Background: The description of deuteron-induced inclusive reactions has been
an important subject in direct nuclear reaction studies and nuclear data
science. For proton-induced inclusive processes, the semiclassical distorted
wave model (SCDW) is one of the most successful models based on quantum
mechanics.

Purpose: We improve SCDW for deuteron-induced inclusive processes and clarify
the importance of the proper treatment of the kinematics of the deuteron
inside a nucleus.

Methods: The double differential cross section (DDX) of the inclusive
deuteron-emission process $(d,d^{\prime}x)$ is described by one-step SCDW. The
changes in the kinematics due to the distortion effect, the refraction effect,
is taken into account by the local semiclassical approximation (LSCA).

Results: The calculated DDXs of $(d,d^{\prime}x)$ reasonably reproduce
experimental data in the small energy-transfer region and at forward and
middle angles with some exceptions. The angular distributions of
$(d,d^{\prime}x)$ are improved by including the refraction effect.

Conclusion: The proper treatment of the changes in the kinematics of the
deuteron inside a nucleus is necessary in describing the ($d$,$d^{\prime}x$)
reaction. The effect of the changes on the DDX of $(d,d^{\prime}x)$ is
significant compared to on the proton-induced inclusive process
$(p,p^{\prime}x)$ because of the stronger distortion effect on the deuteron.

††preprint: APS/123-QED

##  I INTRODUCTION

Deuteron has the smallest binding energy among all stable nuclei. As
originated from the idea of Butler [1], the weakly-bound nature of deuteron
has been utilized for carrying out one-nucleon transfer reactions to study the
single-particle (s.p.) structure of nuclei [2]. Furthermore, deuteron-induced
reactions have opened many physics cases to reveal three-body dynamics of
reaction systems in which a fragile nucleus is involved [3, 4, 5, 6, 7, 8].
Roles of deuteron breakup channels, in which proton and neutron are in
continuum states, have intensively been investigated.

The fragileness of deuteron is also important for nuclear data science. The
international fusion materials irradiation facility (IFMIF) [9], which aims at
using the inclusive $(d,nx)$ reaction at 40 MeV as an intense neutron source,
is one of the most well-known international scientific projects using deuteron
accelerator. The central idea of IFMIF is that the incident deuteron is broken
up by interacting with the target and intense neutron with about half the
deuteron incident energy is emitted; statistical decay after forming a
compound nucleus is also considered to contribute to the neutron emission for
large energy transfer. Quite recently, an integrated code system describing
deuteron-induced reactions, which is designated as DEURACS, has been
constructed and successfully applied to analysis of $(d,nx)$ reaction data
[10, 11, 12]. It was found that the description of deuteron breakup channels
is of crucial importance for accurately evaluating the amount of the emitted
neutron, its angular and energy distribution in particular.

From the viewpoint of direct nuclear reaction study, the most challenging part
for describing ($d,nx$) is the deuteron breakup with exciting the target
nucleus A, which is called the nonelastic breakup (NEB). NEB contains a huge
number of final states of A and it is almost impossible to describe each
nuclear state accurately. DEURACS employs the Glauber model [13] to circumvent
the difficulty; the eikonal and adiabatic approximations allow one to describe
NEB as a combination of neutron elastic and proton nonelastic processes, and
the latter can easily be evaluated by using the closure property of the proton
scattering matrix [14, 15]. The validity of the Glauber model is, however,
rather questionable at low incident energy and/or for large momentum and
energy transfer. In fact, the agreement between the result of DEURACS and
experimental data for ($d,nx$) at middle emission angles is slightly flawed
compared with that at forward angles [12]. Although the neutron emission cross
section is forward-peaked and the “deviation” is not very serious for
practical use, the description of NEB of deuteron without using the eikonal
and adiabatic approximations will be an important subject of nuclear reaction
study. Recently, the Ichimura-Austern-Vincent (IAV) model [16] has
successfully been applied to NEB in several cases [8, 17, 18]. It should be
noted, however, that in the IAV model for ($d,nx$), the kinematics of the
neutron are not affected at all by the nonelastic processes for which the
proton and A undergo. In this sense, the three-body kinematics are not treated
in a fully consistent manner in the IAV model.

On the other hand, for proton-induced inclusive processes, ($p,p^{\prime}x$),
several quantum-mechanical models [19, 20, 21, 22] have been developed and
successfully reproduced experimental data. Among them, the semiclassical
distorted wave model (SCDW) [22, 23, 24, 25, 26, 27] has no free adjustable
parameter and allows a simple intuitive picture of ($p,p^{\prime}x$). The
original SCDW adopted the local Fermi-gas model (LFG) for initial and final
nuclear s.p. states. Although LFG will be totally unrealistic for modeling
specific nuclear states, it will reasonably describe the total response of a
nucleus to which many initial and final states contribute. It should be noted
that, in SCDW, there is no kinematical assumption or restriction for the
reaction particles. This idea for treating processes via a huge number of
nuclear states is expected to work also for deuteron-induced reactions. Note
that the latest version of SCDW adopts the Wigner transform of one-body
density matrices calculated with a s.p. model for nuclei [26] instead of LFG;
for reducing numerical task, we use LFG in this work.

The main purpose of this study is to extend SCDW to deuteron-induced inclusive
processes. Although our ultimate goal is to describe ($d,nx$), as the first
step, we focus on the inclusive deuteron-emission process ($d,d^{\prime}x$).
We assume for simplicity that scattering waves of the incoming and outgoing
deuteron can be described with a phenomenological optical potential, meaning
that deuteron breakup channels are not explicitly treated but net loss of
probability is implicitly taken into account by the absorption effect of the
deuteron optical potential. On the other hand, as in SCDW studies on
($p,p^{\prime}x$), we respect kinematics of deuteron inside a nucleus, by
using the local semiclassical approximation (LSCA) [22] to the deuteron
distorted waves. We clarify how the proper treatment of the “refraction” of
deuteron by the distorting potential is important to describe
($d,d^{\prime}x$) experimental data. We include only the one-step process and
mainly discuss the small energy-transfer region.

The construction of this paper is as follows. In Sec. II we describe SCDW for
the inclusive $(d,d^{\prime}x)$ reaction, applying LFG and LSCA. In Sec. III
we compare the calculated DDXs of the inclusive $(d,d^{\prime}x)$ reaction
with experimental data and demonstrate the effect of nuclear refraction.
Finally, a summary is given in Sec. IV.

##  II FORMALISM

We describe the inclusive $(d,d^{\prime}x)$ reaction by one-step SCDW. The
foundation of SCDW is the DWBA series expansion of the transition matrix (T
matrix). The T-matrix element, for which the target nucleus is excited from
the initial s.p. state $\phi_{\alpha}$ to the final one $\phi_{\beta}$, is
given by

$\displaystyle
T_{\beta\alpha}={\Braket{\chi^{(-)}_{f}(\bm{r}_{0})\phi_{\beta}(\bm{r})}{v(\bm{r}_{0}-\bm{r})}{\chi^{(+)}_{i}(\bm{r}_{0})\phi_{\alpha}(\bm{r})}},$
(1)

where $\bm{r}_{0}$ and $\bm{r}$ are the coordinates of the incident deuteron
and the nucleon inside the target, respectively. $\chi_{i}$ ($\chi_{f}$) is
the distorted wave for the deuteron in the initial (final) state. The
superscripts $(+)$ and $(-)$ denote the outgoing and incoming boundary
conditions for $\chi$, respectively. $v$ is the effective interaction between
the deuteron and the target nucleon. The double differential cross section
(DDX) for the emitted deuteron energy $E_{f}$ and the solid angle $\Omega_{f}$
is given by

$\displaystyle\frac{\partial^{2}\sigma}{\partial
E_{f}\partial\Omega_{f}}=C\frac{k_{f}}{k_{i}}\sum_{\alpha,\beta}|T_{\beta\alpha}|^{2}\delta(E_{i}+\varepsilon_{\alpha}-E_{f}-\varepsilon_{\beta}),$
(2)

where $C=4\mu^{2}/(2\pi\hbar^{2})^{2}$, $E_{i}$ is the deuteron incident
energy, $\mu$ is the reduced mass between the deuteron and the target nucleus
and $k_{i}$ ($k_{f}$) is the asymptotic momentum of the incident (emitted)
deuteron. $\varepsilon_{\gamma}~{}(\gamma=\alpha~{}\mathrm{or}~{}\beta)$ is
the kinetic energy of the target nucleon. The summation is taken over all the
initial and the final s.p. states, $\alpha$ and $\beta$, which are relevant to
the inclusive $(d,d^{\prime}x)$ reaction. On expanding the squared modulus in
Eq. (2), one obtains

$\displaystyle\frac{\partial^{2}\sigma}{\partial
E_{f}\partial\Omega_{f}}=C\frac{k_{f}}{k_{i}}\int
d\bm{r}_{0}d\bm{r}\chi^{\ast(-)}_{f}(\bm{r}_{0})v(\bm{r}_{0}-\bm{r})\chi^{(+)}_{i}(\bm{r}_{0})$
$\displaystyle\times\int
d\bm{r}_{0}^{\prime}d\bm{r}^{\prime}\chi^{(-)}_{f}(\bm{r}^{\prime}_{0})v^{\ast}(\bm{r}^{\prime}_{0}-\bm{r}^{\prime})\chi^{\ast(+)}_{i}(\bm{r}^{\prime}_{0})K(\bm{r},\bm{r}^{\prime}),$
(3)

where the kernel $K(\bm{r},\bm{r}^{\prime})$ is define by

$\displaystyle K(\bm{r},\bm{r}^{\prime})$
$\displaystyle\equiv\sum_{\alpha}\phi_{\alpha}(\bm{r})\phi^{\ast}_{\alpha}(\bm{r}^{\prime})\sum_{\beta}\phi^{\ast}_{\beta}(\bm{r})\phi_{\beta}(\bm{r}^{\prime})$
$\displaystyle\times\delta(E_{i}+\varepsilon_{\alpha}-E_{f}-\varepsilon_{\beta}).$
(4)

When a large number of s.p. states are involved, $K(\bm{r},\bm{r}^{\prime})$
becomes a short-ranged function of $|\bm{r}-\bm{r}^{\prime}|$ [24, 23, 28].
The center-of-mass and relative coordinates of the $d$-$N$ system, $\bm{R}$
and $\bm{s}$, respectively, are given by

$\displaystyle\bm{R}$
$\displaystyle=\frac{A_{d}}{A_{d}+1}\bm{r}_{0}+\frac{1}{A_{d}+1}\bm{r},$ (5)
$\displaystyle\bm{s}$ $\displaystyle=\bm{r}_{0}-\bm{r}.$ (6)

Inversely, $\bm{r}_{0}$ and $\bm{r}$ are written as

$\displaystyle\bm{r}_{0}=\bm{R}+\frac{1}{A_{d}+1}\bm{s},$ (7)
$\displaystyle\bm{r}=\bm{R}-\frac{A_{d}}{A_{d}+1}\bm{s},$ (8)

where $A_{d}$ is the mass number of deuteron, i.e., $A_{d}=2$. With the
coordinates $\bm{R}$ and $\bm{s}$, one can rewrite Eq. (3) as

$\displaystyle\frac{\partial^{2}\sigma}{\partial E_{f}\partial\Omega_{f}}$
$\displaystyle=C\frac{k_{f}}{k_{i}}\int
d\bm{R}\,d\bm{s}\,d\bm{R}^{\prime}\,d\bm{s}^{\prime}$
$\displaystyle\times\chi^{\ast(-)}_{f}(\bm{R}+\bm{s}/3)v(\bm{s})\chi^{(+)}_{i}(\bm{R}+\bm{s}/3)$
$\displaystyle\times\chi^{(-)}_{f}(\bm{R}^{\prime}+\bm{s}^{\prime}/3)v^{\ast}(\bm{s}^{\prime})\chi^{\ast(+)}_{i}(\bm{R}^{\prime}+\bm{s}^{\prime}/3)$
$\displaystyle\times K(\bm{R},\bm{s},\bm{R}^{\prime},\bm{s}^{\prime}),$ (9)

where

$\displaystyle K(\bm{R},\bm{s},\bm{R}^{\prime},\bm{s}^{\prime})$
$\displaystyle=\sum_{\alpha}\phi_{\alpha}(\bm{R}-2\bm{s}/3)\phi^{\ast}_{\alpha}(\bm{R}^{\prime}-2\bm{s}^{\prime}/3)$
$\displaystyle\times\sum_{\beta}\phi^{\ast}_{\beta}(\bm{R}-2\bm{s}/3)\phi_{\beta}(\bm{R}^{\prime}-2\bm{s}^{\prime}/3)$
$\displaystyle\times\delta(E_{i}+\varepsilon_{\alpha}-E_{f}-\varepsilon_{\beta}).$
(10)

Here, we make two approximations to Eq. (9). One is LFG for nuclear states and
the other is LSCA for the distorted waves as mentioned in Sec. I. In LFG,
$\phi_{\gamma}~{}(\gamma=\alpha~{}\mathrm{or}~{}\beta)$ is approximated by the
plane wave with momentum $k_{\gamma}$ within a smaller cell-size of $|\bm{s}|$
than the range of $v$. The summation of $\gamma$ is then expressed as an
integral over $k_{\gamma}$, where the threshold momentum between channels
$\alpha$ and $\beta$ is the local Fermi momentum $k_{F}(\bm{R})$, which is
related to the nuclear density $\rho(\bm{R})$ through

$\displaystyle\rho(\bm{R})=4\frac{4\pi}{3}\frac{k^{3}_{F}(\bm{R})}{(2\pi)^{3}}.$
(11)

In LSCA, the short-range propagation of the distorted wave
$\chi_{c}~{}(c=i~{}\mathrm{or}~{}f)$ from a reference point $\bm{R}$ is
approximated by the plane wave, i.e.,

$\displaystyle\chi_{c}(\bm{R}+\bm{s}/3)\simeq\chi_{c}(\bm{R})e^{i\bm{k}_{c}(\bm{R})\cdot\bm{s}/3}.$
(12)

This approximation is valid because the range of the $d$-$N$ interaction $v$
is short and therefore only a small $\bm{s}$ is relevant to the reaction. In
Eq. (12), $\bm{k}_{c}(\bm{R})$ is the real local momentum of the deuteron. The
direction of $\bm{k}_{c}(\bm{R})$ is taken to be the same as that of the flux
of the distorted wave $\chi_{c}(\bm{R})$. The magnitude of
$\bm{k}_{c}(\bm{R})$ is given by that of the real part of the complex local
momentum $\bm{K}_{c}(\bm{R})$ satisfying the local energy conservation [24]

$\displaystyle\frac{\hbar^{2}k_{c}^{2}}{2\mu}=\frac{\hbar^{2}K^{2}_{c}(\bm{R})}{2\mu}+U_{c}(\bm{R}),$
(13)

where $U_{c}(\bm{R})~{}(c=i~{}\mathrm{or}~{}f)$ is a complex distorting
potential for the deuteron. In present calculation, we use the optical
potential for $U_{c}(\bm{R})$ as mentioned in Sec. III.1.

Using LFG and LSCA, one can rewrite Eq. (9) as

$\displaystyle\frac{\partial^{2}\sigma}{\partial E_{f}\partial\Omega_{f}}$
$\displaystyle=\frac{C}{(2\pi)^{3}}\frac{k_{f}}{k_{i}}\int
d\bm{R}|\chi^{(-)}_{f}(\bm{R})|^{2}|\chi^{(+)}_{i}(\bm{R})|^{2}$
$\displaystyle\times\int_{k_{\alpha}\leq
k_{F}(\bm{R})}d\bm{k}_{\alpha}\int_{k_{\beta}>k_{F}(\bm{R})}d\bm{k}_{\beta}$
$\displaystyle\times\left|\int
d\bm{s}v(\bm{s})e^{-i\bm{q}(\bm{R})\cdot\bm{s}}\right|^{2}$
$\displaystyle\times\delta(\bm{k}_{i}(\bm{R})+\bm{k}_{\alpha}-\bm{k}_{f}(\bm{R})-\bm{k}_{\beta})$
$\displaystyle\times\delta(E_{i}+\varepsilon_{\alpha}-E_{f}-\varepsilon_{\beta}),$
(14)

where $\bm{q}(\bm{R})$ is the local momentum transfer defined by
$\bm{k}_{i}(\bm{R})-\bm{k}_{f}(\bm{R})$. In Eq. (14), Dirac’s delta functions
and the ranges of the integrations, $k_{\alpha}\leq k_{F}(\bm{R})$ and
$k_{\beta}>k_{F}(\bm{R})$, guarantee that the $d$-$N$ elementary process
satisfies the Pauli principle and the energy and local momentum conservation
in the $(d,d^{\prime}x)$ reaction. We make the on-the-energy-shell
approximation to the squared modulus of the matrix element of $v$:

$\displaystyle\frac{\mu^{2}_{dN}}{(2\pi\hbar^{2})^{2}}\left|\int
d\bm{s}v(\bm{s})e^{-i\bm{q}(\bm{R})\cdot\bm{s}}\right|^{2}\simeq\left(\frac{d\sigma_{dN}}{d\Omega}\right)_{\theta_{dN}(\bm{R}),E_{dN}(\bm{R})},$
(15)

where $\mu_{dN}$ is the reduced mass of the $d$-$N$ system.
$\theta_{dN}(\bm{R})$ is the local $d$-$N$ scattering angle between the
initial relative momentum $\bm{\kappa}(\bm{R})$ and the final one
$\bm{\kappa}^{\prime}(\bm{R})$, which are defined by

$\displaystyle\bm{\kappa}(\bm{R})$
$\displaystyle\equiv\frac{1}{A_{d}+1}\bm{k}_{i}(\bm{R})-\frac{A_{d}}{A_{d}+1}\bm{k}_{\alpha},$
(16) $\displaystyle\bm{\kappa}^{\prime}(\bm{R})$
$\displaystyle\equiv\frac{A_{d}}{A_{d}+1}\bm{k}_{f}(\bm{R})-\frac{1}{A_{d}+1}\bm{k}_{\beta}.$
(17)

The local $d$-$N$ scattering energy $E_{\bm{R}}$ is defined by

$\displaystyle E_{dN}(\bm{R})=\frac{\hbar^{2}\kappa^{2}(\bm{R})}{2\mu_{dN}}.$
(18)

By substituting Eq. (15) for Eq. (14) and integrating over $\bm{k}_{\beta}$,
one obtains the following closed form of the DDX of the inclusive
$(d,d^{\prime}x)$ reaction:

$\displaystyle\frac{\partial^{2}\sigma}{\partial E_{f}\partial\Omega_{f}}$
$\displaystyle=\left[\frac{A_{d}A}{A_{d}+A}\right]^{2}\frac{k_{f}}{k_{i}}\int
d\bm{R}$
$\displaystyle\times|\chi^{(-)}_{f}(\bm{R})|^{2}|\chi^{(+)}_{i}(\bm{R})|^{2}\left[\frac{\partial^{2}\sigma}{\partial
E_{f}\partial\Omega_{f}}\right]_{\bm{R}}\rho(\bm{R}),$ (19)

where $A$ is the mass number of the target nucleus. The DDX of the elementary
process averaged over $\bm{k}_{\alpha}$ at $\bm{R}$ in the Fermi sphere
characterized by $k_{F}^{3}(\bm{R})$ is given by

$\displaystyle\left[\frac{\partial^{2}\sigma}{\partial
E_{f}\partial\Omega_{f}}\right]_{\bm{R}}$
$\displaystyle=\frac{1}{(4\pi/3)k^{3}_{F}(\bm{R})}\left[\frac{A_{d}+1}{A_{d}}\right]^{2}$
$\displaystyle\times\int_{k_{\alpha}\leq
k_{F}(\bm{R})}d\bm{k}_{\alpha}\left(\frac{d\sigma_{dN}}{d\Omega}\right)_{\theta_{dN}(\bm{R}),E_{dN}(\bm{R})}$
$\displaystyle\times\delta(E_{i}+\varepsilon_{\alpha}-E_{f}-\varepsilon_{\beta}).$
(20)

LSCA can incorporate the distortion effect on the kinematics of the incident
and emitted particles. This effect can be regarded as refraction due to the
distorting potential because the direction of the local momentum changes
continuously as a function of ${\bm{R}}$. To clarify the refraction effect, we
also consider the asymptotic momentum approximation (AMA), which replaces
$\bm{k}_{c}(\bm{R})$ with $\bm{k}_{c}$, i.e.,
$\bm{k}_{c}(\bm{R})\rightarrow\bm{k}_{c}$ in Eq. (12).

$\displaystyle\chi_{c}(\bm{R}+\bm{s}/3)\simeq\chi_{c}(\bm{R})e^{i\bm{k}_{c}\cdot\bm{s}/3}.$
(21)

If AMA shown in Eq. (21) is used instead of LSCA, $\bm{k}_{c}(\bm{R})$ is
constant with respect to $\bm{R}$ in Eqs. (13)-(20). The effect of the
refraction is discussed in Sec. III.3. The validity of LSCA and AMA is given
in Appendix.

Figure 1: Comparison of the experimental data and calculated DDXs of the
inclusive $(d,d^{\prime}x)$ reaction on ${}^{27}\mathrm{Al}$ and
${}^{58}\mathrm{Ni}$ at 100 MeV and 80 MeV, ${}^{208}\mathrm{Pb}$ at 100 MeV
and 70 MeV, ${}^{93}\mathrm{Nb}$ and ${}^{181}\mathrm{Ta}$ at 100 MeV, and
${}^{90}\mathrm{Zr}$ and ${}^{232}\mathrm{Th}$ at 70 MeV for different
deuteron emission angles. The solid (dashed) lines represent the DDXs with
LSCA (AMA). The experimental data at 100 MeV are taken from Ref. [29] and
those at 80 and 70 MeV are from Ref. [30].

##  III RESULTS AND DISCUSSION

###  III.1 Numerical inputs

We assume the Woods-Saxon shaped global optical potential by An and Cai [31]
for the distorting potential $U_{c}(\bm{R})$ of the deuteron scattering off
target nuclei. The effect of the nonlocality of the deuteron distorting
potentials is taken into account by multiplying the scattering waves by the
Perey factor [32] $F_{c}(R)=[1-\mu\beta^{2}/(2\hbar^{2})U_{c}(R)]^{-1/2}$,
where $\mu$ is the reduced mass of the deuteron and the target. The range of
nonlocality $\beta$ for the deuteron is taken to be 0.54 [33]. We assume the
Woods-Saxon form for the nuclear density as

$\displaystyle\rho(R)=\frac{\rho_{0}}{1+\mathrm{exp}\left(\frac{R-R_{\rho}}{a_{\rho}}\right)},$
(22)

where the radial parameter is given by $R_{\rho}=r_{\rho}A^{1/3}$ with
$r_{\rho}=1.15$ fm, the diffuseness parameter is set to $a_{\rho}=0.5$ fm as
in Ref. [22], and $A$ being the mass number of the target nucleus. The
constant $\rho_{0}$ is determined to normalize the integrated value of
$\rho(R)$ to $A$. The local Fermi momentum is calculated from the nucleon
density as in Eq. (11).

For the free $d$-$N$ scattering cross section, we use the numerical table
fitted with several Gaussian functions to reproduce the experimental data of
$p$-$d$ scattering from 5 to 800 MeV [34]. In this table, the cross section
does not diverge at $0^{\circ}$ because we neglect the Coulomb elastic
scattering. For the free $p$-$N$ scattering cross section used in the
calculation of the $(p,p^{\prime}x)$ process, we use the nucleon-nucleon t
matrix by Franey and Love [35, 36].

###  III.2 Results of the SCDW calculation for $(d,d^{\prime}x)$ reactions and
comparison with data

We show the DDXs of the inclusive $(d,d^{\prime}x)$ reaction calculated with
SCDW using LSCA and AMA, and compare them with experimental data. Below we
discuss the DDXs as a function of the emission energy $E_{f}$ of the deuteron
with fixed $\Omega_{f}$. Figure 1 shows the calculated DDXs as a function of
$E_{f}$ at the several laboratory scattering angles $\theta$. The solid
(dashed) lines represent the DDXs using LSCA (AMA). The DDXs are calculated
for ${}^{27}\mathrm{Al}$ and ${}^{58}\mathrm{Ni}$ at the incident energy
$E_{i}=100$ and $80$ MeV, ${}^{208}\mathrm{Pb}$ at 100 and 70 MeV, and
${}^{93}\mathrm{Nb}$ and ${}^{181}\mathrm{Ta}$ (${}^{90}\mathrm{Zr}$ and
${}^{232}\mathrm{Th}$) at 100 MeV (70 MeV). The experimental data at 100 MeV
are taken from Ref. [29] and those at $E_{i}=80$ and $70$ MeV are taken from
Ref. [30]. In the experimental data at 80 and 70 MeV, the sharp increase at
very large $E_{f}$ is due to the elastic scattering events.

The contribution of multi-step processes becomes larger as the energy transfer
$\omega$ or $\theta$ increases, as discussed in $(p,nx)$ reaction in Ref.
[27]. For this reason, we focus only on the regions:
$\omega=E_{i}-E_{f}\lesssim 15$ MeV and $\theta\lesssim 60^{\circ}$ in Fig. 1,
where one-step process is expected to be dominant. The calculated DDXs with
LSCA reproduce the experimental data reasonably well, although in some cases
undershooting is found. It turns out that the undershooting is pronounced at
smaller $E_{i}$ or for heavier target nuclei and is particularly noticeable in
the DDXs on ${}^{58}\mathrm{Ni}$ at 100 MeV. Further investigation is needed
on this underestimation issue.

In the small energy-transfer region, $\omega\lesssim 15$ MeV, the DDXs with
AMA decrease rapidly as $\theta$ increases and show too strong angular
dependence compared to the data. In contrast to the AMA results, the DDXs with
LSCA show better agreements with the data in the angular dependence. At
$\theta\geq 30^{\circ}$, DDXs with LSCA are much closer to the experimental
data than those with AMA as discussed in Sec. III.3. These results imply that
the inclusion of the nuclear refraction, i.e., the change in the kinematics of
the deuteron inside the nucleus, is necessary to describe the angular
dependence of the experimental data.

###  III.3 The effect of the refraction

Below we discuss the angular dependence of the nuclear refraction on the DDX
of the ${}^{58}\mathrm{Ni}(d,d^{\prime}x)$ at 80 MeV. Figure 2 shows the
calculated DDXs as a function of the center-of-mass scattering angle
$\theta_{\mathrm{c.m.}}$ with $\omega=5$ MeV and the experimental data in Ref.
[30]. The solid (dashed) line represents the DDX with LSCA (AMA).

Figure 2: DDXs of the ${}^{58}\mathrm{Ni}(d,d^{\prime}x)$ at 80 MeV as a
function of the scattering angle with $\omega=5$ MeV. The solid (dashed) line
corresponds to the calculation with LSCA (AMA). The experimental data are
taken from Ref. [30].

It is shown that the DDX with LSCA gives moderate angular dependence and the
results at the middle and backward angles are improved, compared with those
with AMA. However, the gap between the result with LSCA and the experimental
data remains and the multi-step processes will be necessary to fill it.

By comparing the calculated DDXs using LSCA and AMA, one can see that there
are mainly two effects of the nuclear refraction. One is the extension of the
allowed region of $\theta_{\mathrm{c.m.}}$. For the DDX with AMA, the
($d,d^{\prime}x$) reaction is only allowed to about $42^{\circ}$. On the other
hand, the DDX with LSCA extends up to $\theta_{\mathrm{c.m.}}=73^{\circ}$ and
does not drop off at very small $\theta_{\mathrm{c.m.}}$. This is because
kinematics forbidden in AMA become allowed in LSCA by the refraction of the
momentum of the deuteron in the target nucleus. It should be noted that, as
one may find from Eq. (19), the ${\bm{R}}$ dependence of the kinematics of the
deuteron due to the refraction dictates the averaged local cross section of
the $d$-$N$ elementary process. In other words, whether the $d$-$N$ process
can take place or not depend on ${\bm{R}}$ through $\bm{q}({\bm{R}})$ and
$k_{F}({\bm{R}})$.

To see this more clearly, we show in Fig. 3(a) $\bm{q}(\bm{R})$ calculated
with LSCA, and in Figs. 3(b) and (c) the kinematically-allowed reaction
regions of the elementary process with LSCA and AMA, respectively, at
$\theta_{\mathrm{c.m.}}=60^{\circ}$ of the Fig. 2 case.

Figure 3: (a) The local momentum transfer with LSCA of
${}^{58}\mathrm{Ni}(d,d^{\prime}x)$ at 80 MeV with $\omega=5$ MeV at
$\theta_{\mathrm{c.m.}}=60^{\circ}$. (b) The reaction region with LSCA. The
color bar with the value of “1” means the region in which the $d$-$N$ process
is allowed, while with the value of “0” means the region in which the process
is not allowed. (c) Same as (b) but with AMA. (d), (e) and (f) same as (a),
(b) and (c), respectively, but at $\theta_{\mathrm{c.m.}}=15^{\circ}$.

In Figs. 3(b) and (c), a color bar with the value of “1” indicates the regions
where the $d$-$N$ processes are allowed while “0” indicates the regions where
the $d$-$N$ processes are not allowed. In AMA, which does not include the
nuclear refraction, $\bm{q}(\bm{R})$ is the same as the asymptotic momentum
transfer $\bm{q}$. Figure 3(c) shows that there are no kinematically-allowed
reaction regions; it is found that this is because $q$ is too large to allow
the elementary process. On the other hand, with LSCA, there is a region in
which the $d$-$N$ process is allowed because $\bm{q}(\bm{R})$ is dispersed by
the nuclear refraction and may have smaller values.

The other effect of refraction is the decrease in the DDX at forward angles.
This is because LSCA makes kinematically-allowed reaction regions narrower.
Figures 3(d), (e), and (f) are the same as Figs. 3(a), (b), and (c),
respectively, but at $\theta_{\mathrm{c.m.}}=15^{\circ}$. In Fig. 3(f), one
can see that the $d$-$N$ process is kinematically allowed in the very broad
region when AMA is used. On the other hand, in Fig. 3(e), the reaction region
with LSCA becomes narrower than that with AMA. This is because
$\bm{q}(\bm{R})$ is dispersed and can have too large values, as in Fig. 3(d),
to kinematically allow the $d$-$N$ process.

From these results, we conclude that the two effects of the nuclear refraction
can be understood as the changes in the kinematically-allowed reaction regions
associated with the dispersion of $\bm{q}(\bm{R})$. Figure 4 is the same as
Fig. 2 but with $\omega=30$ MeV. In Fig. 4, we can see that the two effects of
refraction remain even when $\omega$ is large. The angular distribution of the
DDX with LSCA gives a better agreement with the experimental data than that
with AMA, as in the case of $\omega=5$ MeV shown in Fig. 2.

Figure 4: Same as Fig. 2 but for $\omega=30$ MeV.

Figures 5(a) and (b) show the DDXs, with $\omega=5$ MeV, of
${}^{58}\mathrm{Ni}(d,d^{\prime}x)$ at 80 MeV (40 MeV per nucleon) and
${}^{58}\mathrm{Ni}(p,p^{\prime}x)$ at 40 MeV, respectively.

Figure 5: (a)Same as Fig. 2 but horizontal axis is the momentum transfer. (b)
Same as (a) but of ${}^{58}\mathrm{Ni}$$(p,p^{\prime}x)$.

One sees that the effects of the refraction are more significant in the
$(d,d^{\prime}x)$ reaction than in the $(p,p^{\prime}x)$ reaction. This is
mainly because the distorting potential between the deuteron and the target is
deeper than that between the proton and the target. Although the importance of
the nuclear refraction has been pointed out in the preceding studies of
$(p,p^{\prime}x)$ reactions with SCDW [22], its effect is found to be not very
significant. For $(d,d^{\prime}x)$, as shown in Fig. 1, the nuclear refraction
completely changes the behavior of the DDX. To analyze the $(d,d^{\prime}x)$
reaction data, inclusion of the nuclear refraction will be necessary.

##  IV SUMMARY

We have improved SCDW to the inclusive $(d,d^{\prime}x)$ reaction. The
calculated DDXs of the $(d,d^{\prime}x)$ were compared with the experimental
data of various targets at several deuteron emission angles. Except for some
cases, the calculated DDXs with LSCA reasonably reproduce the experimental
data in the regions where one-step process is expected to be dominant, i.e.,
for small energy transfer and at forward angles. On the other hand, the DDXs
with AMA, which does not include changes in the kinematics of the deuteron due
to the distorting potential, has too strong angular dependence and
considerably underestimate the experimental data at $\theta\geq 30^{\circ}$.
These results imply that the nuclear refraction effect on the $d$-$N$
elementary process is necessary to reproduce the experimental data.

We have shown two effects of the refraction by comparing the DDXs with LSCA
and AMA as a function of the scattering angle. One is the extension of the
kinematically-allowed scattering angles to the backward region. The other is
the decrease in the DDX at forward angles. Both effects can be understood by
the changes in the regions where the $d$-$N$ elementary processes are allowed.
It is confirmed that the refraction effect is more significant on the
$(d,d^{\prime}x)$ than on the $(p,p^{\prime}x)$ due to the strong distortion
effect on deuteron.

For better agreement with experimental data in the large energy-transfer
region, it will be necessary to modify the present SCDW model for multi-step
processes. Another future work will be to consider the deuteron breakup, which
is not explicitly treated in this study, to describe the inclusive $(d,nx)$
reaction that is important in nuclear data science.

###### Acknowledgements.

The authors thank Y. Chazono for providing us with the $d$-$N$ scattering
cross section. H.N. and K.O. thank Y. Watanabe for fruitful discussions. This
work has been supported in part by Grants-in-Aid of the Japan Society for the
Promotion of Science (Grants No. JP20K14475, No. JP21H00125, and No.
JP21H04975). The computation was carried out with the computer facilities at
the Research Center for Nuclear Physics, Osaka University.

*

##  Appendix A Validity of LSCA and AMA

The validity of LSCA in nucleon scattering has been examined in Refs. [24,
37]. In these papers, it is shown that LSCA works well for the propagation up
to about 1.5 fm at energies above 50 MeV. The validity of LSCA for the
$\alpha$ particle was also verified in Ref. [38].

Figure 6: The validity of LSCA and AMA. The real part of the exact
$\chi^{(+)}_{i}$ (solid line), with LSCA (dashed line), and with AMA (dotted
line) are compared. In Figs. (a) and (b), the propagation in radial direction
from $(6~{}\mathrm{fm},~{}120^{\circ},~{}0^{\circ})$ and
$(2~{}\mathrm{fm},~{}90^{\circ},~{}0^{\circ})$ in the spherical coordinate
representation are shown, respectively.

For the deuteron, however, its validity has not been confirmed. In Fig. 6, we
examine the validity of LSCA and AMA for the $d$-${}^{58}\mathrm{Ni}$
distorted wave $\chi^{(+)}_{i}$ at 50 MeV per nucleon. Fig. 6 shows the
propagation in the radial direction from (a)
$\bm{R}_{a}\equiv(6~{}\mathrm{fm},~{}120^{\circ},~{}0^{\circ})$ and (b)
$\bm{R}_{b}\equiv(2~{}\mathrm{fm},~{}90^{\circ},~{}0^{\circ})$ in the
spherical coordinate representation. The solid, dashed, and dotted lines show,
respectively, the real part of the exact wave function, that with LSCA, and
that with AMA. In Fig. 6(a), both approximations reproduce well the
propagation up to about 0.7 fm. It should be noted that the range of the
interaction between the deuteron and the nucleon is about 2.2 fm, and from the
factor $1/(A_{d}+1)=1/3$ for $\bm{s}$ in Eq. (7), LSCA and AMA are required to
be valid for the propagation up to about 0.7 fm. In Fig. 6(b), on the other
hand, while LSCA reproduces the propagation of the wave function well, AMA
does not. This is because the direction of the propagation direction $\bm{s}$
from $\bm{R}_{b}$ is orthogonal to the asymptotic momentum $\bm{k}_{c}$ of the
deuteron, i.e., $\bm{k}_{c}\cdot\bm{s}=0$ in Eq. (12). These results show that
the kinematics of the deuteron at $\bm{R}_{b}$ are significantly different
from the asymptotic ones due to the distorting potential, thus LSCA is
essential to trace the deuteron momentum inside the target nucleus.

## References

  * Butler [1950] S. T. Butler, Phys. Rev. 80, 1095 (1950).
  * Timofeyuk and Johnson [2020] N. Timofeyuk and R. Johnson, Progress in Particle and Nuclear Physics 111, 103738 (2020).
  * Kamimura  _et al._ [1986] M. Kamimura, M. Yahiro, Y. Iseri, Y. Sakuragi, H. Kameyama, and M. Kawai, Progress of Theoretical Physics Supplement 89, 1 (1986).
  * Austern  _et al._ [1987] N. Austern, Y. Iseri, M. Kamimura, M. Kawai, G. Rawitscher, and M. Yahiro, Physics Reports 154, 125 (1987).
  * Deltuva and Fonseca [2009] A. Deltuva and A. C. Fonseca, Phys. Rev. C 79, 014606 (2009).
  * Upadhyay  _et al._ [2012] N. J. Upadhyay, A. Deltuva, and F. M. Nunes, Phys. Rev. C 85, 054621 (2012).
  * Ogata and Yoshida [2016] K. Ogata and K. Yoshida, Phys. Rev. C 94, 051603(R) (2016).
  * Potel  _et al._ [2017] G. Potel, G. Perdikakis, B. V. Carlson, M. C. Atkinson, W. H. Dickhoff, J. E. Escher, M. S. Hussein, J. Lei, W. Li, A. O. Macchiavelli, A. M. Moro, F. M. Nunes, S. D. Pain, and J. Rotureau, The European Physical Journal A 53, 178 (2017).
  * Moeslang  _et al._ [2006] A. Moeslang, V. Heinzel, H. Matsui, and M. Sugimoto, Fusion Engineering and Design 81, 863 (2006), proceedings of the Seventh International Symposium on Fusion Nuclear Technology.
  * Nakayama  _et al._ [2016] S. Nakayama, H. Kouno, Y. Watanabe, O. Iwamoto, and K. Ogata, Phys. Rev. C 94, 014618 (2016).
  * Nakayama  _et al._ [2020] S. Nakayama, O. Iwamoto, and Y. Watanabe, EPJ Web Conf. 239, 03014 (2020).
  * Nakayama  _et al._ [2021] S. Nakayama, O. Iwamoto, Y. Watanabe, and K. Ogata, Journal of Nuclear Science and Technology 58, 805 (2021).
  * Glauber [1959] R. J. Glauber, in Lectures in Theoretical Physics,Vol. 1, p. 315 (Interscience, New York, 1959).
  * Hussein and McVoy [1985] M. Hussein and K. McVoy, Nuclear Physics A 445, 124 (1985).
  * Hencken  _et al._ [1996] K. Hencken, G. Bertsch, and H. Esbensen, Phys. Rev. C 54, 3043 (1996).
  * Ichimura  _et al._ [1985] M. Ichimura, N. Austern, and C. M. Vincent, Phys. Rev. C 32, 431 (1985).
  * Lei and Moro [2015a] J. Lei and A. M. Moro, Phys. Rev. C 92, 044616 (2015a).
  * Lei and Moro [2015b] J. Lei and A. M. Moro, Phys. Rev. C 92, 061602(R) (2015b).
  * Feshbach  _et al._ [1980] H. Feshbach, A. Kerman, and S. Koonin, Annals of Physics 125, 429 (1980).
  * Nishioka  _et al._ [1988] H. Nishioka, H. Weidenmüller, and S. Yoshida, Annals of Physics 183, 166 (1988).
  * Tamura  _et al._ [1982] T. Tamura, T. Udagawa, and H. Lenske, Phys. Rev. C 26, 379 (1982).
  * Luo and Kawai [1991] Y. L. Luo and M. Kawai, Phys. Rev. C 43, 2367 (1991).
  * Kawai and Weidenmüller [1992] M. Kawai and H. A. Weidenmüller, Phys. Rev. C 45, 1856 (1992).
  * Watanabe  _et al._ [1999] Y. Watanabe, R. Kuwata, S. Weili, M. Higashi, H. Shinohara, M. Kohno, K. Ogata, and M. Kawai, Phys. Rev. C 59, 2136 (1999).
  * Ogata  _et al._ [1999] K. Ogata, M. Kawai, Y. Watanabe, S. Weili, and M. Kohno, Phys. Rev. C 60, 054605 (1999).
  * Weili  _et al._ [1999] S. Weili, Y. Watanabe, M. Kohno, K. Ogata, and M. Kawai, Phys. Rev. C 60, 064605 (1999).
  * Ogata  _et al._ [2002] K. Ogata, Y. Watanabe, S. Weili, M. Kohno, and M. Kawai, Nuclear Physics A 703, 152 (2002).
  * Watanabe and Kawai [1993] Y. Watanabe and M. Kawai, Nuclear Physics A 560, 43 (1993).
  * Förtsch  _et al._ [2002] S. Förtsch, D. Ridikas, W. Mittig, H. Savajols, P. Roussel–Chomaz, G. Steyn, and J. Lawrie, Journal of Nuclear Science and Technology 39, 792 (2002).
  * Wu  _et al._ [1979] J. R. Wu, C. C. Chang, and H. D. Holmgren, Phys. Rev. C 19, 370 (1979).
  * An and Cai [2006] H. An and C. Cai, Phys. Rev. C 73, 054605 (2006).
  * Perey and Buck [1962] F. Perey and B. Buck, Nuclear Physics 32, 353 (1962).
  * [33] TWOFNR, Usermanual, https://people.nscl.msu.edu/ brown/reaction-codes/twofnr/twofnr.pdf.
  * Chazono  _et al._ [2022] Y. Chazono, K. Yoshida, and K. Ogata, Phys. Rev. C 106, 064613 (2022).
  * Love and Franey [1981] W. G. Love and M. A. Franey, Phys. Rev. C 24, 1073 (1981).
  * Franey and Love [1985] M. A. Franey and W. G. Love, Phys. Rev. C 31, 488 (1985).
  * Minomo  _et al._ [2010] K. Minomo, K. Ogata, M. Kohno, Y. R. Shimizu, and M. Yahiro, Journal of Physics G: Nuclear and Particle Physics 37, 085011 (2010).
  * Yoshida  _et al._ [2016] K. Yoshida, K. Minomo, and K. Ogata, Phys. Rev. C 94, 044604 (2016).


# Lagrangian Neural Network with Differentiable Symmetries and Relational
Inductive Bias

Ravinder Bhattoo
Department of Civil Engineering
Indian Institute of Technology Delhi
<EMAIL_ADDRESS>
& Sayan Ranu
Department of Computer Science and Engineering
Indian Institute of Technology Delhi
<EMAIL_ADDRESS>
& N. M. Anoop Krishnan
Department of Civil Engineering
Indian Institute of Technology Delhi
<EMAIL_ADDRESS>

###### Abstract

Realistic models of physical world rely on differentiable symmetries that, in
turn, correspond to conservation laws. Recent works on Lagrangian and
Hamiltonian neural networks show that the underlying symmetries of a system
can be easily learned by a neural network when provided with an appropriate
inductive bias. However, these models still suffer from issues such as
inability to generalize to arbitrary system sizes, poor interpretability, and
most importantly, inability to learn translational and rotational symmetries,
which lead to the conservation laws of linear and angular momentum,
respectively. Here, we present a momentum conserving Lagrangian neural network
(MCLnn) that learns the Lagrangian of a system, while also preserving the
translational and rotational symmetries. We test our approach on linear and
non-linear spring systems, and a gravitational system, demonstrating the
energy and momentum conservation. We also show that the model developed can
generalize to systems of any arbitrary size. Finally, we discuss the
interpretability of the MCLnn, which directly provides physical insights into
the interactions of multi-particle systems.

_Keywords_ Machine Learning, Dynamical Systems

##  1 Introduction and Related Work

Realistic modeling of the time-evolution of multi-particle physical systems
lie at the heart of several fields of science and engineering [1, 2, 3, 4].
Examples include multiple bodies connected by springs, bodies interacting
under gravitational forces, or even systems at the atomic and mesoscale such
as proteins or colloidal gels, which range several orders ($\approx 10^{20}$)
of length and timescales from atomic systems to planetary systems [2, 5].
Traditionally, the trajectory of these systems are obtained by solving the
associated differential equations numerically. These differential equations,
in turn, are derived based on the invariant or conserved quantities of a
system. Recent studies on Lagrangian (Lnn) and Hamiltonian neural networks
show that one of these invariant quantities, energy, can be learned directly
from the data enabling realistic simulation of systems [6, 7, 4, 8, 9, 10,
11]. However, these approaches fail to model multi-particle interactions
effectively due to the following limitations.

$\bullet$ Conservation of momentum: Invoking Noether’s theorem [12], “a
differentiable symmetry in action results in a conserved quantity”.
Conversely, “every conserved quantity can be derived from an underlying
symmetry in action”. It can be shown that an interacting multi-particle
system, that is closed, respects three conservation laws, namely, energy,
linear momentum and angular momentum [12, 1]. While the first is a consequence
of symmetry of the system with respect to time, the second and third is a
consequence of symmetry with respect to translation and rotation in space,
respectively. Thus, the Lagrangian of such a system remains invariant under
the translation and rotation operations on the Cartesian coordinates. For
example, three balls connected by a spring will have the same Lagrangian (and
hence the interaction forces) when the system as a whole is subjected to
rotation or translation. That is, as long as the relative positions of the
balls are preserved, translation or rotation by any magnitude will not have an
effect on the system dynamics. While Lnns have been demonstrated to learn the
symmetry in time by conserving the energy [6], they have not been designed to
respect the other two conservation laws with respect to momentum. Learning
these symmetries directly from the data is an extremely hard problem and may
require large amounts of data. Furthermore, as we will show in Sec. 4, even
when trained on large volumes of training data, the Lnns fail to generalize
well on unseen data. This questions the realistic nature of the dynamics
simulated by Lnns.

$\bullet$ Generalizability to unseen system sizes: Lnns lack the ability to
generalize to system sizes beyond the training set. For instance, an Lnn
trained on three balls connected by a spring cannot model a system of four
balls connected by a spring despite having the same conserved quantities and
interactions. This limits the broader applicability of Lnns to realistic
physical systems.

$\bullet$ Interpretability: Lnn [6] directly learns (and predicts) the
Lagrangian of a system as a function of the position and velocities of all the
particles in the system. Owing to this design, it is not possible to recover
the inter-particle dynamics within a multi-particle system in the form of
pairwise potential energy or forces. This limits the interpretability of the
Lnn as the contribution of positions and velocities towards the total
Lagrangian of the system is represented using a black-box deep neural network.

In this paper, we address the above limitations. Specifically, our
contributions are as follows:

  * •

Problem Formulation and Architecture Design: We reformulate Lnn [6] with a
relational inductive bias by applying a transformation on the Cartesian
coordinates of the system and decoupling the terms of the Lagrangian. Armed
with this reformulation, we develop a momentum conserving Lnn called MCLnn.
MCLnn introduces several key innovations. (1) First, MCLnn preserves the
rotational and translation symmetries. (2) Second, MCLnn generalizes to unseen
system sizes. (3) Third, MCLnn generates interpretable models with the output
characterizing pair-wise interactions of the bodies in the system, and thereby
providing physical insights into the system dynamics.

  * •

Theoretical Characterization: We rigorously establish that MCLnn conserves
energy, linear and angular momentum. Furthermore, in contrast to Lnn [6], we
show that MCLnn can directly learn trajectories (positions) of systems without
requiring training on acceleration. This is a desirable property. In many
cases including experimental systems, access to acceleration of each entity at
every time instant may not be possible. While position of particles is a
direct observable, acceleration is a derived quantity.

  * •

Empirical Evaluation: We perform in-depth evaluation across three multi-
particle systems namely, balls connected by linear and non-linear springs, and
gravitational system. Our experiments establish that MCLnn performs
significantly better than Lnn with long-term stability, generalizes to
arbitrary-sized multi-particle systems, and is capable of learning interaction
dynamics with $3$-$4$ orders of magnitude lower volume of training data.

##  2 Preliminaries and Problem Formulation

In this section, we introduce the preliminary concepts central to our work and
formulate the problem. As notational convention, vectors are represented with
an overhead arrow (Ex. $\vec{v}$) and higher-order tensors, such as matrices,
in bold. A summary notations used is provided in the appendix A.7.

######  Definition 1 (Multi-particle System)

An $\mathcal{N}$-body multi-particle system contains $\mathcal{N}$ particles
$P=\\{n_{1},\cdots,n_{\mathcal{N}}\\}$. At any given time $t$, particle
$n_{i}$ is characterized by its position $\vec{q}_{i}^{\,t}$ and velocity
$\dot{\vec{q}}_{i}^{\,t}$.

The positions of particles in an $\mathcal{N}$-body system are not static.
They change due to various interaction forces at play (Ex: a gravitational
system). These positional changes are captured in the form of trajectories.

######  Definition 2 (Multi-particle Trajectory)

The trajectory of an $\mathcal{N}$-body multi-particle system $P$ over a time
horizon $T=[t_{s},t_{e}]$ is the sequence of position and velocity vectors
$({\mathbf{q}}^{\,t},{\mathbf{\dot{q}}}^{\,t})\mid t\in T\\}$, where
$\mathbf{{q}}^{\,t}=\\{\vec{q}_{i}^{\,t}\mid n_{i}\in P\\}$ and
$\dot{\mathbf{q}}^{t}=\\{\dot{\vec{q}}_{i}^{\,t}\mid n_{i}\in P\\}$.

In physics, the time-evolution or trajectory of interacting particles is
obtained by solving the differential equations of motion. The numerical
solution to these equations provide the acceleration of the particle, which
can then be used to obtain the updated velocity and positions. Indeed, the
acceleration of particles can be directly learned by neural networks as a
function of its position and velocities. However, it fails to learn the
underlying symmetries, which in turn leads to the laws of conservation of
energy and momenta.

Recent approaches to predict trajectories by learning the Lagrangian through
an Lnn [6, 7, 9, 11] has been shown to be effective in learning the symmetry
in time [6, 9]. Lagrangian, defined as
$L(\mathbf{{q}},\mathbf{\dot{q}})=T(\mathbf{\dot{q}})-V(\mathbf{q})$, is a
scalar functional of kinetic energy ($T(\mathbf{\dot{q}})$) and potential
energy ($V(\mathbf{{q}})$), that depends on both the set of positions
$\mathbf{q}$, and velocities $\mathbf{\dot{q}}$ of the particles in the system
[1]. Lnn bypasses the necessity to learn accelerations for each particle and
directly predicts the scalar Lagrangian for a system. Additionally, the
trajectory predicted using Euler-Lagrange (EL) equation conserves energy
resulting in an overall better prediction of the system dynamics.

Theoretically, although learning the Lagrangian of a system is enough to
predict its dynamics, the Lagrangian itself exhibits some symmetries that are
not imposed by the EL equation. Lagrangian of a system exhibits translational
and rotational symmetry. Specifically, the Lagrangian of a system
$L(\mathbf{q},\mathbf{\dot{q}})$ and
$L(\mathbf{q^{\prime}},\mathbf{\dot{q}})$, where
$\mathbf{q^{\prime}}=\\{\vec{q}_{i}+\vec{\epsilon}\mid n_{i}\in P\\}$ for any
finite value of $\vec{\epsilon}$ are equivalent. Similarly,
$L(\mathbf{q},\mathbf{\dot{q}})$ and
$L(\mathbf{Q}\mathbf{q},\mathbf{\dot{q}})$, where $\mathbf{Q}$ is an
orthogonal tensor representing a rotation, are equivalent. These symmetries
are not provided as an inductive bias in the Lnn framework.

Further, Lnns are trained by giving the coordinates as input and minimizing
the loss between the actual acceleration and acceleration predicted by Lnn
through EL equations. This requires a priori access to the actual
acceleration, which may not be available in many cases, for example,
experimental trajectory of colloidal gels visualized by video camera.

Here, we aim to address these open challenges by reformulating the Lnn with a
relational inductive bias, resulting in Momentum Conserving Lnn (MCLnn).
Specifically, we aim to learn the dynamics of multi-particle interacting
systems purely from the trajectory of the constituent particles.

######  Problem 1 (Momentum Conserving Lnn for Trajectory Prediction)

Let $\mathbb{T}$ be a set of trajectories of $\mathcal{N}$-body systems.
Furthermore, let there be a hidden joint distribution of physics-constrained
configurational and temporal space from which $\mathbb{T}$ has been sampled.
Our goal is to learn this hidden distribution. Towards that end, we want to
learn:

  1. i.

A generative model $p(\mathbb{T})$ that maximizes the likelihood of generating
$\mathbb{T}$.

  2. ii.

The generative model must respect the laws of physics such as conservation of
energy and momenta.

Once learned, this generative model can be used to predict trajectories of
unseen $\mathcal{N}$-body systems.

Figure 1: MCLnn framework.

##  3 Theory

The flowchart of MCLnn is shown in Fig. 1. In this approach, the positions and
velocities at any time $t$ is used to predict the trajectory of the system,
while respecting the conservation laws of energy and momenta. We analyze the
performance of MCLnn on three systems, namely, linear spring, non-linear
spring, and gravitational system to predict the accurate dynamics of these
multi-particle interacting systems.

###  3.1 Euler-Lagrange Equation

Consider a system of $\mathcal{N}$-particles interacting with each other under
a potential. Let us define functional, namely, “action”, $S$ as:

$S=\int_{t_{0}}^{t_{1}}L\ dt$ (1)

where $L$ is the Lagrangian of the system. Then trajectory taken by this
system to move from the position $\mathbf{q^{0}}$ to $\mathbf{q^{1}}$ in time
$t_{0}$ to $t_{1}$ is the one that makes the “action” $S$ stationary. This
leads to the equation of motion governing the dynamics of system, namely, the
EL equation as:

$\frac{d}{dt}\frac{\partial L}{\partial\mathbf{\dot{q}}}=\frac{\partial
L}{\partial\mathbf{q}}$ (2)

The acceleration of particles in the system,
$\mathbf{\ddot{q}}=\\{\ddot{\vec{q}}_{i}\mid n_{i}\in P\\}$ can be computed
directly from the EL equation (see Appendix A.3) as:

$\mathbf{\ddot{q}}=\left(\nabla_{\mathbf{\dot{q}}\mathbf{\dot{q}}}L\right)^{-1}\left[\nabla_{\mathbf{q}}L-\left(\nabla_{\mathbf{\dot{q}}\mathbf{{q}}}L\right)\mathbf{\dot{q}}\right]$
(3)

This acceleration can then be used to compute the updated positions and
velocities of particles. Note that EL equation leads to the conservation of
Hamiltonian, defined as:
$H(\mathbf{q},\mathbf{\dot{q}})=T(\mathbf{\dot{q}})+V(\mathbf{q})$, which also
represents the total energy of the system.

###  3.2 Momentum Conserving Lagrangian Neural Network

As mentioned earlier, Lagrangians exhibit translational and rotational
symmetry. Learning this directly from the positions of particles is extremely
challenging. To address this challenge, we reformulate the Lnn. Consider five
balls connected by linear springs. We will use this running example to explain
MCLnn. The input for the MCLnn is the position and velocities at any time $t$.
Further, we identify the set of all the pairs of particles in the system as
($i,j$). For a given set, we compute the difference between the positions
followed by computing the $l^{2}$ norm to obtain the pairwise distances
$q_{ij}$, which corresponds to the length of the spring. Thus,

$q_{ij}=\sqrt{(\vec{q}_{i}-\vec{q}_{j})\cdot(\vec{q}_{i}-\vec{q}_{j})}$ (4)

The $q_{ij}$ is given as the input to a neural network which outputs a scalar
$V_{ij}$. Note that the $V_{ij}$ may be considered to be the pair-wise
potential energy of a spring connecting two particles $i$ and $j$. $V_{ij}$
summed over all the pairs of particles gives the $V(\mathbf{q})$, the total
potential energy of the system. For multi-particle systems, the kinetic energy
of a particle $n_{i}$ is the function of its velocity $\dot{\vec{q}}_{i}$, can
defined as [1, 9]:

$T_{i}(\dot{\vec{q}}_{i})=\frac{1}{2}m_{i}\dot{\vec{q}}_{i}^{\,2}.$ (5)

Here, we invoke this expression of kinetic energy, which when summed over all
the particles provide the total kinetic energy, $T$. Then, $T$ and $V$ are
used to obtain the $L$ as in Eq. 7. The $L$ when substituted into the EL
equation (Eq. 2), provides the acceleration and the updated position and
velocities of the system using the velocity-Verlet integration [13]. The loss
function of MCLnn is on the predicted and actual positions at time
$t+1,t+2,\ldots,t+n$, which is then back-propagated to train the neural
network. Specifically, the loss function is as follows.

$\mathcal{L}=\frac{1}{|\mathbb{T}|}\sum_{\forall\mathcal{T}\in\mathbb{T}}\left(\frac{1}{|\mathcal{T}|}\sum_{t=2}^{|\mathcal{T}|}\frac{1}{\mathcal{N}}\left(\sum_{i=1}^{\mathcal{N}}\left(\mathcal{T}.q_{i}^{\,t}-p\left(\mathcal{T}.q_{i}^{\,t}\right)\right)^{2}\right)\right)$
(6)

Here, $p(\mathcal{T}.q_{i}^{\,t})$ is the predicted position for the $i^{th}$
particle in $\mathcal{T}$ at time $t$ and $\mathcal{T}.q_{i}^{\,t}$ is the
true position. $|\mathcal{T}|$ denotes the last time step in trajectory
$\mathcal{T}$. The updated expression for Lagrangian for an
$\mathcal{N}$-particle system can be written as:

$L(\mathbf{q},\mathbf{\dot{q}})=\sum_{i=1}^{n}\frac{1}{2}m_{i}\dot{\vec{q}}_{i}^{\,2}-\sum_{i=1}^{n-1}\sum_{j=i+1}^{n}V_{ij}(q_{ij})$
(7)

It should be noted that the specific form of kinetic energy $T$ does not limit
the applicability of the present approach [9]. Since, the focus of the work is
to model multi-particle systems and it is well-known that the kinetic energy
of the system generally follows Eq. 5, we choose this expression to simplify
the problem. In cases where the expression for kinetic energy is not known or
is different, the function can be replaced by an additional neural network
with an approach similar to that for $V$.

###  3.3 Translational and rotational symmetry of MCLnn

###### Theorem 1

MCLnn exhibits translational and rotational symmetry.

Proof. Consider, the Lagrangian of the system of $\mathcal{N}$-particles as
$L(\mathbf{q},\mathbf{\dot{q}})=T(\mathbf{\dot{q}})-V(\mathbf{q})$. First, we
focus on the translational symmetry and next on rotational symmetry.

###### Lemma 1

MCLnn exhibits translational symmetry.

Proof. When the particles are subjected to a translation $\vec{\epsilon}$, the
updated positions are given by
$\mathbf{q^{\prime}}==\\{\vec{q}_{i}+\vec{\epsilon}\mid n_{i}\in P\\}$. The
updated Lagrangian of the system is given by
$L(\mathbf{q^{\prime}},\mathbf{\dot{q}})$. To prove
$L(\mathbf{q^{\prime}},\mathbf{\dot{q}})=L(\mathbf{q},\mathbf{\dot{q}})$, we
need to prove $V(\mathbf{q^{\prime}})=V(\mathbf{q})$ as the other term in the
Lagrangian remains unaffected. Applying the transformation on the positions
and computing the $l^{2}$ norm,

$\displaystyle q_{ij}^{\prime}=$
$\displaystyle\sqrt{[(\vec{q}_{i}+\vec{\epsilon})-(\vec{q}_{j}+\vec{\epsilon})]\cdot[(\vec{q}_{i}+\vec{\epsilon})-(\vec{q}_{j}+\vec{\epsilon})])}$
$\displaystyle=$
$\displaystyle\sqrt{(\vec{q}_{i}-\vec{q}_{j})\cdot(\vec{q}_{i}-\vec{q}_{j})}$
$\displaystyle=$ $\displaystyle q_{ij}$ (8)

which implies $V(\mathbf{q^{\prime}})=V(\mathbf{q})$. Therefore,
$L(\mathbf{q^{\prime}},\mathbf{\dot{q}})=L(\mathbf{q},\mathbf{\dot{q}})$.

###### Lemma 2

MCLnn exhibits rotational symmetry.

Proof. When the system is subjected to pure rotation, the positions are
transformed as $\mathbf{Q}\mathbf{q}$, where $\mathbf{Q}$ is an orthogonal
tensor representing a rotation. Correspondingly, the Lagrangian is modified as
$L(\mathbf{Q}\mathbf{q},\mathbf{\dot{q}})$. As in the case of translation, we
need only to prove that $V(\mathbf{Q}\mathbf{q})=V(\mathbf{q})$. To prove
this, it is worth recalling that, for an orthogonal tensor $\mathbf{Q}$,
$\mathbf{Q}\mathbf{Q}^{\,t}=\mathbf{Q}^{\,t}\mathbf{Q}=I$ . Now, the proof can
be obtained by computing $q_{ij}^{\prime}$ as:

$\displaystyle q_{ij}^{\prime}$
$\displaystyle=\sqrt{(\mathbf{Q}\vec{q_{i}}-\mathbf{Q}\vec{q_{j}})\cdot(\mathbf{Q}\vec{q_{i}}-\mathbf{Q}\vec{q_{j}})}$
$\displaystyle=\sqrt{\mathbf{Q}(\vec{q_{i}}-\vec{q_{j}})\cdot\mathbf{Q}(\vec{q_{i}}-\vec{q_{j}})}$
$\displaystyle=\sqrt{\mathbf{Q}\mathbf{Q}^{\,T}(\vec{q_{i}}-\vec{q_{j}})\cdot(\vec{q_{i}}-\vec{q_{j}})}\hskip
28.45274pt$
$\displaystyle=\sqrt{(\vec{q}_{i}-\vec{q}_{j})\cdot(\vec{q}_{i}-\vec{q}_{j})}$
$\displaystyle=q_{ij}$ (9)

This implies that $V(\mathbf{Q}\mathbf{q})=V(\mathbf{q})$ and hence,
$L(\mathbf{q},\mathbf{\dot{q}})=L(\mathbf{Q}\mathbf{q},\mathbf{\dot{q}})$.
Thus, we demonstrate that the Lagrangian in MCLnn exhibits both translational
and rotational symmetry. Invoking Noether’s theorem, it can be proven that the
Lagrangian preserving these symmetries will consequently conserve of linear
and angular momenta, respectively. (For proof, see App. A.4 and A.5).

Combining Lemma 1 and Lemma 2, we get that the MCLnn will respect the laws of
conservation energy, linear momentum, and angular momentum. $\hfill\square$

###  3.4 Generalizability and interpretability

The MCLnn exhibits a granular structure, which applies the computations on
individual or pair-wise entities that are aggregated to obtain the total
Lagrangian of the system. This structure, by design, allows model
generalizability to unseen system sizes. Once trained on an $n$-body system,
the MCLnn can be applied on an $m$-body system ($m\neq n$) directly by
replacing $n$ with $m$ in Eq. 7.

Finally, we discuss interpretability of MCLnn. The main function learned by
MCLnn is the $V_{ij}$, which depends on $q_{ij}$. Since, $V_{ij}$ adds to form
the total potential energy $V(\textbf{q})$ of the system, it is reasonable to
assume that $V_{ij}$ represents the pair-wise potential energy. As we
demonstrate later in Sec. 4.4, we show that $V_{ij}$ indeed corresponds to the
pair-wise potential energy of the ground truth.

##  4 Experiments

In this section, we benchmark MCLnn and establish:

  * •

Accuracy: MCLnn is more accurate in modeling trajectories of
$\mathcal{N}$-body systems than Lnn [6].

  * •

Conservation of Physics: Consistent with the theoretical analysis in Sec. 3.3,
MCLnn preserves the laws of physics more comprehensively than Lnn

  * •

Generalizability and Interpretability: MCLnn generalizes to unseen data better
than Lnn and enables higher interpretability due to predicting the potential
energy instead of the Lagrangian directly.

###  4.1 Experimental Setup

We use JAX [14], a high-performance numerical computing python package and JAX
MD [15], a JAX based molecular dynamics package for the $n$-body simulations.

####  4.1.1 Tasks

To benchmark the performance of MCLnn, we conduct experiments on three multi-
body systems, namely, linear and non-linear springs, and gravitational
systems.

$\bullet$ Linear spring. We simulate the dynamics of three particles (e.g.:
balls with volume tending to zero) connected with each other by linear
springs. The initial conditions of the simulations include the positions and
velocities of the three particles and the stiffness of the spring. The
potential energy of two balls connected by a spring is given by
$V_{ij}=\frac{1}{2}k(q_{ij}-q_{0})^{2}$, where $q_{ij}$ is the instantaneous
distance between the springs and $q_{0}$ is the equilibrium distance, and $k$
is the stiffness of the spring. The force between connected pair of particles
is computed as $F_{ij}=-\frac{dV_{ij}}{dq_{ij}}=-k(q_{ij}-q_{0})$. Note that
the value of $k$ is kept constant for all the springs.

$\bullet$ Non-linear spring. We simulate the dynamics of three particles
connected by non-linear springs. The initial conditions of the simulations
include the positions and velocities of the three particles. The potential
energy of two particles connected by a non-linear spring is given by
$V_{ij}=\frac{1}{2}k(q_{ij}-q_{0})^{4}$, where $q_{ij}$ is the instantaneous
distance between the springs and $q_{0}$ is the equilibrium distance, and $k$
is the stiffness of the spring. The force between connected pair of particles
is computed as $F_{ij}=-2k(q_{ij}-q_{0})^{3}$. The value of $k$ is kept
constant for all the springs.

$\bullet$ Gravitational system. We simulate the dynamics of four particles
interacting with each other under the gravitational force. The initial
conditions of the simulations include the positions and velocities of the four
masses. The potential energy of a pair of particles with mass $m_{1}$ and
$m_{2}$ if given by $V_{ij}=-\frac{Gm_{1}m_{2}}{q_{ij}}$, where $G$ is
universal gravitational constant and $q_{ij}$ is the instantaneous distance
between the pair of particles. The force between the pair of particles is
computed as $F_{ij}=\frac{Gm_{1}m_{2}}{q_{ij}^{2}}$. Note that the mass $m$ of
all the particles are kept constant.

Data Generation: We use forward simulation of $\mathcal{N}$ particles to
generate training data corresponding to each task. First, an
$\mathcal{N}$-body configuration with given initial positions and velocities
are considered. The potential energy of the system is computed from the
positions using analytical expressions associated with each task as described
above. The force on each particle is then calculated using as the gradient of
potential energy, which is used to compute the acceleration. Finally, the
velocity-Verlet algorithm is used for time integration to obtain the updated
positions of the particles. Please see the appendix (A.6) for configurations
(including the number of particles, initial configuration, and time step) of
each of the tasks.

####  4.1.2 Baseline

The original Lnn takes both the positions and velocities of all the particles
as input and gives the total Lagrangian as the output. In the present case, we
make the assumption on the functional form of kinetic energy (see Eq.5). As
such, to ensure a fair comparison, a slightly modified version of the original
Lnn is considered here as the baseline. For this baseline, the positions are
given as input to the neural network which gives the potential energy
$V(\vec{q})$ as output. The kinetic energy is computed using the functional
form, from which the Lagrangian is computed. All the remaining training
procedures are maintained same as in the original Lnn. It should be noted that
decoupling the $T$ and $V$ terms and providing the expression for $T$ should
make the learning easier for the Lnn and hence should give improved
performance in comparison to the original Lnn. This model will be referred to
as baseline or baseline Lnn, henceforth.

####  4.1.3 Metrics

To analyze the performance of the MCLnn in comparison to the baseline, the
choice of the appropriate metrics is crucial. Since the systems considered
here are chaotic, slight differences in the initial conditions or predictions
will lead to differences in the trajectories. More importantly, given the
current state of a system, multiple “correct” trajectories may exist as long
as they represent the degenerate or equivalent configurational states of the
system. As such, purely computing the error in trajectory with respect to the
ground truth may not be representative of the performance of the model.
Learning and predicting the Lagrangian ensures that the state of the system is
represented accurately, resulting in the prediction of the same or equivalent
trajectories of the system. Further, the trajectories should respect the
spatial and temporal symmetries resulting in conservation of Hamiltonian
(total energy), linear, and angular momenta. Thus, we focus on the Mean
Absolute Error (MAE) between predicted (1) Lagrangian, (2) Hamiltonian (total
energy), (3) linear and (4) angular momenta, with respect to the ground truth.

####  4.1.4 Training

For training MCLnn, 100 different trajectories, each having 20 points with
fixed time intervals, are used for all three tasks. Each trajectory has
different initial conditions and hence represent different configurations that
are accessible to the system. The model performance is evaluated by comparing
the predicted trajectory with respect to the ground truth for the same initial
state.

For training Lnn, we sample 10,000 data points (set of
($\mathbf{q},\mathbf{\dot{q}})$ as input) and $\mathbf{\ddot{q}}$ as output)
from forward simulations. The training is performed using the mean squared
error between the acceleration predicted by the EL equation and the ground
truth. The detailed configurations and parameters associated with simulation
and training of each task is provided in App. A.6.

####  4.1.5 Inference

To asses the long term stability of models, we perform forward simulations
using trained Lnn and MCLnn models for a prolonged time interval for each
task. The initial condition and other parameters are kept same for both Lnn
and MCLnn. Furthermore, the initial conditions of the trajectories during
inference are different from those encountered during the training.
Consequently, all trajectories during evaluation are unseen. We record
Lagrangian, Hamiltonian, linear momentum and angular momentum for the whole
trajectory. The predicted forward simulation is compared with the ground truth
data as generated in Sec. 4.1.1.

Figure 2: Results of the forward simulation of linear spring system from the
same initial configuration using ground truth (continuous line), baseline Lnn
(dashed and dotted line) and MCLnn (dashed lines). (a) The evolution of
Lagrangian, Hamiltonian (total energy), linear and momentum for all the
systems. (b) The error in Lagrangian, Hamiltonian (total energy), linear and
momentum predicted by the baseline Lnn and MCLnn with respect to the ground
truth.

###  4.2 Accuracy

First, we focus on the linear spring system with three particles. Figure
2(a.1)-(a.4) shows the Lagrangian, Hamiltonian (total energy), linear, and
angular momenta of the system. The error in these quantities predicted by the
baseline Lnn and MCLnn are shown in Figure 2(b.1)-(b.4). We observe that the
Lagrangian predicted by the baseline Lnn starts diverging after 30 time steps
of the forward simulation. We notice that the divergence in the Lagrangian is
accompanied by a divergence in the linear momentum as well. Up on further
simulation, the angular momentum also starts diverging. In contrast, the
Lagrangian predicted by the MCLnn follows the ground truth with very little
error. Further, the Hamiltonian, and linear and angular momenta remain
conserved in the MCLnn. This suggests that the trajectory predicted by the
MCLnn is accurate and stable with no long term drift or divergence. It is
worth noting that, for the baseline Lnn, despite the poor predictions of the
Lagrangian at higher values of time steps, the total energy remain conserved.
This points to the fact that the energy conservation is a constraint enforced
by the EL equation and is not a reflection on the quality Lagrangian. Any
value of Lagrangian can still yield a constant Hamiltonian (energy), provided
the constraints as imposed by the EL equations on the Lagrangian are
satisfied. As such, Lnns should be evaluated considering multiple metrics as
demonstrated here.

Now, we focus on the non-linear spring and the gravitational system. Figure
3(a) and (c) shows the Lagrangian, Hamiltonian (total energy), linear, and
angular momenta of the non-linear spring and gravitational systems. The error
in these quantities predicted by the baseline Lnn and MCLnn are shown in
Figure 2(b) and (d). As in the case of the linear spring system, we observe
that there is a long term drift in the Lagrangian predicted by the baseline
Lnn. Similarly, the linear and angular momenta of the baseline also diverge in
a few time steps. In contrast, the Lagrangian predicted by the MCLnn exhibits
a good match with the ground truth. It should be noted that the multi-particle
interacting systems are chaotic in nature and hence the exact trajectory of
the ground truth and simulated system may vary. However, the low values of
error in the Lagrangian along with the long-term stability of both non-linear
spring and gravitational system suggests the realistic nature of the
trajectory predicted by MCLnn. In addition, we observe that the linear and
angular momenta of the system remain conserved in the MCLnn. Altogether, the
results confirm that the MCLnn can successfully learn the generative model
$p(\mathbb{T})$ that can sample the trajectory effectively from the physics-
constrained configurational and temporal space.

Figure 3: Results of the forward simulation of non-linear spring, and
gravitational system from the same initial configuration using ground truth
(continuous line) and MCLnn (dashed lines). (a) The evolution of Lagrangian,
Hamiltonian (total energy), linear and momentum for all the systems. (b) The
error in Lagrangian, Hamiltonian (total energy), linear and momentum predicted
by the baseline Lnn and MCLnn with respect to the ground truth.  Figure 4:
Results of the forward simulation of (a) linear spring with 6 particles, (b)
non-linear spring with 6 particles, and (b) gravitational system with 8
particles from the same initial configuration using ground truth (continuous
line) and MCLnn (dashed lines). The evolution of Lagrangian for all the
systems are plotted.

###  4.3 Generalizability to unseen system sizes

To demonstrate the generalizability of MCLnn, we use the generative model
$p(\mathbb{T})$ for each of the tasks learned during the training. The models
are then used to predict trajectories of systems with different number of
particles. Figure 4 shows the trajectories for linear and non-linear springs
with six particles (MCLnn trained on three particles), and gravitational
system with eight particles (MCLnn trained on four particles), predicted by
the MCLnn in comparison to the ground truth. We observe that the Lagrangian
predicted by the MCLnn exhibits excellent match with the ground truth. In
addition, all other quantities such as energy and momenta remain conserved
(see Appendix Fig. 7). This suggests that once the generative model
$p(\mathbb{T})$ has been learned, MCLnn can generalize it to predict the
trajectory of systems with any number of particles. Note that the baseline Lnn
cannot simulate systems of different sizes as the positions of all the
particles are given simultaneously as the input to the neural network.

###  4.4 Interpretability

Figure 5 shows the $V_{ij}$ predicted by the MCLnn in comparison to the pair-
wise potential energy obtained analytically from ground truth. The shaded
region represents the $q_{ij}$ values that were present in the training set.
We observe that the MCLnn is able to learn the pairwise potential energy
function accurately for the $q_{ij}$ values in the training set. Further, the
MCLnn is able to interpolate values excellently and extrapolate reasonably. It
is worth noting that the MCLnn is trained only on the trajectories of the
system. Hence, no information regarding the potential energy, force, or even
acceleration is given to the neural network during training. This suggests
that MCLnn is indeed able to learn the physics of the problem by training
purely on the trajectory.

Figure 5: Variation of $V_{ij}$ with respect to $q_{ij}$ as learned by the
MCLnn (dashed line) in comparison to the ground truth (continuous line) for
(a) linear spring, (b) non-linear spring, and (c) gravitational potential. The
shaded region represents the range of $q_{ij}$ values in the training set.
Non-shaded region between two shaded region corresponds to interpolation;
otherwise, it represents extrapolation.

##  5 Conclusion

We introduced a new framework, namely, momentum conserving Lagrangian neural
networks (MCLnn), for incorporating physics-based priors in neural networks
for accurate simulations of multi-particle systems. We demonstrated that the
MCLnn respects the symmetries in space and time, leading to conservation of
energy, and linear and angular momentum. We showed that the incorporation of
these additional conservation laws of momenta makes the Lagrangian of the
system stable by avoiding any long-term drift in it. This, in turn, results in
a realistic simulation of multi-body systems. Further, we showed that the
MCLnn once trained, can generalize to systems of any size. Finally, we
demonstrated that the MCLnn is highly interpretable and provides direct
insights into the interaction laws governing the dynamics of multi-particle
systems. This, in turn, allows one to verify the realistic nature of the
function learned by the MCLnn.

At this juncture, it is worth discussing some of the open questions and
shortcomings of MCLnn that can be addressed as part of future works. (i) MCLnn
assumes interaction between all the particles in the system. A graph-based
MCLnn with contrastive loss can potentially address this challenge by
incorporating the topology of the particle system along with differential
importance for nearby particles. (ii) MCLnn considers only pair-wise
interactions for computing the $V_{ij}$. However, there could be additional
interactions involving three ($V_{ijk}$), four ($V_{ijkl}$) or even higher
particles. This could be addressed by feature engineering or graph-based
MCLnn. (iii) Similarly, simulations involving different types of interacting
particles is challenging in MCLnn as it may lead to non-unique solutions. (iv)
In addition, present work can be extended to address more challenging
problems, for instance, to learn the Lagrangian of system with non-
conservative forces, and to learn generalized kinetic energy functions, while
maintaining the granularity and generalizability.

## References

  * [1] Herbert Goldstein.  Classical mechanics.  Pearson Education India, 2011.
  * [2] Cheol Woo Park, Mordechai Kornbluth, Jonathan Vandermause, Chris Wolverton, Boris Kozinsky, and Jonathan P Mailoa.  Accurate and scalable graph neural network force field and molecular dynamics with direct force architecture.  npj Computational Materials, 7(1):1–9, 2021.
  * [3] Manuel A Roehrl, Thomas A Runkler, Veronika Brandtstetter, Michel Tokic, and Stefan Obermayer.  Modeling system dynamics with physics-informed neural networks based on lagrangian mechanics.  IFAC-PapersOnLine, 53(2):9195–9200, 2020.
  * [4] Michael Lutter, Christian Ritter, and Jan Peters.  Deep lagrangian networks: Using physics as model prior for deep learning.  In International Conference on Learning Representations, 2019.
  * [5] Miles Cranmer, Alvaro Sanchez Gonzalez, Peter Battaglia, Rui Xu, Kyle Cranmer, David Spergel, and Shirley Ho.  Discovering symbolic models from deep learning with inductive biases.  Advances in Neural Information Processing Systems, 33, 2020.
  * [6] Miles Cranmer, Sam Greydanus, Stephan Hoyer, Peter Battaglia, David Spergel, and Shirley Ho.  Lagrangian neural networks.  In ICLR 2020 Workshop on Integration of Deep Neural Models and Differential Equations, 2020.
  * [7] Marc Finzi, Ke Alexander Wang, and Andrew G Wilson.  Simplifying hamiltonian and lagrangian neural networks via explicit constraints.  In H. Larochelle, M. Ranzato, R. Hadsell, M. F. Balcan, and H. Lin, editors, Advances in Neural Information Processing Systems, volume 33, pages 13880–13889. Curran Associates, Inc., 2020.
  * [8] Samuel Greydanus, Misko Dzamba, and Jason Yosinski.  Hamiltonian neural networks.  Advances in Neural Information Processing Systems, 32:15379–15389, 2019.
  * [9] Yaofeng Desmond Zhong, Biswadip Dey, and Amit Chakraborty.  Benchmarking energy-conserving neural networks for learning dynamics from data.  In Learning for Dynamics and Control, pages 1218–1229. PMLR, 2021\.
  * [10] Yaofeng Desmond Zhong, Biswadip Dey, and Amit Chakraborty.  Symplectic ode-net: Learning hamiltonian dynamics with control.  In International Conference on Learning Representations, 2019.
  * [11] Yaofeng Desmond Zhong and Naomi Leonard.  Unsupervised learning of lagrangian dynamics from images for prediction and control.  Advances in Neural Information Processing Systems, 33, 2020.
  * [12] Emmy Noether.  Invariant variation problems.  Transport theory and statistical physics, 1(3):186–207, 1971.
  * [13] Dennis C Rapaport.  The art of molecular dynamics simulation.  Cambridge university press, 2004.
  * [14] James Bradbury, Roy Frostig, Peter Hawkins, Matthew James Johnson, Chris Leary, Dougal Maclaurin, and Skye Wanderman-Milne.  Jax: composable transformations of python+ numpy programs, 2018.  URL http://github. com/google/jax, 4:16, 2020.
  * [15] Samuel Schoenholz and Ekin Dogus Cubuk.  Jax md: a framework for differentiable physics.  Advances in Neural Information Processing Systems, 33, 2020.
  * [16] Diederik P. Kingma and Jimmy Ba.  Adam: A method for stochastic optimization.  In ICLR (Poster), 2015.

##  Appendix A Appendix

###  A.1 Baseline Lnn

Figure 6 shows the predicted force with respect to the actual force for the
training data of baseline for the three tasks.

Figure 6: Predicted force for (a) linear spring, (b) non-linear spring, and
(c) gravitational system with respect to the measured forces on the test
dataset of baseline Lnn.

###  A.2 Generalizability to unseen system sizes

Figure 7 shows the generalizability of the MCLnn to unseen system sizes for
the three tasks considered. For all the three tasks, the number of particles
considered are twice that of the training system size. We observe that the
Lagrangian for the unseen system is predicted in excellent agreement to the
ground truth for all the three tasks. In addition, the Hamiltonian (total
energy), and momenta are conserved in all the three tasks, confirming the
realistic nature of the simulations by MCLnn.

Figure 7: Results of the forward simulation of (a) linear spring, (b) non-
linear spring, and (b) gravitational system from the same initial
configuration using ground truth (continuous line) and MCLnn (dashed lines).
The evolution of Lagrangian, Hamiltonian (total energy), linear and angular
momenta for all the systems with different number of particles than that was
present in the training set are plotted.

###  A.3 Computing acceleration from Euler-Lagrange equations

The acceleration of particles can be obtained from the Euler-Lagrange
equations as follows.

$\displaystyle\frac{\partial}{\partial\mathbf{\dot{q}}}\frac{dL}{dt}$
$\displaystyle=\frac{\partial L}{\partial\mathbf{q}}$
$\displaystyle\frac{\partial}{\partial\mathbf{\dot{q}}}\left(\frac{\partial
L}{\partial\mathbf{q}}\frac{dq}{dt}+\frac{\partial
L}{\partial\mathbf{\dot{q}}}\frac{d\mathbf{\dot{q}}}{dt}\right)$
$\displaystyle=\frac{\partial L}{\partial\mathbf{q}}$
$\displaystyle\frac{\partial}{\partial\mathbf{\dot{q}}}\left(\frac{\partial
L}{\partial\mathbf{q}}\mathbf{\dot{q}}+\frac{\partial
L}{\partial\mathbf{\dot{q}}}\mathbf{\ddot{q}}\right)$
$\displaystyle=\frac{\partial L}{\partial\mathbf{q}}$
$\displaystyle\frac{\partial}{\partial\mathbf{\dot{q}}}\frac{\partial
L}{\partial\mathbf{q}}\mathbf{\dot{q}}+\frac{\partial}{\partial\mathbf{\dot{q}}}\frac{\partial
L}{\partial\mathbf{\dot{q}}}\ddot{\mathbf{q}}$ $\displaystyle=\frac{\partial
L}{\partial\mathbf{q}}$
$\displaystyle\frac{\partial}{\partial\mathbf{\dot{q}}}\frac{\partial
L}{\partial\mathbf{\dot{q}}}\ddot{\mathbf{q}}$ $\displaystyle=\frac{\partial
L}{\partial\mathbf{q}}-\frac{\partial}{\partial\mathbf{\dot{q}}}\frac{\partial
L}{\partial\mathbf{q}}\mathbf{\dot{q}}$ $\displaystyle\ddot{\mathbf{q}}$
$\displaystyle=\left(\frac{\partial}{\partial\mathbf{\dot{q}}}\frac{\partial
L}{\partial\mathbf{\dot{q}}}\right)^{-1}\left[\frac{\partial
L}{\partial\mathbf{q}}-\frac{\partial}{\partial\mathbf{\dot{q}}}\frac{\partial
L}{\partial\mathbf{q}}\mathbf{\dot{q}}\right]$
$\displaystyle\ddot{\mathbf{q}}$
$\displaystyle=\left(\nabla_{\mathbf{\dot{q}}\mathbf{\dot{q}}}L\right)^{-1}\left[\nabla_{q}L-\left(\nabla_{\mathbf{\dot{q}}\mathbf{q}}L\right)\mathbf{\dot{q}}\right]$

###  A.4 Translational Symmetry

The conservation of linear momentum can be derived from translational symmetry
of Lagrangian as follows.

$\displaystyle L(\mathbf{q},\mathbf{\dot{q}})$
$\displaystyle=L(\mathbf{q}+\epsilon,\mathbf{\dot{q}})$
$\displaystyle\mathbf{q^{\,{}^{\prime}}}$ $\displaystyle=\mathbf{q}+\epsilon$
$\displaystyle\implies\frac{\partial\mathbf{q^{\,{}^{\prime}}}}{\partial\mathbf{q}}$
$\displaystyle=1$ $\displaystyle\sum\frac{\partial L}{\partial\mathbf{q}}$
$\displaystyle=\sum\frac{\partial
L}{\partial\mathbf{q^{\,{}^{\prime}}}}=\sum\frac{\partial
L}{\partial\mathbf{q}}\frac{\partial\mathbf{q}}{\partial\mathbf{q^{\,{}^{\prime}}}}=\sum\frac{\partial
L}{\partial\mathbf{q}}.1$ $\displaystyle\delta L$
$\displaystyle=\sum\delta\mathbf{q}\frac{\partial L}{\partial\mathbf{q}}$
$\displaystyle\delta L$ $\displaystyle=\sum\epsilon\frac{\partial
L}{\partial\mathbf{q}}$ $\displaystyle\delta L$
$\displaystyle=\epsilon\sum\frac{\partial L}{\partial\mathbf{q}}$

since,

$\displaystyle\delta L$ $\displaystyle=0$
$\displaystyle\implies\sum\frac{\partial L}{\partial\mathbf{q}}$
$\displaystyle=0$ $\displaystyle\sum\frac{\partial L}{\partial\mathbf{q}}$
$\displaystyle=0$ $\displaystyle\sum\frac{d}{dt}\frac{\partial
L}{\partial\mathbf{\dot{q}}}$ $\displaystyle=0$
$\displaystyle\frac{d}{dt}\sum\frac{\partial L}{\partial\mathbf{\dot{q}}}$
$\displaystyle=0$ $\displaystyle\sum\frac{\partial
L}{\partial\mathbf{\dot{q}}}$ $\displaystyle=constant$ $\displaystyle\sum
p_{i}$ $\displaystyle=\textrm{constant (linear momentum)}$

where $p_{i}=\frac{\partial L}{\partial\dot{q}_{i}}$

###  A.5 Rotational Symmetry

The conservation of angular momentum can be derived from rotational symmetry
of Lagrangian as follows.

Under an infinitesimal rotation $\delta\theta$

$\displaystyle\delta\vec{q}_{i}$ $\displaystyle=\delta\theta\times\vec{q}_{i}$

and

$\displaystyle\delta\vec{\dot{q}}_{i}$
$\displaystyle=\delta\theta\times\vec{\dot{q}}_{i}$

since, Lagrangian does not change with infinitesimal rotation $\delta\theta$

$\displaystyle\delta L$ $\displaystyle=\sum_{i}\frac{\partial
L}{\partial\vec{q}_{i}}\delta\vec{q}_{i}+\sum_{i}\frac{\partial
L}{\partial\dot{\vec{q}}_{i}}\delta\dot{\vec{q}}_{i}=0$

using generalized momentum

$\displaystyle\frac{\partial L}{\partial\dot{\vec{q}}_{i}}=\vec{p}_{i}$

then EL equation gives

$\displaystyle\frac{d}{dt}\vec{p}_{i}-\frac{\partial L}{\partial\vec{q}_{i}}$
$\displaystyle=0$ $\displaystyle\dot{\vec{p}}_{i}$
$\displaystyle=\frac{\partial L}{\partial\vec{q}_{i}}$
$\displaystyle\implies\delta L$
$\displaystyle=\sum_{i}\dot{\vec{p}}_{i}\cdot\delta\vec{q}_{i}+\sum_{i}\vec{p}_{i}\cdot\delta\dot{\vec{q}}_{i}=0$
$\displaystyle\sum_{i}\dot{\vec{p}}_{i}\cdotp(\delta\theta\times\vec{q}_{i})+\vec{p}_{i}\cdot(\delta\theta\times\vec{\dot{q}}_{i})$
$\displaystyle=0$
$\displaystyle\sum_{i}\delta\theta\cdotp(\vec{q}_{i}\times\dot{\vec{p}}_{i})+\delta\theta\cdotp(\dot{\vec{q}}_{i}\times\vec{p}_{i})$
$\displaystyle=0$
$\displaystyle\sum_{i}\delta\theta\cdotp[(\vec{q}_{i}\times\vec{\dot{p}}_{i})+(\vec{\dot{q}}_{i}\times\vec{p}_{i})]$
$\displaystyle=0$
$\displaystyle\sum_{i}\delta\theta\cdotp\frac{d}{dt}(\vec{q}_{i}\times\vec{p}_{i})$
$\displaystyle=0$

because $\delta\theta$ is arbitrary, then

$\displaystyle\sum_{i}\frac{d}{dt}(\vec{r}_{i}\times\mathbf{p}_{i})$
$\displaystyle=0$
$\displaystyle\frac{d}{dt}(\sum_{i}\vec{r}_{i}\times\mathbf{p}_{i})$
$\displaystyle=0$ angular momentum
$\displaystyle=\sum_{i}(\vec{r}_{i}\times\vec{p}_{i})=\text{constant}$

###  A.6 Initial Conditions and Task Configuration

The number of hidden units are chosen from hyperparamter search given in Table
1. The activation function for all hidden units are square plus
$f(x)=\frac{(x+\sqrt{x^{2}+4})}{2}$ which is similar to soft plus [6]. We use
ADAM [16] optimiser for model training with learning rate given in respective
configuration tables. We use velocity-Verlet for time integration during
trajectory evolution. Mass of all the particles are maintained as 1.0 units.

Table 1: Train and validation loss for each NN architecture for the linear spring task for MCLnn. Note that for hyperparamter search (i.e. number of hidden units), minimum train loss was set to 1.0e-8 as stopping criteria. Hidden Layers | Linear Spring
---|---
2, 2 | 0.0012 , 0.0015
4, 4 | 7.98e-09, 2.36e-08
8, 8 | 9.86e-09, 3.80e-08
16, 16 | 9.76e-09, 3.02e-08

####  A.6.1 Linear Spring System

$\text{Initial
position}=\begin{bmatrix}0.486657678894505&0.755041888583519&0.0\\\
-0.681737994414464&0.293660233197210&0.0\\\
-0.022596327468640&-0.612645601255358&0.0\\\ \end{bmatrix}$
$\text{Initial
velocity}=\begin{bmatrix}-0.182709864466916&0.363013287999004&0.0\\\
-0.579074922540872&-0.748157481446087&0.0\\\
0.761784787007641&0.385144193447218&0.0\\\ \end{bmatrix}$
$\text{dt}=0.01$
$\text{mass}=1.0$
$\text{stride}=10$
$\text{runs}=20=\text{points per trajectory}$
$\text{lr}=1.0e-3$
$\text{layers}=[10,10]$
$\text{epochs}=100000$
$\text{samples}=100=\text{number of trajectories}$
$\text{seed}=100$
$\text{time step}=dt\times stride=0.1$

####  A.6.2 Non-linear Spring System

$\text{Initial
position}=\begin{bmatrix}0.486657678894505&0.755041888583519&0.0\\\
-0.681737994414464&0.293660233197210&0.0\\\
-0.022596327468640&-0.612645601255358&0.0\\\ \end{bmatrix}$
$\text{Initial
velocity}=\begin{bmatrix}-0.182709864466916&0.363013287999004&0.0\\\
-0.579074922540872&-0.748157481446087&0.0\\\
0.761784787007641&0.385144193447218&0.0\\\ \end{bmatrix}$
$\text{dt}=0.01$
$\text{mass}=1.0$
$\text{stride}=10$
$\text{runs}=20=\text{points per trajectory}$
$\text{lr}=1.0e-3$
$\text{layers}=[10,10]$
$\text{epochs}=100000$
$\text{samples}=100=\text{number of trajectories}$
$\text{seed}=100$
$\text{time step}=dt\times stride=0.1$

####  A.6.3 Gravitational System

$\text{Initial position}=\begin{bmatrix}1.0&0.0&0.0\\\ 9.0&0.0&0.0\\\
11.0&0.0&0.0\\\ -1.0&0.0&0.0\\\ \end{bmatrix}$
$\text{Initial velocity}=\begin{bmatrix}0.0&0.05&0.0\\\ 0.0&-0.05&0.0\\\
0.0&0.65&0.0\\\ 0.0&-0.65&0.0\\\ \end{bmatrix}$
$\text{dt}=0.01$
$\text{mass}=1.0$
$\text{stride}=10$
$\text{runs}=20=\text{points per trajectory}$
$\text{lr}=1.0e-3$
$\text{layers}=[10,10]$
$\text{epochs}=100000$
$\text{samples}=100=\text{number of trajectories}$
$\text{seed}=100$
$\text{time step}=dt\times stride=0.1$

###  A.7 Default Notation

$\displaystyle a$ |  A scalar (integer or real)
---|---
$\displaystyle\vec{a}$ |  A vector
$\displaystyle{\bm{\mathsfit{A}}}$ |  A matrix
$\displaystyle{\bm{I}}_{n}$ |  Identity matrix with $n$ rows and $n$ columns
$\displaystyle{\bm{I}}$ |  Identity matrix with dimensionality implied by context
$\displaystyle{\bm{e}}^{(i)}$ |  Standard basis vector $[0,\dots,0,1,0,\dots,0]$ with a 1 at position $i$
$\displaystyle\vec{a}_{i}$ |  Vector corresponding to particle $i$
$\displaystyle{A}_{i,j}$ |  Element $i,j$ of matrix ${\bm{A}}$
$\displaystyle{\bm{A}}_{i,:}$ |  Row $i$ of matrix ${\bm{A}}$
$\displaystyle{\bm{A}}_{:,i}$ |  Column $i$ of matrix ${\bm{A}}$
$\displaystyle{\mathsfit{A}}_{i,j,k}$ |  Element $(i,j,k)$ of a 3-D tensor ${\bm{\mathsfit{A}}}$
$\displaystyle{\bm{\mathsfit{A}}}_{:,:,i}$ |  2-D slice of a 3-D tensor
$\displaystyle\frac{dy}{dx}$ |  Derivative of $y$ with respect to $x$
$\displaystyle\frac{\partial y}{\partial x}$ |  Partial derivative of $y$ with respect to $x$
$\displaystyle\nabla_{\vec{x}}L$ |  Gradient of $L$ with respect to $\vec{x}$
$\displaystyle\nabla_{\bm{X}}y$ |  Matrix derivatives of $y$ with respect to ${\bm{X}}$
$\displaystyle\frac{\partial f}{\partial{\bm{x}}}$ |  Jacobian matrix ${\bm{J}}\in\mathbb{R}^{m\times n}$ of $f:\mathbb{R}^{n}\rightarrow\mathbb{R}^{m}$
$\displaystyle\int f({\bm{x}})d{\bm{x}}$ |  Definite integral over the entire domain of ${\bm{x}}$
$\displaystyle\int_{\mathbb{S}}f({\bm{x}})d{\bm{x}}$ |  Definite integral with respect to ${\bm{x}}$ over the set ${\mathbb{S}}$
$\displaystyle P({\textnormal{a}})$ |  A probability distribution over a discrete variable
$\displaystyle p({\textnormal{a}})$ |  A probability distribution over a continuous variable, or over a variable whose type has not been specified