Current address: ]Elettra-Sincrotrone Trieste S.C.p.A., 34149 Basovizza,
Trieste, Italy

# Magnetization switching in the inertial regime

Kumar Neeraj  Department of Physics, Stockholm University, 106 91 Stockholm,
Sweden  Matteo Pancaldi  [  Department of Physics, Stockholm University, 106
91 Stockholm, Sweden  Valentino Scalera  Faculty of Computer Science, Free
University of Bozen-Bolzano, Bolzano, Italy  Salvatore Perna  DIETI,
University of Naples Federico II, Naples, Italy  Massimiliano d’Aquino  DIETI,
University of Naples Federico II, Naples, Italy  Claudio Serpico  DIETI,
University of Naples Federico II, Naples, Italy  Stefano Bonetti
<EMAIL_ADDRESS>Department of Physics, Stockholm University, 106
91 Stockholm, Sweden  Department of Molecular Sciences and Nanosystems, Ca’
Foscari University of Venice, 30172 Venezia-Mestre, Italy

###### Abstract

We have numerically solved the Landau-Lifshitz-Gilbert (LLG) equation in its
standard and inertial forms to study the magnetization switching dynamics in a
$3d$ thin film ferromagnet. The dynamics is triggered by ultrashort magnetic
field pulses of varying width and amplitude in the picosecond and Tesla range.
We have compared the solutions of the two equations in terms of switching
characteristic, speed and energy analysis. Both equations return qualitatively
similar switching dynamics, characterized by regions of slower _precessional_
behavior and faster _ballistic_ motion. In case of inertial dynamics,
ballistic switching is found in a 25% wider region in the parameter space
given by the magnetic field amplitude and width. The energy analysis of the
dynamics is qualitatively different for the standard and inertial LLG
equations. In the latter case, an extra energy channel, interpreted as the
kinetic energy of the system, is available. Such extra channel is responsible
for a resonant energy absorption at THz frequencies, consistent with the
occurence of spin nutation.

††preprint: APS/123-QED

##  I Introduction

The traditional method of writing information in magnetic hard disk drives
consists in reversing the magnetization direction via the application of
magnetic fields produced by external currents and localized via a so-called
”write-head”. In order to achieve efficient switching, the magnetic field is
applied nearly anti-parallel to the direction of the initial magnetization
state, and the switching process thus obtain is often referred to as “damped”
switching Mallinson (2000); Bertotti  _et al._ (2003a); Stöhr and Siegmann
(2006). The switching time in this process is limited by the macroscopic
relaxation time of the magnetization of the order of 100 ps, and correctly
described by the standard Landau Lifhitz Gilbert (LLG) equation Gilbert
(2004). Geritts et. al. Gerrits  _et al._ (2002) demonstrated instead a
technique by which ultrafast magnetization reversal can be achieved using
picosecond-long magnetic field pulses _transverse_ to the magnetization
direction. Such a switching technique was also reported in other studies Kaka
and Russek (2002); Hiebert  _et al._ (2002). The switching times reported in
these works depend on the amplitude of the magnetic field and on the duration
of the pulse, with a general trend that a pulse of larger amplitude will
reduce the switching time. Coincidentally, at the same time, a series of
experiments performed at the Stanford linear accelerator demonstrated the
ultrafast switching of magnetization by intense magnetic fields created by
relativistic electron bunches. In the paper published by Tudosa et al. Tudosa
_et al._ (2004), it was argued that deterministic magnetization switching
cannot occur faster than 2 ps, setting this as the ultimate speed for magnetic
reversal. However, due to the complexity of such an accelerator based
experiment, direct observation of such an ultrafast switching in the time
domain was not possible.

A few years later, it was shown that a novel idea, i.e. using magnetic
inertia, can greatly enhance the switching speed in antiferromagnets Kimel
_et al._ (2009), up to 10 times faster than that reported in the above-
mentioned studies. So far, nutation-type magnetization motions in
ferromagnetic systems have been studied mostly theoretically, in the framework
of classical LLG dynamics at GHz frequencies Serpico  _et al._ (2004), as well
as in other numerical studies considering the inertial version of the LLG,
often refered to as iLLG equation Böttcher  _et al._ (2011); Bhattacharjee
_et al._ (2012a); Bastardis  _et al._ (2018); Makhfudz  _et al._ (2020). These
works predicted the appearance of a spin nutation with an intrinsic resonance
in the $10^{11}-10^{15}$ Hz range. The direct detection of spin nutation in
ferromagnets recently achieved experimentally Neeraj  _et al._ (2021) has
allowed to narrow down the rather broad frequency range to the $10^{12}$ Hz
one, i.e. in the THz region, for typical $3d$ ferromagnetic alloys such as
NiFe and CoFeB.

In this work, we explore the role of inertia in the magnetization switching of
a thin film ferromagnet triggered by magnetic field pulses in the picosecond
range. Using macrospin simulations, we create a map of the magnetization
dynamics as a function of the pulse duration and amplitude. The dynamics is
obtained and analyzed solving both the standard LLG equation and its inertial
form, the so-called iLLG equation. We used realistic material parameters for
an archetypal ferromagnet, namely Ni81Fe19 (permalloy), but they are valid for
all thin film ferromagnets with similar characteristics. We show that the iLLG
simulations predict a larger stability region than the LLG one for one the two
characteristic switching processes, and that the effect of the nutation
resonance can be detected in switching experiments. We discuss our results
analyzing the dynamics of the classical energy terms.

##  II Methods

Figure 1: (a) Spherical coordinates used for the numerical simulations
described in the main text. (b) Geometry of the thin film system considered.
The easy magnetization axis lies along the y-direction, and the in-plane hard
magnetization axis is along the x-direction. (c) Main plot: magnetization
state diagram for different magnetic pulse amplitude and FWHM. Side plots:
magnetization precession trajectories in selected points of the diagram. The
color bar shows the y component of the magnetization vector at the end of the
simulation. The magnetization starts always from the positive y direction,
i.e. aligned parallel to the easy magnetization axis.

Until recently, the dynamics of magnetization switching in ferromagnetic
systems has been mostly described by the conventional LLG equation Gilbert
(2004), which can be written as

$\frac{d{\textbf{M}}}{dt}=-|\gamma|{\textbf{M}}\times\left({\textbf{H}_{\textrm{eff}}}-\frac{\alpha}{|\gamma|M_{s}}\frac{d\textbf{M}}{dt}\right),$
(1)

where $|\gamma|/(\mu_{0}2\pi)\approx 28$ GHz/T is the gyromagnetic ratio,
$\textbf{H}_{\textrm{eff}}$ is the effective magnetic field, calculated as the
variational derivative of energy with respect to the magnetization, $M_{s}$ is
the saturation magnetization, and $\alpha$ is the Gilbert damping. The first
term on the right-hand side describes the precession motion while the second
describes the damping of precession Gilbert (2004). However, according to
recent theoretical studies the inclusion of an extra inertial term has been
suggested Ciornei  _et al._ (2011); Olive  _et al._ (2012); Fähnle  _et al._
(2011); Mondal  _et al._ (2017) as

$\frac{d{\textbf{M}}}{dt}=-|\gamma|{\textbf{M}}\times\left[{\textbf{H}_{\textrm{eff}}}-\frac{\alpha}{|\gamma|M_{s}}\left(\frac{d{\textbf{M}}}{dt}+\tau\frac{d^{2}{\textbf{M}}}{dt^{2}}\right)\right],$
(2)

where $\tau$ is the angular momentum relaxation time. The extra term is the
second derivative of the magnetization vector which leads to spin nutation, in
addition to precession and damping. The microscopic origin of inertia is still
not clearly understood and several theoretical studies have been presented
Mondal  _et al._ (2017); Makhfudz  _et al._ (2020). Mondal et al Mondal  _et
al._ (2017) have suggested that the inclusion of higher order spin-orbit
coupling terms in the Hamiltonian of the system will produce an equation of
this kind.

In order to understand the dynamical behavior of magnetization shown by Eqs.
(1) and (2) we perform simulations on an infinite thin magnetic film, assuming
that the magnetization is spatially homogeneous (macrospin approximation). The
effective magnetic field $\textbf{H}_{\textrm{eff}}$ appearing in Eq. (2)
takes into account the different interactions occurring among elementary
magnetic moments, namely

$\textbf{H}_{\textrm{eff}}=\textbf{H}_{\textrm{a}}+\textbf{H}_{\textrm{ani}}+\textbf{H}_{\textrm{m}},$
(3)

where $\textbf{H}_{\textrm{a}}$ is the applied magnetic field,
$\textbf{H}_{\textrm{ani}}=2K_{1}/(\mu_{0}M_{s})\bm{e}_{y}$ is the uniaxial
magneto-crystalline anisotropy field with $K_{1}$ being the anisotropy
constant and $\textbf{H}_{\textrm{m}}=-{N}\cdot\bm{M}$ the magnetostatic
(demagnetizing) field. $N$ is the demagnetization tensor, with
${N=\text{diag}(N_{x},N_{y},N_{z})}$ when referred to the principal axes of
the system. For the case of a thin film infinite in the $x,y$ directions, the
demagnetizing field takes the form $\textbf{H}_{\textrm{m}}=-M_{z}\bm{e}_{z}$
Devolder  _et al._ (2006); d’Aquino (2004); Bauer  _et al._ (2000). We
consider a thin film with an in-plane easy and hard axis, with an anisotropy
field of $\mu_{0}\textbf{H}_{\textrm{ani}}=0.1$ T, as shown in Fig. 1(b). This
anisotropy field defines two equilibrium states of magnetization in the system
along the y-axis.

We solve the LLG and iLLG differential equations using two independent
numerical codes. The first one is based on a classical Runge-Kutta fourth-
order method in spherical coordinates defined by a set of three unit vectors
{$\textbf{e}_{\textrm{x}},\textbf{e}_{\textrm{y}},\textbf{e}_{\textrm{y}}$} as
shown in Fig. 1(a). Details are given in Appendix A. The second code relies on
the appropriate extension to the iLLG dynamics of the implicit midpoint rule
time-stepping d’Aquino  _et al._ (2005) in Cartesian coordinates. The
considered magnetic thin film has material parameters similar to those of
polycrystalline permalloy, i.e. $\alpha=0.023$, $\mu_{0}M_{s}=0.92$ T,
$\tau=11.3$ ps Neeraj  _et al._ (2021). The applied magnetic field has a
Gaussian-like shape with a Full-Width Half-Maximum (FWHM) varying from 0 to 3
ps in steps of 10 fs, while the amplitude ($\mu_{0}\textbf{H}_{\textrm{a}}$)
is varied in steps of 100 mT from 0 to 10 T. The field is always applied along
the $x$-axis, i.e. perpendicular to the easy magnetization axis and
equilibrium states.

##  III Results and discussion

For a thin film system, magnetization switching can be understood as a three
step process Stöhr and Siegmann (2006). (_i_) An applied pulse
H${}_{\textrm{a}}$ lying in the film plane and perpendicular to M will exert a
Zeeman torque, and will lead M to precess out of the in-plane easy axis.
(_ii_) This torque will simultaneously create a demagnetizing field
perpendicular to the film plane. The demagnetizing field will further lead to
precession of M around the axis perpendicular to the film plane. (_iii_)
Eventually, once the applied magnetic field is turned off, the magnetization
relaxes along the direction of the effective magnetic field defined in Eq.
(3).

Figure 1(c) shows the diagram of the magnetization switching obtained by
numerically solving the iLLG equation. Similar diagrams have been extensively
used for studying magnetization switching in different case studies Bertotti
_et al._ (2003b); Bhattacharjee  _et al._ (2012b); Wienholdt  _et al._ (2012);
Shutyi and Sementsov (2020), but not yet for the case of inertial dynamics in
ferromagnets. The diagram shows the state of magnetization with respect to the
applied Gaussian pulse width and amplitude, with the yellow and blue regions
representing the non-switched and, respectively, the switched state of
magnetization. It should be noted that the presence of an in-plane hard axis
creates an energy landscape characterized by two stable states along the
$y$-axis, two saddle-type equilibria along the $x$-axis, and two unstable
equilibria along the $z$-axis.

In all the regions of the diagram, we can extract the full magnetization
dynamics, as shown in the trajectories on the unitary sphere in Fig. 1(c) for
a few characteristics cases. Case 1 is representative of the region where the
magnetization state after excitation relaxes back to its initial equilibrium
state. In case 2, the applied pulse deposits enough energy into the system for
the magnetization to cross the barrier created by the presence of the in-plane
hard axis, and to switch to the opposite magnetization state. We denote this
type of switching as the precessional switching. For cases 3 and 4, the
magnetization reaches just close to the energy barrier (i.e., the top of the
sphere) of the out-of-plane hard axis where it precesses several times before
relaxing to the initial and, respectively opposite state. It is worthwhile to
point out that this type of switching exhibits quasi-random relaxation
behavior as a result of the multi-stability and the low dissipation of
magnetization dynamicsSerpico  _et al._ (2009). The most peculiar region is
the one represented by case 5, where the magnetization is able to cross the
energy barrier of the out-of-plane hard axis and switches to the opposite
magnetization state with negligible precession. We denote this type of
switching as ballistic switching. Note that in literature the terms
_precessional_ and _ballistic_ are sometimes used interchangeably d’Aquino
_et al._ (2004); Nozaki and Matsuyama (2006); Bazaliy (2011). A reason for
this situation is the fact that previous studies considered only slower time
scales, where switching such as the one in case 5 was never observed. Here,
considering much faster dynamics, we are able to identify the qualitative
difference between the two switching processes.

Figure 2:  Magnetization switching times calculated for different pulse
amplitude and FWHM width using (a) the LLG and (b) the iLLG equations. All
axes and amplitudes are in logarithmic scale. (c) Switching time along the
diagonal line cuts (i.e. perpendicular to the lines of constant pulse energy)
shown by the white solid lines in (a) and (b). The dashed vertical lines
indicate the boundaries of the ballistic switching region for the two
simulations.

From the dynamical response of the magnetization to different Gaussian pulses,
the magnetization switching time can be obtained. The exact definition of
switching time reported in literature is ambiguous and no convention is
reported Stöhr and Siegmann (2006). In some studies the switching time is
defined to be the time taken by the magnetization to completely reach
equilibrium along the opposite easy axis direction. In this study, we define
the magnetization switching time as the time taken for magnetization to cross
the energy barrier for the last time before relaxing on either side of the
easy axis (+y or -y axis in our geometry). Based on this definition (see
Appendix B for the details) we calculate a diagram of the magnetization
switching time for different FWHM and amplitude of the externally applied
magnetic field pulse as shown in Fig. 2. Fig. 2(a) and (b) were obtained by
solving the LLG and, respectively, the iLLG equation.

In both plots, the initial black regions (marked as 1) correspond to the
regions where the magnetization never crosses the energy barrier, and
therefore there is no switching. On the other hand, region 2 indicates the
area of relatively longer switching times, where the magnetization precesses
several times around the unitary magnetization sphere before relaxing to
either one of the two equilibrium states. Once again, as in Fig. 1, the most
peculiar observation is found in region 3, which corresponds to the case of
ballistic switching. Here, we notice that the switching time is the fastest
and is largely independent of the magnetic pulse width and amplitude. We
highlight the choice of a logarithmic scale to properly highlight the
different orders of magnitude between the switching times in regions 2 and 3.
A closer look at the difference in the two simulation results shows that the
width of region 3 for the case of the iLLG dynamics is larger than that of the
LLG one, as highlighted by the dashed vertical lines in Fig. 2(c). The figure
is obtained taking a diagonal cut of Fig. 2(a),(b). The width of the ballistic
switching region is approximately 3.34 T ps for the LLG dynamics, and 4.16 T
ps for the iLLG one, i.e. approximately 25% larger in the latter case. We also
notice that, for the LLG dynamics, the ballistic region starts at lower
$\mu_{0}$H${}_{\textrm{a}}\times$FWHM values than for the iLLG case.

We can interpret these observations with analogies taken from classical
mechanics. For a given driving force, a system at rest with larger inertia
will react less to that force, but then inertia helps in preserving the motion
of the system once the dynamics has started.

Figure 3: Temporal evolution of the kinetic, potential, dissipated and
deposited energy terms for a few selected simulation parameters. (a) LLG and
(b) iLLG dynamics of the different energy terms for an applied magnetic field
of 2 T amplitude and 1 ps FWHM, i.e. in the _precessional_ switching region.
(c) LLG and (d) iLLG energy dynamics for 8 T, 2.1 ps magnetic field pulse,
i.e. in the _ballistic_ switching region.

In addition to the speed of the switching process, it is equally important to
understand the energy associated with it. In order to do this, we calculated
the work per unit volume $\Delta W$ performed by the field pulse on the
magnetic system as

$\Delta
W=\int_{0}^{t_{p}}\mu_{0}\textbf{H}_{\textrm{a}}\cdot\frac{d\textbf{M}}{dt}dt=\frac{\mu_{0}}{\gamma
M_{s}}\int_{0}^{\infty}\alpha\bigg{|}\frac{d\bm{M}}{dt}\bigg{|}^{2}\,dt\,,$
(4)

where $t_{p}$ is the time instant at which the external pulse amplitude goes
back to zero after reaching its maximum. The latter formula comes from the
energy balance equation (see the derivation in Appendix C), i.e.

$\int_{0}^{t}\mu_{0}\textbf{H}_{\textrm{a}}\cdot\frac{d\textbf{M}}{dt}dt=\int_{0}^{t}\bigg{(}\frac{dA}{dt}+\frac{dK}{dt}+\frac{\alpha\,\mu_{0}}{\gamma
M_{s}}\bigg{|}\frac{d\bm{M}}{dt}\bigg{|}^{2}\bigg{)}dt\,.$ (5)

The left-hand side of Eq. (5) is the energy deposited by the pulse as function
of time $t$ after turning the field on. This is equal to the energy density
absorbed by the system, which in turn is the sum of the potential energy
$A=K_{1}(1-(M_{y}/M_{s})^{2})+\dfrac{1}{2}\mu_{0}M_{z}^{2}$, kinetic energy
$K=\dfrac{\mu_{0}}{2\gamma
M_{s}}\alpha\tau\displaystyle\left|\frac{d\textbf{M}}{dt}\right|^{2}$, and
dissipated energy $\displaystyle\int_{0}^{t}\dfrac{\mu_{0}}{\gamma
M_{s}}\alpha\left|\dfrac{d\textbf{M}}{dt}\right|^{2}\,dt$. Eq. (5) is derived
for the iLLG case but it is directly applicable to LLG equation without the
kinetic term, since $\tau=0$ in this case.

Using Eq. (5), we plot in Fig. 3 the temporal evolution of the energy transfer
from the magnetic field pulse into the system without (LLG) and with inertia
(iLLG), and for the case of precessional (magnetic field with amplitude 2 T,
and 1 ps FWHM duration) and ballistic switching (8 T, 2.1 ps FWHM). In all
cases, at long enough times, the dissipated energy converges to the value of
the deposited energy, indicated with solid and, respectively, dashed black
lines. However, the way such values is reached is remarkably different for the
four cases.

In the case of precessional switching in the standard LLG dynamics, Fig. 3(a),
the energy is dissipated rather slowly, consistent with the comparatively long
switching times shown in Fig. 2. In the energy picture, the potential energy
of the precession is dissipated away by the Gilbert damping in the hundreds of
ps range. For the precessional switching in the case of iLLG dynamics, Fig.
3(b,) the dissipation of the deposited energy is again a slow process, but
there are two key differences with the LLG case. First, the total deposited
energy is approximately twice as large, owing to the fact that an additional
energy channel, i.e. the kinetic one, is now available, which also creates a
fast overshoot at the beginning of the dynamics. Second, in addition to the
slow dissipation of energy, there is now a relatively fast periodic exchange
between potential and kinetic energies. This takes place at THz rates and it
is the signature of spin nutation in the energy dynamics. Interestingly, the
kinetic energy relaxes faster than the potential energy.

Figure 4: Total energy deposited by a 8 T magnetic field pulse as a function
of pulse FWHM. The plot is computed for different values of the angular
momentum relaxation time constant $\tau$ used in the iLLG dynamics. Inset:
amplitude of peak energy absorbed as a function of $\tau$.

We now turn to the case of the ballistic switching, i.e. Figs. 3(c),(d). Here,
the dissipated energy for the LLG and iLLG dynamics is almost identical, but
in the case of the iLLG dynamics in Fig. 3(d) there is, similarly to Fig.
3(b), a larger overshoot than for the LLG case. This is again due to the
additional contribution of the kinetic energy term. This overshoot is
intriguing and it is a rather general feature of the iLLG dynamics. We
demonstrate this by calculating the LLG and iLLG diagrams for all magnetic
field amplitudes and widths considered in Figs. 1 and 2, and we plot this in
Fig. 6 in Appendix C. In order to investigate this observation further, we
take a horizontal line cut in Fig. 6(b) in correspondence of the 8 T value,
and we plot it in Fig. 4. A peak is present in the deposited energy value for
a pulse FWHM of around 0.5 ps and $\tau=10$ ps. By varying the value of
$\tau$, we observe that the maximum work is performed by the applied field
when the pulse duration 6$\sigma$, with $\sigma\approx\text{FWHM}/2.4$,
approximately matches the period of the nutation resonance $2\pi\alpha\tau$.
Hence, we can attribute the occurrence of the peak to the resonant absorption
of the nutation resonance. In the inset of Fig. 4 we also plot the peak energy
value for the different $\tau$, observing a monotonic increase. This is
generally consistent with a torque-driven dynamics with a constant
gyromagnetic ratio $\gamma$ where, for a given magnetic field amplitude, the
larger torque on the magnetization, hence its maximum displacement, is
achieved when the magnetic field varies more slowly.

##  IV Summary and conclusion

From our LLG and iLLG simulations we explored the different switching dynamics
triggered by ultrafast magnetic field pulses of different FWHM and amplitude
in the ps and Tesla ranges. Depending on those parameters, the magnetization
exhibits precessional or ballistic switching.The ballistic switching is always
much faster than the precessional switching, for both the LLG and iLLG
equations. We also observed that the width of the ballistic switching region
for the iLLG case is 25% larger than the LLG one. From an application
perspective, this difference in width of the ballistic region may turn out to
be useful for the reliability of ultrafast switching, and guide the choice of
materials for magnetic storage devices with larger inertia.

We further showed how the energy is deposited into the system by the external
magnetic field pulse. For the iLLG dynamics, the external pulse drives the
spin nutation with a characteristic resonant feature, and the deposited energy
has an additional kinetic energy channel not available in the standard LLG
dynamics, where instead only the potential energy channel exists.

Our results will be useful for the design of devices relying on ultrafast
magnetization switching using picosecond magnetic field pulses. The
realization of such short and intense magnetic fields have been suggested
recently designing magnetic-field enhancing metamaterials Polley  _et al._
(2018), by the use of vector laser beams Sederberg  _et al._ (2020) or thanks
to ultrafast electronic switches Jhuria  _et al._ (2020). Finally, we
anticipate that our findings will be relevant to much of the recent works in
the magnetism community aimed at the understanding of inertial spin dynamics
Thonig  _et al._ (2017); Giordano and Déjardin (2020); Mondal  _et al._
(2021); Mondal (2021); Titov  _et al._ (2021a); Thibaudeau and Nicolis (2021);
Lomonosov  _et al._ (2021); Ruggeri (2021); Cherkasskii  _et al._ (2021);
Rahman and Bandyopadhyay (2021); Titov  _et al._ (2021b); Anders  _et al._
(2020); Lou  _et al._ (2021); Gupta  _et al._ (2021).

###### Acknowledgements.

K.N., M.P. and S.B. acknowledge support from the European Research Council,
Starting Grant 715452 MAGNETIC-SPEED-LIMIT.

##  Appendix A Magnetization dynamics in spherical coordinates

The dynamical Eqs. (1) and (2) can be decomposed in the reference frame of a
spherical coordinate system defined by a set of three unit vectors
{$\textbf{e}_{\textrm{r}},\textbf{e}_{\mathrm{\theta}},\textbf{e}_{\mathrm{\phi}}$}.
The magnetization vector in this reference frame is defined by the set of
coordinates $M=(M_{s},\theta,\phi)$, where the unit vector
$\textbf{e}_{\textrm{r}}$ is aligned with M while eθ and
$\textbf{e}_{\mathrm{\phi}}$ point in the direction of increasing $\theta$ and
$\phi$ respectively. They are both tangent to the unit sphere, as shown in Fig
1(a). The expression for the net magnetic field in angular components can be
written in terms of the Cartesian coordinates as

$\displaystyle\textbf{H}_{\textrm{r}}$
$\displaystyle=\textbf{H}_{\textrm{x}}\cos\phi\sin\theta+\textbf{H}_{\textrm{y}}+\sin\phi\cos\theta+\textbf{H}_{\textrm{z}}\cos\theta$
(6) $\displaystyle\textbf{H}_{\mathrm{\theta}}$
$\displaystyle=\textbf{H}_{\textrm{x}}\cos\phi\sin\theta+\textbf{H}_{\textrm{y}}+\sin\phi\cos\theta-\textbf{H}_{\textrm{z}}\sin\theta$
(7) $\displaystyle\textbf{H}_{\mathrm{\phi}}$
$\displaystyle=-\textbf{H}_{\textrm{x}}\sin\phi+\textbf{H}_{\textrm{y}}\cos\phi\,.$
(8)

###  A.1 LLG equation in spherical coordinates

Since $\displaystyle{\frac{d\mathbf{M}}{dt}}$ is perpendicular to M, a vector
multiplication of Eq. (1) by M will give

$(-\gamma\eta
M_{\textrm{s}}^{2}+\textbf{M}\times)\frac{d\mathbf{M}}{dt}=-M_{s}^{2}\textbf{H}_{\mathrm{eff\perp}},$
(9)

where $\textbf{H}_{\mathrm{eff\perp}}$ is the component of effect field
perpendicular to M. The time derivative of the magnetization in spherical
coordinates is written as,
$\displaystyle\frac{d\mathbf{M}}{dt}=M_{s}\frac{d\theta}{dt}\textbf{e}_{\mathrm{\theta}}+M_{s}\frac{d\phi}{dt}\sin\theta\textbf{e}_{\mathrm{\phi}}$.
Plugging this into Eq. (9) and reordering the terms will give the time
derivative of angular components in the following canonical form

$\displaystyle\frac{d\theta}{dt}$ $\displaystyle=$
$\displaystyle\kappa(-H_{\phi}+\alpha H_{\theta}),$ (10)
$\displaystyle\frac{d\phi}{dt}$ $\displaystyle=$
$\displaystyle\kappa\dfrac{2}{\sin\theta}(-H_{\theta}+\alpha H_{\phi}),$ (11)

where $\kappa=\dfrac{\gamma}{1+\alpha^{2}}$.

###  A.2 Inertial LLG equation in spherical coordinates

Similarly, the iLLG equation can be written in spherical coordinates. As in
the last section, vector multiplication of Eq. (2) by M will give

$\gamma\eta\tau
M_{\textrm{s}}^{2}\frac{d^{2}\mathbf{M}_{\perp}}{dt^{2}}=\gamma
M_{s}^{2}\textbf{H}_{\mathrm{eff\perp}}-\gamma\eta
M_{\textrm{s}}^{2}\frac{\mathbf{M}}{dt}+\textbf{M}\times\frac{\mathbf{M}}{dt}$
(12)

where the subscript $\perp$ indicates the perpendicular component of M. Now
the second derivative of the magnetization in spherical coordinates can be
written as

$\frac{d^{2}\mathbf{M}}{dt^{2}}=-M_{s}\bigg{(}\frac{d\phi^{2}}{dt}\sin^{2}\theta+\frac{d\theta^{2}}{dt}\bigg{)}\textbf{e}_{\textrm{r}}+M_{s}\bigg{(}\frac{d^{2}\phi}{dt^{2}}-\frac{d\phi^{2}}{dt}\sin\theta\cos\theta\bigg{)}\textbf{e}_{\mathrm{\theta}}+M_{s}\bigg{(}2\frac{d\theta}{dt}\frac{d\phi}{dt}\cos\theta+\frac{d^{2}\phi}{dt^{2}}\sin\theta\bigg{)}\textbf{e}_{\mathrm{\phi}}.$
(13)

Substituting this relation in Eq. (12) gives

$\displaystyle\frac{d^{2}\theta}{dt^{2}}$
$\displaystyle=\frac{1}{\alpha\tau}\bigg{(}\alpha\tau\frac{d\phi^{2}}{dt}+\gamma
H_{\theta}-\alpha\frac{d\theta}{dt}-\frac{d\phi}{dt}\sin\theta\bigg{)}$ (14)
$\displaystyle\frac{d^{2}\phi}{dt^{2}}$
$\displaystyle=\frac{1}{\alpha\tau\sin\theta}\bigg{(}-2\alpha\tau\frac{d\theta}{dt}\frac{d\phi}{dt}+\gamma
H_{\phi}-\alpha\frac{d\phi}{dt}\sin\theta$
$\displaystyle-\frac{d\theta}{dt}\sin\theta\bigg{)}$ (15)

The equilibrium condition in our simulation is defined as
{$\theta=\pi/2$,$\phi=\pi/2$, $d\theta/dt=0$, $d\phi/dt=0$}.

##  Appendix B Magnetization dynamics in time domain

Figure 5: Temporal evolution of the $y$-component of the magnetization for the
LLG and iLLG dynamics for a case of precessional switching. The red circle
indicates the point used to calculate the switching time, the black circle the
time where the magnetization state (switched or not switched) was determined.
The inset is the magnetic field pulse used to trigger the dynamics in the main
figure.

Fig. 5 shows the temporal evolution of the y component of the magnetization
when triggered by the magnetic field pulses represented in the inset. The
black circular marker shows the point in time when the magnetization state is
recorded to plot the diagram shown in Fig. 1(c). The switching time in this
study is defined as the time when the magnetization crosses the equator for
the last time before reaching the steady state, red circle in Fig. 5. In the
following, the derivation of Eq. (4) in the main text is reported.

##  Appendix C  Energy diagram

In Fig. 6, we plot the energy deposited by the pulse for the LLG and iLLG
equations as a function of applied pulse width and amplitude. The energy is
calculated using Eq. (4), by taking a direct integral over a time duration of
the full pulse width 6$\sigma$, with $\sigma\approx\text{FWHM}/2.4$.

Figure 6:  Diagram of the energy density deposited by a magnetic field pulse
as a function of its amplitude and width, calculated using (a) the LLG
equation and (b) the iLLG equation. The horizontal white dashed line in
correspondence of the 8 T value is the line cut plotted in Fig. 4 in the main
text.

###  C.1 Energy balance in iLLG dynamics

For the sake of simplicity of notation, we consider the iLLG dynamics in its
dimensionless form:

$\frac{d\bm{m}}{dt}=-\bm{m}\times\left(\bm{h}_{\mathrm{eff}}-\alpha\frac{d\bm{m}}{dt}-\xi\frac{d^{2}\bm{m}}{dt^{2}}\right)\,\,,$
(16)

where $\bm{m}(t)$ is the magnetization unit-vector (normalized by the
saturation magnetization $M_{s}$), time is measured in units of $(\gamma
M_{s})^{-1}$, $\gamma=2.21\times 10^{5}$ A-1s-1m is the absolute value of the
gyromagnetic ratio, $\alpha$ is the dimensionless Gilbert damping parameter,
and $\xi$ measures the strength of the inertial effects in the magnetization
dynamics. The effective field is

$\bm{h}_{\mathrm{eff}}=-\frac{\partial
g}{\partial\bm{m}}=-D\cdot\bm{m}+\bm{h}_{a}\,\,,$ (17)

which is expressed as the gradient of the free energy
$g(\bm{m},\bm{h}_{a})=\frac{1}{2}\bm{m}\cdot
D\cdot\bm{m}-\bm{h}_{a}\cdot\bm{m}$, with $D=\mathrm{diag}(D_{x},D_{y},D_{z})$
being the effective demagnetizing tensor (referred to principal axes and also
including anisotropy) and $\bm{h}_{a}$ the external applied field. $\xi$ is
related to the time scale of inertial dynamics and can be expressed as
$\xi=\gamma M_{s}\alpha\tau$ to be consistent with Eq. (2) and Ref. Neeraj
_et al._ (2021).

In order to investigate the energy aspects of the inertial magnetization
dynamics, we dot-multiply Eq. (16) by the term in brackets, so that the right-
hand side vanishes

$\left(\bm{h}_{\mathrm{eff}}-\alpha\frac{d\bm{m}}{dt}-\xi\frac{d^{2}\bm{m}}{dt^{2}}\right)\cdot\frac{d\bm{m}}{dt}=0\,\,.$
(18)

By recalling Eq. (17) and after some algebra, we obtain

$-\frac{dg}{dt}-\bm{m}\cdot\frac{d\bm{h}_{a}}{dt}-\alpha\left|\frac{d\bm{m}}{dt}\right|^{2}-\frac{d}{dt}\left(\frac{\xi}{2}\left|\frac{d\bm{m}}{dt}\right|^{2}\right)=0\,\,.$
(19)

Recasting terms, one arrives at the following equation for the energy balance
in iLLG dynamics

$\frac{d}{dt}\left(g+\frac{\xi}{2}\left|\frac{d\bm{m}}{dt}\right|^{2}\right)=-\bm{m}\cdot\frac{d\bm{h}_{a}}{dt}-\alpha\left|\frac{d\bm{m}}{dt}\right|^{2}\,\,.$
(20)

It is apparent that, for constant-in-time applied field, and when the Gilbert
damping is zero ($\alpha=0$), the magnetization dynamics preserves the
quantity:

$\tilde{g}(\bm{m},d_{t}\bm{m},\bm{h}_{a})=g(\bm{m},\bm{h}_{a})+k(d_{t}\bm{m})\,,$
(21)

where

$\displaystyle
k(d_{t}\bm{m})=\frac{\xi}{2}\left|\frac{d\bm{m}}{dt}\right|^{2}\,,$ (22)

and where $d_{t}\bm{m}$ is a compact notation for $d\bm{m}/dt$. We note that
the form of $\tilde{g}$ is analogous to the the total mechanical energy
(potential + kinetic). For later use, we also introduce the following total
energy defined by using, as potential energy, the Helmholtz free energy
$a(\bm{m})$:

$\tilde{a}(\bm{m},d_{t}\bm{m})=a(\bm{m})+k(d_{t}\bm{m})\,,$ (23)

where we recall that $g(\bm{m},\bm{h}_{a})=a(\bm{m})-\bm{m}\cdot\bm{h}_{a}$.
On the basis of the the above considerations, we can conclude that Eq. (16)
with $\alpha=0$ describes a conservative dynamical system with energy equal to
$\tilde{g}(\bm{m},d_{t}\bm{m},\bm{h}_{a})$.

###  C.2 Energy dissipation during magnetization switching for LLG and iLLG
dynamics

In order to compare the work to be done by the external field during
magnetization switching in the presence or absence of inertia, one has to use
Eq. (20). Thus, the rate of change of the free energy $g$ during the switching
process becomes

$\frac{d\tilde{g}}{dt}=-\bm{m}\cdot\frac{d\bm{h}_{a}}{dt}-\alpha\left|\frac{d\bm{m}}{dt}\right|^{2}\,.$
(24)

By using in Eq. (24) the expression of $g$ in terms of Helmholtz free energy,
one can be easily convinced that Eq. (24) can be recast in a form which allows
the interpretation of the involved terms in the framework of the laws of
thermodynamics, i.e.

$\frac{da}{dt}=-\alpha\left|\frac{d\bm{m}}{dt}\right|^{2}-\frac{d}{dt}\left(\frac{\xi}{2}\left|\frac{d\bm{m}}{dt}\right|^{2}\right)+\bm{h}_{a}\cdot\frac{d\bm{m}}{dt}\,\,.$
(25)

Let us now consider the work performed by the field during a generic evolution
of magnetization in the interval $[t_{1},t_{2}]$. By integrating Eq. (25) with
respect to time in that interval, we obtain

$\displaystyle\Delta
W=\int_{t_{1}}^{t_{2}}\bm{h}_{a}\cdot\frac{d\bm{m}}{dt}dt=\Delta
a+\left[\frac{\xi}{2}\left|\frac{d\bm{m}}{dt}\right|^{2}\right]_{t_{2}}$
$\displaystyle-\left[\frac{\xi}{2}\left|\frac{d\bm{m}}{dt}\right|^{2}\right]_{t_{1}}+\int_{t_{1}}^{t_{2}}\alpha\left|\frac{d\bm{m}}{dt}\right|^{2}dt\,,$
(26)

where

$\displaystyle\Delta{a}=a(\bm{m}(t_{2}))-a(\bm{m}(t_{1}))\,,$ (27)

and where we have denoted with $\Delta W$ the work performed on the magnetic
particle during the process.

In computing $\Delta W$ for a switching process, it is convenient to use the
initial time instant $t_{0}$, i.e. the instant before which the applied
magnetic field is zero. At $t_{0}$ the magnetization is in an equilibrium
point $\bm{m}_{0}$, which corresponds to a minimum of the Helmholtz free
energy $a(\bm{m})$ and $d\bm{m}/dt=0$. Now, if $t_{p}$ is the time instant
when the external field pulse is again zero after reaching a finite value, one
could choose $t_{2}=t_{p}$, and compute $\Delta W$ as

$\displaystyle\Delta
W=\int_{t_{0}}^{t_{p}}\bm{h}_{a}\cdot\frac{d\bm{m}}{dt}dt\,.$ (28)

This computation has the difficulty that from Eq. (26) one is required to know
the values of $\bm{m}$ and $d\bm{m}/dt$ in both $t_{0}$ and $t_{p}$. A
simplified formula can be obtained by taking into account that

$\displaystyle\Delta
W=\int_{t_{0}}^{t_{p}}\bm{h}_{a}\cdot\frac{d\bm{m}}{dt}dt\approx\int_{t_{0}}^{\infty}\bm{h}_{a}\cdot\frac{d\bm{m}}{dt}dt\,,$
(29)

assuming that the applied field is negligible for $t>t_{p}$, which is a good
approximation for a Gaussian-like magnetic field pulse if one takes
$t_{p}\approx 6\sigma$. In addition, in the typical switching process, the
asymptotic magnetization state ($t\rightarrow\infty$) and the initial
magnetization state are rest states (with $d\bm{m}/dt=0$) and they have the
same value of the Helmholtz free energy. By using this fact in Eq. (26) and
Eq. (29), we arrive to the equation

$\displaystyle\Delta
W=\int_{t_{0}}^{t_{p}}\bm{h}_{a}\cdot\frac{d\bm{m}}{dt}dt\approx\int_{t_{0}}^{\infty}\alpha\left|\frac{d\bm{m}}{dt}\right|^{2}dt\,\,.$
(30)

Eq. (30) is the dimensionless version of Eq. (4) used for computations in this
paper.

An analogous result can be obtained in cyclic transformations starting and
ending in the same magnetization state (e.g. as in the case of an hysteresis
loop). In this case we obtain

$\displaystyle 0$
$\displaystyle=-\oint\left[\alpha\left|\frac{d\bm{m}}{dt}\right|^{2}+\frac{d}{dt}\left(\frac{\xi}{2}\left|\frac{d\bm{m}}{dt}\right|^{2}\right)\right]dt+\oint\bm{h}_{a}\cdot\frac{d\bm{m}}{dt}dt\,,$
(31) $\displaystyle\Delta W$
$\displaystyle=\oint\alpha\left|\frac{d\bm{m}}{dt}\right|^{2}dt\geq 0\,\,,$
(32)

The right hand side represents the amount of energy dissipated into heat
through intrinsic entropy production.

## References

  * Mallinson (2000) J. Mallinson, IEEE Transactions on Magnetics 36, 1976 (2000).
  * Bertotti  _et al._ (2003a) G. Bertotti, I. Mayergoyz, C. Serpico, and M. Dimian, Journal of applied physics 93, 6811 (2003a).
  * Stöhr and Siegmann (2006) J. Stöhr and H. C. Siegmann, Solid-State Sciences. Springer, Berlin, Heidelberg 5 (2006).
  * Gilbert (2004) T. L. Gilbert, IEEE transactions on magnetics 40, 3443 (2004).
  * Gerrits  _et al._ (2002) T. Gerrits, H. Van Den Berg, J. Hohlfeld, L. Bär, and T. Rasing, Nature 418, 509 (2002).
  * Kaka and Russek (2002) S. Kaka and S. E. Russek, Applied Physics Letters 80, 2958 (2002).
  * Hiebert  _et al._ (2002) W. Hiebert, G. Ballentine, and M. Freeman, Physical Review B 65, 140404 (2002).
  * Tudosa  _et al._ (2004) I. Tudosa, C. Stamm, A. Kashuba, F. King, H. Siegmann, J. Stöhr, G. Ju, B. Lu, and D. Weller, Nature 428, 831 (2004).
  * Kimel  _et al._ (2009) A. Kimel, B. Ivanov, R. Pisarev, P. Usachev, A. Kirilyuk, and T. Rasing, Nature Physics 5, 727 (2009).
  * Serpico  _et al._ (2004) C. Serpico, M. d’Aquino, G. Bertotti, and I. D. Mayergoyz, Journal of Applied Physics 95, 7052 (2004).
  * Böttcher  _et al._ (2011) D. Böttcher, A. Ernst, and J. Henk, Journal of Physics: Condensed Matter 23, 296003 (2011).
  * Bhattacharjee  _et al._ (2012a) S. Bhattacharjee, L. Nordström, and J. Fransson, Physical review letters 108, 057204 (2012a).
  * Bastardis  _et al._ (2018) R. Bastardis, F. Vernay, and H. Kachkachi, Physical Review B 98, 165444 (2018).
  * Makhfudz  _et al._ (2020) I. Makhfudz, E. Olive, and S. Nicolis, Applied Physics Letters 117, 132403 (2020).
  * Neeraj  _et al._ (2021) K. Neeraj, N. Awari, S. Kovalev, D. Polley, N. Z. Hagström, S. S. P. K. Arekapudi, A. Semisalova, K. Lenz, B. Green, J.-C. Deinert, _et al._ , Nature Physics 17, 245 (2021).
  * Ciornei  _et al._ (2011) M.-C. Ciornei, J. Rubí, and J.-E. Wegrowe, Physical Review B 83, 020410 (2011).
  * Olive  _et al._ (2012) E. Olive, Y. Lansac, and J.-E. Wegrowe, Applied Physics Letters 100, 192407 (2012).
  * Fähnle  _et al._ (2011) M. Fähnle, D. Steiauf, and C. Illg, Physical Review B 84, 172403 (2011).
  * Mondal  _et al._ (2017) R. Mondal, M. Berritta, A. K. Nandy, and P. M. Oppeneer, Physical Review B 96, 024425 (2017).
  * Devolder  _et al._ (2006) T. Devolder, H. W. Schumacher, and C. Chappert, Spin Dynamics in Confined Magnetic Structures III , 1 (2006).
  * d’Aquino (2004) M. d’Aquino, http://wpage. unina. it/mdaquino/PhD_thesis/main/main. html (2004).
  * Bauer  _et al._ (2000) M. Bauer, J. Fassbender, B. Hillebrands, and R. Stamps, Physical Review B 61, 3410 (2000).
  * d’Aquino  _et al._ (2005) M. d’Aquino, C. Serpico, G. Miano, I. D. Mayergoyz, and G. Bertotti, Journal of Applied Physics 97, 10E319 (2005).
  * Bertotti  _et al._ (2003b) G. Bertotti, I. Mayergoyz, C. Serpico, and M. d’Aquino, IEEE Transactions on Magnetics 39, 2501 (2003b).
  * Bhattacharjee  _et al._ (2012b) S. Bhattacharjee, A. Bergman, A. Taroni, J. Hellsvik, B. Sanyal, and O. Eriksson, Physical Review X 2, 011013 (2012b).
  * Wienholdt  _et al._ (2012) S. Wienholdt, D. Hinzke, and U. Nowak, Physical review letters 108, 247207 (2012).
  * Shutyi and Sementsov (2020) A. Shutyi and D. Sementsov, Physics of the Solid State 62, 1850 (2020).
  * Serpico  _et al._ (2009) C. Serpico, M. d’Aquino, G. Bertotti, and I. D. Mayergoyz, IEEE Transactions on Magnetics 45, 5224 (2009).
  * d’Aquino  _et al._ (2004) M. d’Aquino, W. Scholz, T. Schrefl, C. Serpico, and J. Fidler, Journal of Applied Physics 95, 7055 (2004).
  * Nozaki and Matsuyama (2006) Y. Nozaki and K. Matsuyama, Journal of applied physics 100, 053911 (2006).
  * Bazaliy (2011) Y. B. Bazaliy, Journal of Applied Physics 110, 063920 (2011).
  * Polley  _et al._ (2018) D. Polley, M. Pancaldi, M. Hudl, P. Vavassori, S. Urazhdin, and S. Bonetti, Journal of Physics D: Applied Physics 51, 084001 (2018).
  * Sederberg  _et al._ (2020) S. Sederberg, F. Kong, and P. B. Corkum, Physical Review X 10, 011063 (2020).
  * Jhuria  _et al._ (2020) K. Jhuria, J. Hohlfeld, A. Pattabi, E. Martin, A. Y. A. Córdova, X. Shi, R. L. Conte, S. Petit-Watelot, J. C. Rojas-Sanchez, G. Malinowski, _et al._ , Nature Electronics 3, 680 (2020).
  * Thonig  _et al._ (2017) D. Thonig, O. Eriksson, and M. Pereiro, Scientific reports 7, 1 (2017).
  * Giordano and Déjardin (2020) S. Giordano and P.-M. Déjardin, Physical Review B 102, 214406 (2020).
  * Mondal  _et al._ (2021) R. Mondal, S. Großenbach, L. Rózsa, and U. Nowak, Physical Review B 103, 104404 (2021).
  * Mondal (2021) R. Mondal, Journal of Physics: Condensed Matter 33, 275804 (2021).
  * Titov  _et al._ (2021a) S. Titov, W. Coffey, Y. P. Kalmykov, M. Zarifakis, and A. Titov, Physical Review B 103, 144433 (2021a).
  * Thibaudeau and Nicolis (2021) P. Thibaudeau and S. Nicolis, arXiv preprint arXiv:2103.04787 (2021).
  * Lomonosov  _et al._ (2021) A. M. Lomonosov, V. V. Temnov, and J.-E. Wegrowe, arXiv preprint arXiv:2105.07376 (2021).
  * Ruggeri (2021) M. Ruggeri, arXiv preprint arXiv:2103.09888 (2021).
  * Cherkasskii  _et al._ (2021) M. Cherkasskii, M. Farle, and A. Semisalova, arXiv preprint arXiv:2103.07864 (2021).
  * Rahman and Bandyopadhyay (2021) R. Rahman and S. Bandyopadhyay, Journal of Physics: Condensed Matter (2021).
  * Titov  _et al._ (2021b) S. Titov, W. Coffey, Y. P. Kalmykov, and M. Zarifakis, Physical Review B 103, 214444 (2021b).
  * Anders  _et al._ (2020) J. Anders, C. Sait, and S. Horsley, arXiv preprint arXiv:2009.00600 (2020).
  * Lou  _et al._ (2021) P. C. Lou, A. Katailiha, R. G. Bhardwaj, W. P. Beyermann, D. M. Juraschek, and S. Kumar, Nano Letters 21, 2939 (2021).
  * Gupta  _et al._ (2021) R. Gupta, S. Husain, A. Kumar, R. Brucas, A. Rydberg, and P. Svedlindh, Advanced Optical Materials 9, 2001987 (2021).


\vdots&\ddots&\vdots\\\
x_{m,1}&\dots&x_{m,n}\end{bmatrix}\right)&\coloneqq\begin{bmatrix}\frac{\exp(x_{1,1})}{\sum_{i=j}^{n}\exp(x_{1,j})}&\dots&\frac{\exp(x_{1,n})}{\sum_{j=1}^{n}\exp(x_{1,j})}\\\
\vdots&\ddots&\vdots\\\
\frac{\exp(x_{m,1})}{\sum_{i=1}^{m}\exp(x_{m,j})}&\dots&\frac{\exp(x_{m,n})}{\sum_{j=1}^{n}\exp(x_{m,j})}\end{bmatrix}.\\\
\end{split}$

##  Appendix J Additional Empirical Results

In this section, we provide additional empirical results for the degenerate
parabolic Generalized Porous Medium (GPME) family of conservation laws as well
as for hyperbolic conservation laws.

###  J.1 GPME Family of Conservation Laws

Here, we include additional solution profiles and conservation profiles over
time for the GPME family of equations, ranging from the “easy" diffusion
(heat), “medium" PME, to the “hard" Stefan equations.

####  J.1.1 Diffusion (Heat) Equation

Figure 9:  Solution profiles for the diffusion (heat) equation at time $t=0.5$
for diffusivity (conductivity) test-time parameter $k=1$ in the top row and
$k=5$ in the bottom row. Each model is trained on samples of
$k\in\mathcal{A}=[1,5]$. The shaded region illustrates $\pm 3$ standard
deviation uncertainty intervals. ProbConserv-ANP and HardC-ANP both display
tighter uncertainty bounds than the baseline ANP, while SoftC-ANP is more
diffuse. The uncertainty is relatively homoscedastic on this “easy” case.

##### Solution profiles.

Figure 9 shows the solution profiles for the “easy" diffusion equation, at
time $t=0.5$, where a sine curve is damped over time for test-time parameter
$k=1,5\in\mathcal{A}=[1,5]$. Table 9 shows the corresponding metrics.

Table 9:  Mean and standard error for CE $\times 10^{-3}$ (should be zero), LL (higher is better) and MSE $\times 10^{-4}$ (lower is better) over $n_{\text{test}}=50$ for the (“easy") diffusion equation at time $t=0.5$ with variable diffusivity constant $k$ parameter in the range $\mathcal{A}=[1,5]$ and test-time parameter values $k=1,5$.  | $k=1$ |  |  | $k=5$ |  |
---|---|---|---|---|---|---
| CE | LL | MSE | CE | LL | MSE
ANP | 4.68 (0.10) |  $2.72$ (0.02) |  $1.71$ (0.41) | 1.76 (0.04) |  $3.28$ (0.02) |  $0.547$ (0.08)
SoftC-ANP | 3.47 (0.17) |  $2.40$ (0.02) |  $2.24$ (0.78) | 2.86 (0.05) |  $2.83$ (0.02) |  $1.75$ (0.24)
HardC-ANP |  0 (0.00) |  3.08 (0.04) |  $\mathbf{1.37}$ (0.33) |  0 (0.00) |  3.64 (0.03) |  $\mathbf{0.461}$ (0.07)
ProbConserv-ANP |  0 (0.00) |  2.74 (0.02) |  1.55 (0.33) |  0 (0.00) |  3.30 (0.02) |  0.485 (0.07)

####  J.1.2 Porous Medium Equation (PME)

##### Results for different $\lambda$ for SoftC-ANP.

As is the case with PINNs (Raissi et al., 2019), the SoftC-ANP method has a
hyper-parameter $\lambda$ that controls the balance in the training loss
between the reconstruction and differential term. A higher value of $\lambda$
places more emphasis on the residual of the PDE term and less emphasis on the
evidence lower bound (ELBO) from the ANP.

To investigate whether tuning $\lambda$ will lead to significantly different
results, we report results for different values of $\lambda$ for the SoftC-ANP
on the Porous Medium Equation (PME). Since these results are presented on the
same test dataset used in Table 2, it provides an optimistic case on how
tuning $\lambda$ could improve the results for SoftC-ANP. Table 10 shows that
the predictive performance is roughly the same across different values of
$\lambda$, with both MSE and LL worse than the original ANP across the board
and the conservation error (CE) $G\mu-b$ at the final time worse for $m=6$.

Table 10: Investigation of the effect of the soft constraint penalty parameter $\lambda$ in the SoftC-ANP baseline. The metrics CE $\times 10^{-3}$ (should be zero), LL (higher is better) and MSE $\times 10^{-4}$ (lower is better) are reported for the (“medium") PME at time $t=0.5$ with variable $m$ parameter in the range $\mathcal{A}=[0.99,6]$ and test-time parameters $m\in\\{1,3,6\\}$. We see that the performance is not significantly changed as a function of $\lambda$, and, surprisingly, that the unconstrained ANP ($\lambda=0)$ performs better in most metrics than SoftC-ANP. | $m=1$ |  |  | $m=3$ |  |  | $m=6$ |  |
---|---|---|---|---|---|---|---|---|---
| CE | LL | MSE | CE | LL | MSE | CE | LL | MSE
ANP $(\lambda=0)$ | $6.67$ | $\bf{3.49}$ | $\bf{0.94}$ | $-1.23$ | $\bf{3.67}$ | $\bf{1.90}$ | $\bf{-2.58}$ | $\bf{3.81}$ | $\bf{7.62}$
SoftC-ANP $(\lambda=0.01)$ | $5.58$ | $3.11$ | $1.11$ | $-0.61$ | $3.46$ | $2.03$ | $-3.00$ | $3.49$ | $7.76$
SoftC-ANP $(\lambda=0.1)$ | $5.58$ | $3.11$ | $1.11$ | $-0.67$ | $3.46$ | $2.07$ | $-3.01$ | $3.49$ | $7.87$
SoftC-ANP $(\lambda=1)$ | $5.62$ | $3.11$ | $1.11$ | $-0.65$ | $3.46$ | $2.06$ | $-3.03$ | $3.49$ | $7.82$
SoftC-ANP $(\lambda=10)$ | $\bf{5.52}$ | $3.11$ | $1.08$ | $\bf{-0.56}$ | $3.46$ | $2.04$ | $-3.02$ | $3.49$ | $7.76$
SoftC-ANP $(\lambda=100)$ | $5.62$ | $3.11$ | $1.11$ | $-0.59$ | $3.46$ | $2.03$ | $-3.03$ | $3.49$ | $7.69$

#####  ProbConserv-ANP with diffusion.

As described in subsection D.2, we explore adding numerical diffusion for
eliminating artificial small-scale noises when enforcing conservation. Table
11 shows that adding artificial diffusion improves both MSE and LL compared to
the conservation constraint alone. Figures 10-11 illustrate that by removing
such artificial noises, ProbConserv-ANP with diffusion leads to tighter
uncertainty bounds as well as higher LL than the other baselines.

Table 11: Mean and standard error for CE $\times 10^{-3}$ (should be zero), LL
(higher is better) and MSE $\times 10^{-4}$ (lower is better) over
$n_{\text{test}}=50$ runs for the (“medium”) PME at time $t=0.5$ with variable
$m$ parameter in the range $\mathcal{A}=[0.99,6]$. We see that ProbConserv-ANP
(w/diff) improves the performance on ProbConserv-ANP by applying smoothing at
the sharp boundary as the test-time parameter $m$ is increased.

| $m=1$ |  |  | $m=3$ |  |  | $m=6$ |  |
---|---|---|---|---|---|---|---|---|---
| CE | LL | MSE | CE | LL | MSE | CE | LL | MSE
ANP |  $6.67$ (0.39) |  $3.49$ (0.01) |  $0.94$ (0.09) |  $-1.23$ (0.29) |  $3.67$ (0.00) |  $1.90$ (0.04) |  $-2.58$ (0.23) |  $3.81$ (0.01) |  7.67 (0.09)
SoftC-ANP |  $5.62$ (0.35) |  $3.11$ (0.01) |  $1.11$ (0.14) |  $-0.65$ (0.30) |  $3.46$ (0.00) |  $2.06$ (0.03) |  $-3.03$ (0.26) |  $3.49$ (0.00) |  $7.82$ (0.09)
HardC-ANP |  0 (0.00) | 3.16 (0.04) | 0.43 (0.04) |  0 (0.00) |  $3.44$ (0.03) | 1.86 (0.03) |  0 (0.00) | 3.40 (0.05) |  7.61 (0.09)
ProbConserv-ANP |  0 (0.00) | 3.56 (0.01) |  0.17 (0.02) |  0 (0.00) | 3.68 (0.00) | 2.10 (0.07) |  0 (0.00) | 3.83 (0.01) | 10.4 (0.04)
ProbConserv-ANP (w/diff) |  0 (0.00) |  4.04 (0.02) |  0.15 (0.02) |  0 (0.00) |  3.96 (0.00) |  1.43 (0.05) |  0 (0.00) |  4.03 (0.01) | 7.91 (0.03)

##### Solution and error profiles.

Figures 10-11 illustrate the differing solution profiles and errors for the
PME for various values of $m\in\\{1,3,6\\}$, respectively. As expected, we see
a gradient for $m>1$ that becomes sharper and approaches infinity for $m=6$.
Increasing $m$ results in smaller values of the PDE parameter denoting the
pressure $k(u)=u^{m}$, which increases the degeneracy for smaller values of
$k(u)$, i.e., larger values of $m$. In this case the problem also becomes more
challenging. For $m=1$, we have a piecewise linear solution, and for $m=3,6$
we see sharper oscillatory uncertainty bounds at the front or free boundary,
resulting in some negative values at this boundary as well. We see the value
of the uncertainty quantification to reflect that the model is certain in the
parabolic regions to the left and right of the sharp boundary especially in
the zero (degeneracy) region, and is most uncertain at the boundary
(degeneracy) point.

Figure 10: Solution profiles and uncertainty intervals for the PME predicted
by our ProbConserv-ANP and other baselines. The solutions are obtained for
three scenarios with increasing sharpness in the profile as $m$ is increased
from $m=1$ to $m=6$ from left to right, respectively. The HardC-ANP model,
which assumes constant variance for the whole domain, results in too high
uncertainty in the zero (degenerate) region, unlike our proposed ProbConserv-
ANP approach that incorporates the variance information to effectively handle
this heteroscedasticity. Adding diffusion to ProbConserv-ANP removes the
oscillations locally at the degeneracy, as desired. Figure 11: The solution
errors with uncertainty bounds as a function of $x$ for our ProbConserv-ANP
and other baselines for the PME with parameter $m\in\\{1,3,6\\}$ after
training on $m\in\mathcal{A}=[0.99,6]$. The shaded region indicates $\pm 3$
standard deviations as estimated by each model. For $m=1$, both ProbConserv-
ANP with and without diffusion result in solutions with smaller errors. While
HardC-ANP model reduces the error scale, it underestimates the zero portion of
the solution, which is nonphysical, as the solution quantity cannot be
negative. For $m\in\\{3,6\\}$, while the error magnitude becomes dominant at
the shock position for all methods, ProbConserv-ANP with diffusion provides
the lowest errors with the tightest confidence interval.

####  J.1.3 Stefan

Figure 12: True mass over time for each model. The true mass conservation
profile over time is matched exactly by our proposed method ProbConserv-ANP
and the hard-constrained HardC-ANP by design. The unconstrained ANP and
surprisingly even the differential form soft-constrained SoftC-ANP have a non-
physical linear mass profile over time.

Figure 12 shows ProbConserv-ANP follows the true profile of conserved mass in
the system over time by design. We also see that the unconstrained ANP and
surprisingly the soft-constrained SoftC-ANP that applied the differential form
as a soft constraint does not result in conservation being satisfied since it
does not enforce it exactly. For these baselines, the mass profile over time
is linear and does not match the true profile which is proportional to
$\sqrt{t}$.

###  J.2 Hyperbolic Equations

Here, we demonstrate that our approach ProbConserv-ANP also works for
hyperbolic conservation laws by considering the linear advection problem
(“medium”) and Burgers’ equation (“hard”), which are both introduced in Table
5 in Appendix C.

####  J.2.1 Linear Advection

Figure 13 displays the system total mass, $U(t)=\int_{\Omega}u(t,x)d\Omega$ as
a function of time, obtained by our ProbConserv-ANP model and the other
baselines and compared against the true curve. The results are obtained for
two values of test-time parameter $\beta=1,3$ denoting the velocity with
training range $\beta\in\mathcal{A}=[1,5]$. The unconstrained ANP contradicts
the system true mass at all times including $t=0$. By proper incorporation of
the conservation constraint, both ProbConserv-ANP and HardC-ANP methods are
able to predict the system mass and capture the actual trend over time exactly
while the soft-constrained differential form SoftC-ANP baseline results in
little improvement.

Figure 14 shows the predicted solution profiles and corresponding uncertainty
intervals for time $t=0.1$ and test-time parameter $\beta=1,3$. Our
ProbConserv-ANP model predicts sharper shock profile centered around the
actual shock position. On the contrary, both ANP and HardC-ANP lead to highly
diffusive profiles which are shifted toward the left of actual shock
interface, leading to the under-estimation of the shock position on this
downstream task. This under-estimation becomes more evident in Figure 15,
which indicates the corresponding histograms of shock position. The histograms
associated with the ANP and HardC-ANP models are skewed to the left and both
result in the averaged shock positions which are lower than the actual value
depicted by the solid vertical line. By proper leveraging of our finite volume
based physical constraint, our ProbConserv-ANP results in proper uncertainty
quantification which leads to accurate prediction of shock location compared
to the other baseline models. Table 12 also shows this accuracy improvement
with a maximum improvement of $2.86\times$ in MSE for $\beta=1$.

Figure 13: System total mass as a function of time $t$ for the linear
advection problem with test-time parameter $\beta=1,3$ and training parameter
range $\mathcal{A}=[1,5]$. Both ProbConserv-ANP and HardC-ANP satisfy
conservation of mass while the unconstrained ANP and soft-constrained SoftC-
ANP baselines deviate from the actual trend completely at all times. Figure
14: Solution profiles and uncertainty intervals for linear advection problem
at time $t=0.1$ for test-time parameter $\beta=1,3$ and training parameter
range $\mathcal{A}=[1,5]$. Despite satisfaction of conservation constraint,
HardC-ANP predicts a highly diffusive profile and remarkable underestimation
of shock interface region especially for $\beta=1$. The prediction error is
even higher for the unconstrained ANP model which does not enforce the
conservation, and the shock interface is shifted further away from the true
solution. ProbConserv-ANP results in a sharper profile than other the
baselines and the predicted shock interface is around the actual shock
position leading to more accurate shock position estimation. Figure 15: The
histogram of shock position for the linear advection problem, computed as the
mean plus or minus 3 standard deviations. Due to the shift in the shock
interface, both the ANP and HardC-ANP models underestimate the position of the
shock, and the underestimation is more significant for $\beta=1$. The
ProbConserv-ANP model provides a histogram distributed almost symmetrically
around the true shock interface and thus leads to an accurate estimate of the
shock position. Table 12:  Mean and standard error for CE (should be zero), LL
(higher is better) and MSE $\times 10^{-2}$ (lower is better) over
$n_{\text{test}}=50$ runs for the hyperbolic linear advection problem at time
$t=0.1$ with variable $\beta$ parameter in the range $\mathcal{A}=[1,5]$.

| $\beta=1$ |  |  | $\beta=3$ |  |
---|---|---|---|---|---|---
| CE | LL | MSE | CE | LL | MSE
ANP | -0.136 (0.004) | 0.96 (0.01) | 5.72 (0.25) | 0.042 (0.003) | 0.51 (0.01) | 2.03 (0.01)
SoftC-ANP | -0.137 (0.004) |  1.58 (0.03) | 7.64 (0.34) | 0.013 (0.003) |  2.31 (0.02) | 2.87 (0.20)
HardC-ANP |  0 (0.00) | -2.96 (0.34) | 4.59 (0.17) |  0 (0.00) | 1.34 (0.21) | 1.93 (0.07)
ProbConserv-ANP |  0 (0.00) | 1.06 (0.01) |  2.00 (0.06) |  0 (0.00) | 0.52 (0.01) |  1.62 (0.01)

####  J.2.2 Burgers’ Equation

Figure 16 illustrates that the total mass is linear over time, and in this
case is approximately satisfied by our ProbConserv-ANP and the baselines.
Figure 17 shows the waiting time phenomenon, where the piecewise linear
initial condition self-sharpens until the breaking time $t_{b}=1/a$, where it
forms a rightward moving shock. We see that the breaking time is inversely
proportional to the slope, and that the shock forms sooner for larger values
of $a$.

Figure 16: True mass as a function of time $t$ for the Burgers’ equation with
test-time parameter $a=1,3$ and training parameter range $\mathcal{A}=[1,4]$.
Figure 17: Solution profiles and uncertainty intervals for Burgers’ equation
at time $t=0.5$ for test-time parameter $a=1,3$ and training parameter range
$\mathcal{A}=[1,4]$.