# Tackling the Curse of Dimensionality in Fractional and Tempered Fractional
PDEs with Physics-Informed Neural Networks

Zheyuan Hu Department of Computer Science, National University of Singapore,
Singapore, 119077 (e0792494<EMAIL_ADDRESS>Kenji
Kawaguchi11footnotemark: 1 Zhongqiang Zhang Department of Mathematical
Sciences, Worcester Polytechnic Institute, Worcester, MA 01609 USA
<EMAIL_ADDRESS>George Em Karniadakis Division of Applied Mathematics, Brown
University, Providence, RI 02912, USA (george_karniadakis@brown.edu)Advanced
Computing, Mathematics and Data Division, Pacific Northwest National
Laboratory, Richland, WA, United States

###### Abstract

Fractional and tempered fractional partial differential equations (PDEs) are
effective models of long-range interactions, anomalous diffusion, and non-
local effects. Traditional numerical methods for these problems are mesh-
based, thus struggling with the curse of dimensionality (CoD). Physics-
informed neural networks (PINNs) offer a promising solution due to their
universal approximation, generalization ability, and mesh-free training. In
principle, Monte Carlo fractional PINN (MC-fPINN) estimates fractional
derivatives using Monte Carlo methods and thus could lift CoD. However, this
may cause significant variance and errors, hence affecting convergence; in
addition, MC-fPINN is sensitive to hyperparameters. In general, numerical
methods and specifically PINNs for tempered fractional PDEs are under-
developed. Herein, we extend MC-fPINN to tempered fractional PDEs to address
these issues, resulting in the Monte Carlo tempered fractional PINN (MC-
tfPINN). To reduce possible high variance and errors from Monte Carlo
sampling, we replace the one-dimensional (1D) Monte Carlo with 1D Gaussian
quadrature, applicable to both MC-fPINN and MC-tfPINN. We validate our methods
on various forward and inverse problems of fractional and tempered fractional
PDEs, scaling up to 100,000 dimensions. Our improved MC-fPINN/MC-tfPINN using
quadrature consistently outperforms the original versions in accuracy and
convergence speed in very high dimensions.

##  1 Introduction

Fractional partial differential equations (PDEs) extend the concept of
differentiation to non-integer orders, allowing more flexibility in modeling
complex phenomena. They have found applications in various fields, such as
physics, biology, finance, and engineering. Fractional PDEs usually involve
time-fractional derivatives or fractional Laplacians. These non-integer order
operators capture non-stadard behavior such as long-range interactions,
anomalous diffusion, and scale invariance, which are not adequately described
by integer-order derivatives. Furthermore, tempered fractional PDEs introduce
a tempering parameter that relaxes the long-range interactions of the
fractional operator. This tempering parameter allows a transition between
standard integer-order PDEs and fully non-local fractional PDEs. Tempered
fractional PDEs are particularly useful for modeling phenomena of a balance
between local interactions and non-local effects. By tuning the tempering
parameter, one can adjust the influence of non-locality in the system, making
tempered fractional PDEs versatile tools for capturing a wide range of
behaviors in various applications.

Given their importance in science and engineering, numerous methods have been
proposed to obtain their numerical solutions, see e.g. in [8, 19] and many
others. In addition to traditional mesh-based methods, Monte Carlo methods are
also developed for fractional operators, e.g., [22, 32]. The physics-informed
neural network (PINN) [30] has emerged as a mesh-free approach, benefiting
primarily from the powerful generalization capabilities of neural networks
[20], their universal approximation properties [2], and robust optimization
using stochastic gradients [21]. Recently, PINN has been validated as highly
effective in solving high-dimensional PDEs [16, 15, 18], lifting the curse of
dimensionality. Regarding fractional PDEs, fractional PINN (fPINN) [29] and
Monte Carlo fPINN (MC-fPINN) [12] represent seminal works for solving
fractional PDEs with PINNs. Despite demonstrating high accuracy and
scalability across various low-dimensional and high-dimensional fractional PDE
forward/inverse problems, their application to particularly high dimensions,
such as problems in thousands of dimensions, remains limited. Moreover, the
PINN approach is lacking for the equally significant and more practical
tempered fractional PDEs.

Furthermore, fPINN [29] suffers from the curse of dimensionality as it uses a
grid-based method to discretize the fractional derivatives. Although MC-fPINN
[12] employs Monte Carlo to overcome this issue, the computation of the
fractional derivative’s singularity-integrated integral using Monte Carlo may
lead to significant variance and errors, hindering PINN convergence. Moreover,
MC-fPINN relies on many hyperparameters, such as setting a radius to divide
the integral of the fractional derivative into internal and external parts,
and it further requires a truncation threshold for the former integral to
avoid numerical overflow. Therefore, the accuracy of MC-fPINN is highly
sensitive to these hyperparameters, thereby limiting its practical
applicability. Additionally, while recent efforts have scaled up PINN to
thousands of dimensions [16, 15], these endeavors primarily address integer-
order PDEs and do not apply to the fractional-order PDEs considered in this
work.

To address these issues, this paper makes two core contributions: (1) We
extend MC-fPINN [12] to tempered fractional PDEs, resulting in Monte Carlo
tempered fractional PINN, denoted as MC-tfPINN. (2) Furthermore, to overcome
the high variance and error associated with Monte Carlo in MC-fPINN, as well
as its dependency on many hyperparameters, we propose replacing the 1D Monte
Carlo computation for singularity-integrated parts with a 1D Gaussian
quadrature. This quadrature method is also applied to MC-tfPINN, resulting in
improved versions.

Specifically, MC-fPINN primarily relies on Monte Carlo sampling from a Beta
distribution. For tempered fractional PDEs, we use Monte Carlo sampling from a
Gamma distribution. Additionally, due to the practicality introduced by the
tempering effect, MC-tfPINN does not need to split the integral of the
tempered fractional derivative into two parts for calculation, making its
algorithm more concise.

For the second part of improving MC-fPINN and MC-tfPINN, we observe that
simulating a $d$-dimensional fractional PDE operator requires computing a
$d$-dimensional integral, which can be reformulated as a product integral of a
1D radial singular integral and an integral over a $d$-dimensional unit
sphere. For the former, the singularity of the integral introduces significant
variance and error to Monte Carlo, leading us to propose using Gaussian
quadrature to compute this 1D radial singular integral. Compared to Monte
Carlo sampling, integration using the Gauss quadrature provides theoretical
guarantees of accuracy when computing singular integrals, making it more
precise. Additionally, with Monte Carlo, resampling at the beginning of each
optimization epoch is required, which requires additional time, whereas
quadrature points remain fixed and do not require resampling. Moreover, as
quadrature points are fixed, we do not need a truncation threshold to avoid
numerical overflow. In the Monte Carlo method, such a threshold is necessary
due to the randomness of the sampling to prevent sampling extreme values.
Furthermore, regarding MC-fPINN’s need to set a radius to divide the integral
of the fractional derivative into internal and external parts, we choose this
radius as the diameter of the PDE support domain, thereby reducing the
dependence on hyperparameters and simplifying the computation in the second
part, thus further accelerating PINN.

We validate the effectiveness of our proposed methods on several fractional
and tempered fractional PDEs, scaling up to 100,000 dimensions, far exceeding
previous works. Firstly, for the fractional Poisson equation, we test a
complex anisotropic solution, confirming the higher accuracy and faster
training speed of the improved MC-fPINN over the regular MC-fPINN. Next, we
compare MC-tfPINN and the quadrature-based improved MC-tfPINN on tempered
fractional Poisson and diffusion equations, once again confirming better
efficiency and higher accuracy of our methods. Additionally, we vary the
fractional order and tempering coefficient in the tempered fractional PDE
operator, demonstrating the stability of our models across various dimensions.
Finally, we test inverse problems and show that our proposed models are
versatile and applicable.

In summary, the novelty and contribution of the work lies in the following
aspects: 1) extending MC-fPINN to tempered fractional PDEs; and 2)introducing
improved versions based on quadrature over Monte Carlo for fPINN, rendering
simulations of high-dimensional fractional and tempered fractional PDEs more
accurate and efficient.

##  2 Related Work

For fractional PDE solvers, Pang et al. [29] introduce a fractional Physics-
Informed Neural Network (fPINN), utilizing neural networks as proxies for
solutions and employing a discretized approach to fractional derivatives for
training. However, the discretization in fPINN is grid-based and thus suffers
from the curse of dimensionality. Guo et al. [12] develop a Monte Carlo
fractional PINN (MC-fPINN), which approximates Caputo-type time-fractional
derivatives and fractional Laplacian within a hyper-singular integral
framework, leveraging Monte Carlo methods based on Beta distributions. Ma et
al. [27] further propose bi-orthogonal fPINN to deal with stochastic
fractional PDEs. Chen et al. [3] use PINNs to address both forward and inverse
problems associated with Fokker-Planck-Lévy equations and use the
discretization scheme in [11] to estimate fractional derivatives. Firoozsalari
et al. [10] employ Gaussian quadrature for calculating numerical integrals
associated with fractional derivatives, concentrating exclusively on time
derivatives and neglecting to address scenarios involving high-dimensional
fractional Laplacians. Leonenko and Podlubny [22] introduce a Monte Carlo
method to estimate the Grünwald–Letnikov fractional derivative. Based on the
approach by Leonenko and Podlubny [22], Wang and Karniadakis [34] propose a
general quasi Monte Carlo PINN for solving fractional PDEs on irregular
domains. Sheng et al. [32] propose a novel Monte Carlo method for fractional
Laplacian in high dimensions.

For tempered fractional PDEs, a pedagogical guide for tempered fractional
calculus is provided in [31]. The pedagogical guide for the tempered factional
diffusion process is also presented in [1], where the tempered stable
distribution is discussed for the tempered Lévy process. Duo et al. [6]
propose an accurate finite difference method for tempered fractional Laplacian
and research the tempered effect. Sun et al. [33] provide algorithm
implementation and analysis for 2D tempered fractional Laplacian. Deng et al.
[5] solve boundary problems of fractional and tempered fractional PDEs. Li and
Deng [23] propose several high-order numerical schemes for the tempered
fractional diffusion equations. Hanert and Piret [13] propose a Chebyshev
pseudospectral method for the tempered fractional diffusion equations. There
also exists a spectral method for tempered fractional cases [35].

##  3 Preliminary

We consider high-dimensional tempered fractional PDEs on a bounded spatial
domain $\Omega\subset\mathbb{R}^{d}$:

$\displaystyle\frac{\partial^{\gamma,\lambda_{t}}u(\boldsymbol{x},t)}{\partial
t^{\gamma}}+c\left(-\Delta_{\lambda_{\boldsymbol{x}}}\right)^{\alpha/2}u(\boldsymbol{x},t)+\boldsymbol{v}\cdot\nabla_{\boldsymbol{x}}u(\boldsymbol{x},t)$
$\displaystyle=f(\boldsymbol{x},t),\quad(\boldsymbol{x},t)\in\Omega\times(0,T].$
(1) $\displaystyle u(\boldsymbol{x},t=0)$
$\displaystyle=g(\boldsymbol{x}),\quad\boldsymbol{x}\in\Omega.$ (2)
$\displaystyle u(\boldsymbol{x},t)$
$\displaystyle=0,\quad(\boldsymbol{x},t)\in\Omega^{c}\times(0,T].$ (3)

Here, $f(\boldsymbol{x},t)$ is the forcing term, $c\in\mathbb{R}^{+}$ is the
diffusion coefficient, $\lambda_{t}\geq 0$ and $\lambda_{\boldsymbol{x}}\geq
0$ are the tempering coefficients, and $\boldsymbol{v}\in\mathbb{R}^{d}$ is
the flow velocity. The time-tempered fractional derivative is defined as [9]:

$\displaystyle\frac{\partial^{\gamma,\lambda}f(t)}{\partial
t^{\gamma}}:=\frac{1}{\Gamma(1-\gamma)}\int_{0}^{t}\exp\left(-\lambda(t-\tau)\right)(t-\tau)^{-\gamma}\frac{\partial
f(\tau)}{\partial\tau}d\tau=e^{-\lambda
t}\frac{\partial^{\gamma}\left(e^{\lambda t}f(t)\right)}{\partial
t^{\gamma}},\quad 0<\gamma<1.$ (4)

When $\lambda=0$, the tempered fractional derivative becomes the Caputo time-
fractional derivative. The tempered fractional Laplacian is defined [5]:

$\displaystyle(-\Delta_{\lambda})^{\alpha/2}u:=C_{d,\alpha,\lambda}\text{P.V.}\int_{\mathbb{R}^{d}}\frac{u(\boldsymbol{x})-u(\boldsymbol{y})}{\exp\left(\lambda\|\boldsymbol{x}-\boldsymbol{y}\|\right)\|\boldsymbol{x}-\boldsymbol{y}\|^{d+\alpha}}d\boldsymbol{y},\quad\alpha\in(0,2),$
(5)

where P.V. denotes the principle value of the integral and
$C_{d,\alpha,\lambda}=\frac{\Gamma\left(\frac{d}{2}\right)}{2\pi^{\frac{d}{2}}\left|\Gamma\left(-{\alpha}\right)\right|}$.
Note that when $\lambda=0$, the tempered fractional Laplacian becomes the
fractional Laplacian e hyper-singular integral (13). When
$\lambda_{t}=\lambda_{\boldsymbol{x}}=0$, the equation (1) becomes the
fractional PDEs considered in [12]. For simplicity, we denote

$\displaystyle\mathcal{L}^{\phi}\left[u(\boldsymbol{x},t)\right]:=\frac{\partial^{\gamma,\lambda_{t}}u(\boldsymbol{x},t)}{\partial
t^{\gamma}}+c\left(-\Delta_{\lambda_{\boldsymbol{x}}}\right)^{\alpha/2}u(\boldsymbol{x},t)+\boldsymbol{v}\cdot\nabla_{\boldsymbol{x}}u(\boldsymbol{x},t),$
(6)

and $\phi=\\{\gamma,c,\alpha,\boldsymbol{v}\\}$ is a collection of
coefficients and parameters in the problem above.

###  3.1 Tempered fractional PDEs: forward and inverse problems

We consider forward and inverse tempered fractional PDE problems. In the
forward problem, we aim to solve the PDE solution $u(\boldsymbol{x},t)$ given
all coefficients/parameters $\phi$), the forcing term $f(\boldsymbol{x},t)$
and initial and boundary conditions. In the inverse problem, we have sparse
observation data over the exact solution $u$ while we are given initial and
boundary conditions. The goal is to find the unknown coefficients/parameters
in the PDE operator, e.g., $\lambda,\alpha,\gamma,\boldsymbol{v}$, etc.

###  3.2 Physics-Informed Neural Networks (PINNs)

We use PINNs [30] to solve both forward and inverse problems of the fractional
and tempered fractional PDEs defined in the previous subsections.

PINN is a neural network-based PDE solver that parameterizes the PDE surrogate
solution using the model $u_{\theta}(\boldsymbol{x},t)$ with parameters
$\theta$. PINN optimizes the neural network model
$u_{\theta}(\boldsymbol{x},t)$ to satisfy the PDE operator, initial, and
boundary conditions. As long as PINN approximates these three conditions well
enough, it can approximate the PDE solution accurately. See PINN’s theoretical
guarantees and foundations in [4, 14, 26, 28].

Boundary condition. We hard constraint the vanishing boundary condition using
a function $\rho(\boldsymbol{x})$, such that $\rho(\boldsymbol{x})>0$ for
$\boldsymbol{x}\in\Omega$ and $\rho(\boldsymbol{x})=0$ for
$\boldsymbol{x}\in\Omega^{c}$. i.e., the surrogate model will be

$\displaystyle u_{\theta}(\boldsymbol{x},t)=\rho(\boldsymbol{x})\cdot
u_{\theta}^{\text{NN}}(\boldsymbol{x},t),$ (7)

where $u_{\theta}^{\text{NN}}$ is a full-connected network and $u_{\theta}$
will be the final surrogate model. Here, hard constraints are essential for
stable training and successful optimization of PINNs.

Initial condition. Given the initial points
$\\{\boldsymbol{x}_{n,\text{initial}}\\}_{n=1}^{N_{\text{initial}}}\subset\Omega$,
we minimize the $L_{2}$ distance between PINN prediction and exact initial
condition:

$\displaystyle\mathcal{L}_{\text{initial}}(\theta)=\frac{1}{N_{\text{initial}}}\sum_{n=1}^{N_{\text{initial}}}{|u_{\theta}(\boldsymbol{x}_{n,\text{initial}},0)-g(\boldsymbol{x})|}^{2}.$
(8)

PDE operator. Given the residual/collocation points
$\\{\boldsymbol{x}_{n,\text{residual}},t_{n,\text{residual}}\\}_{n=1}^{N_{\text{residual}}}\subset\Omega\times(0,T]$,
we minimize the PINN residual:

$\mathcal{L}_{\text{residual}}(\theta,\phi)=\frac{1}{N_{\text{residual}}}\sum_{n=1}^{N_{\text{residual}}}|\mathcal{L}^{\phi}u_{\theta}(\boldsymbol{x}_{n,\text{residual}},t_{n,\text{residual}})-f(\boldsymbol{x}_{n,\text{residual}},t_{n,\text{residual}})|^{2}.$
(9)

In forward problems, we have full knowledge of PDE coefficients $\phi$, and
thus, we only optimize $\theta$ while in inverse problems, we optimize both
PINN model parameters $\theta$ and PDE unknown coefficients $\phi$.

Data loss. Inverse problems allow access to sparse observation over the exact
solution
$\\{\boldsymbol{x}_{n,\text{data}},t_{n,\text{data}},u^{*}_{n}\\}_{n=1}^{N_{\text{data}}}$,
which further leads to the following data loss for the PINN surrogate model to
capture these observation data to recover the unknown PDE coefficients $\phi$:

$\displaystyle\mathcal{L}_{\text{data}}(\theta)=\frac{1}{N_{\text{data}}}\sum_{n=1}^{N_{\text{data}}}{|u_{\theta}(\boldsymbol{x}_{n,\text{data}},t_{n,\text{data}})-u^{*}_{n}|}^{2}.$
(10)

In addition to these loss functions, we usually weighted sum them to form the
final loss. The weights are denoted $\lambda_{\text{initial}}$ for the initial
loss, $\lambda_{\text{residual}}$ for the residual loss, and
$\lambda_{\text{data}}$ for the data loss, respectively. Given the different
backgrounds of forward and inverse problems, we summarize their respective
loss functions for PINNs.

  * •

In the forward problem, we optimize the initial and residual losses:

$\displaystyle\min_{\theta}\lambda_{\text{initial}}\mathcal{L}_{\text{initial}}(\theta)+\lambda_{\text{residual}}\mathcal{L}_{\text{residual}}(\theta).$
(11)

  * •

In the inverse problem, we optimize all losses:

$\displaystyle\min_{\theta,\phi}\lambda_{\text{initial}}\mathcal{L}_{\text{initial}}(\theta)+\lambda_{\text{residual}}\mathcal{L}_{\text{residual}}(\theta,\phi)+\lambda_{\text{data}}\mathcal{L}_{\text{data}}(\theta).$
(12)

###  3.3 Monte Carlo fPINN (MC-fPINN)

Although PINN, combined with modern automatic differentiation libraries, can
utilize the above loss function to solve integer-order PDEs, automatic
differentiation does not cover fractional derivatives. In Guo et al. [12],
Monte Carlo fPINN is proposed to simulate the fractional Laplacian and
derivatives to enable PINN’s application to fractional PDEs.

####  3.3.1 Fractional Laplacian

In Monte Carlo-based approaches, e.g., Monte-Carlo fPINN [12], the following
fractional Laplacian is defined as the hyper-singular integral [24]:

$\displaystyle(-\Delta)^{\alpha/2}u:=C_{d,\alpha}\text{P.V.}\int_{\mathbb{R}^{d}}\frac{u(\boldsymbol{x})-u(\boldsymbol{y})}{\|\boldsymbol{x}-\boldsymbol{y}\|^{d+\alpha}}d\boldsymbol{y},\quad\alpha\in(0,2),$
(13)

where P.V. denotes the principle value of the integral and
$C_{d,\alpha}=\frac{2^{\alpha}\Gamma\left(\frac{d+\alpha}{2}\right)}{\pi^{\frac{d}{2}}\left|\Gamma\left(-\frac{\alpha}{2}\right)\right|}$.
According to [12], the integral is divided into two parts for Monte Carlo
simulation,

$\displaystyle(-\Delta)^{\alpha/2}u=C_{d,\alpha}\left(\int_{\boldsymbol{y}\in
B_{r_{0}}(\boldsymbol{x})}\frac{u(\boldsymbol{x})-u(\boldsymbol{y})}{\|\boldsymbol{x}-\boldsymbol{y}\|^{d+\alpha}}d\boldsymbol{y}+\int_{\boldsymbol{y}\notin
B_{r_{0}}(\boldsymbol{x})}\frac{u(\boldsymbol{x})-u(\boldsymbol{y})}{\|\boldsymbol{x}-\boldsymbol{y}\|^{d+\alpha}}d\boldsymbol{y}\right),$
(14)

where $B_{r_{0}}(\boldsymbol{x})=\\{\boldsymbol{y}\ |\
\|\boldsymbol{y}-\boldsymbol{x}\|\leq r_{0}\\}$. [12] shows that Monte Carlo
can simulate the two parts:

$\displaystyle\int_{\boldsymbol{y}\in
B_{r_{0}}(\boldsymbol{x})}\frac{u(\boldsymbol{x})-u(\boldsymbol{y})}{\|\boldsymbol{x}-\boldsymbol{y}\|^{d+\alpha}}d\boldsymbol{y}=\frac{|\mathcal{S}^{d-1}|r_{0}^{2-\alpha}}{2(2-\alpha)}\mathbb{E}_{\boldsymbol{\xi},r\sim
f_{1}(r)}\left[\frac{2u(\boldsymbol{x})-u(\boldsymbol{x}-\boldsymbol{\xi}r)-u(\boldsymbol{x}+\boldsymbol{\xi}r)}{r^{2}}\right],$
(15)

where we sample $r$ based on $r/r_{0}\sim\operatorname{Beta}(2-\alpha,1)$ and
$\boldsymbol{\xi}$ uniformly from the unit sphere. To prevent the numerical
instability due to $r\approx 0$, [12] adopts a truncation to ensure the
sampled $r$ is larger than some fixed positive $\epsilon$. The second integral
defined outside of the ball $B_{r_{0}}$ is rewritten as

$\displaystyle\int_{\boldsymbol{y}\notin
B_{r_{0}}(\boldsymbol{x})}\frac{u(\boldsymbol{x})-u(\boldsymbol{y})}{\|\boldsymbol{x}-\boldsymbol{y}\|^{d+\alpha}}d\boldsymbol{y}=\frac{|\mathcal{S}^{d-1}|r_{0}^{-\alpha}}{2\alpha}\mathbb{E}_{\boldsymbol{\xi},r\sim
f_{2}(r)}\left[2u(\boldsymbol{x})-u(\boldsymbol{x}-\boldsymbol{\xi}r)-u(\boldsymbol{x}+\boldsymbol{\xi}r)\right],$
(16)

where we sample $r$ based on $r_{0}/r\sim\operatorname{Beta}(\alpha,1)$ and
$\boldsymbol{\xi}$ uniformly from the unit sphere.

####  3.3.2 Time-fractional Derivative

Monte Carlo on the Beta distribution can also evaluate the time-fractional
derivative, see Guo et al. [12] or (24) for the detailed derivation:

$\displaystyle\frac{\partial^{\gamma}f(t)}{\partial
t^{\gamma}}=\frac{\gamma}{1-\gamma}t^{1-\gamma}\mathbb{E}_{\tau\sim\text{Beta}(1-\gamma,1)}\left[\frac{f(t)-f(t-t\tau)}{t\tau}\right]+\frac{f(t)-f(0)}{t^{\gamma}}.$
(17)

Hence, we can simulate the differential operators in the fractional PDEs using
Monte Carlo simulation. Then, we can plug these PDE operator simulators into
the PINN loss to enable PINN to solve them.

##  4 Proposed Method

This section presents the main proposed methods. Specifically, we first extend
MC-fPINN [12] to tempered fractional cases in subsection 4.1 using a pure
Monte Carlo-based approach, enabling the simulation of tempered fractional
Laplacian to enable PINN’s application to tempered fractional PDEs. Then, in
subsection 4.2 we introduce our improved MC-fPINN that utilizes a more
accurate and efficient quadrature to estimate the 1D integral about $r$ in the
fractional derivatives. Finally, subsection 4.3 extends the improving
techniques to tempered fractional cases.

###  4.1 Monte Carlo Tempered fPINN (MC-tfPINN)

We first extend MC-fPINN [12] to tempered fractional cases using a Monte
Carlo-based approach. The tempered fractional Laplacian can be similarly
defined [5]:

$\displaystyle(-\Delta_{\lambda})^{\alpha/2}u:=C_{d,\alpha,\lambda}\text{P.V.}\int_{\mathbb{R}^{d}}\frac{u(\boldsymbol{x})-u(\boldsymbol{y})}{\exp\left(\lambda\|\boldsymbol{x}-\boldsymbol{y}\|\right)\|\boldsymbol{x}-\boldsymbol{y}\|^{d+\alpha}}d\boldsymbol{y},\quad\alpha\in(0,2),$
(18)

where P.V. denotes the principle value of the integral and
$C_{d,\alpha,\lambda}=\frac{\Gamma\left(\frac{d}{2}\right)}{2\pi^{\frac{d}{2}}\left|\Gamma\left(-{\alpha}\right)\right|}$.
First, we represent the fractional Laplacian as an expectation:

$\displaystyle\int_{\mathbb{R}^{d}}\frac{u(\boldsymbol{x})-u(\boldsymbol{y})}{\exp\left(\lambda\|\boldsymbol{x}-\boldsymbol{y}\|\right)\|\boldsymbol{x}-\boldsymbol{y}\|^{d+\alpha}}d\boldsymbol{y}$
$\displaystyle=\int_{\mathbb{R}^{d}}\frac{u(\boldsymbol{x})-u(\boldsymbol{x}-\boldsymbol{y})}{\exp\left(\lambda\|\boldsymbol{y}\|\right)\|\boldsymbol{y}\|^{d+\alpha}}d\boldsymbol{y}$
$\displaystyle=\frac{1}{2}\int_{\mathbb{R}^{d}}\frac{2u(\boldsymbol{x})-u(\boldsymbol{x}-\boldsymbol{y})-u(\boldsymbol{x}+\boldsymbol{y})}{\exp\left(\lambda\|\boldsymbol{y}\|\right)\|\boldsymbol{y}\|^{d+\alpha}}d\boldsymbol{y}$
$\displaystyle=\frac{1}{2}\int_{\mathcal{S}^{d-1}}\int_{0}^{\infty}\frac{2u(\boldsymbol{x})-u(\boldsymbol{x}-\boldsymbol{\xi}r)-u(\boldsymbol{x}+\boldsymbol{\xi}r)}{\exp\left(\lambda
r\right)r^{1+\alpha}}drd\boldsymbol{\xi}$
$\displaystyle=\frac{1}{2}\int_{\mathcal{S}^{d-1}}\int_{0}^{\infty}\frac{2u(\boldsymbol{x})-u(\boldsymbol{x}-\boldsymbol{\xi}r)-u(\boldsymbol{x}+\boldsymbol{\xi}r)}{r^{2}}\exp\left(-\lambda
r\right)r^{1-\alpha}drd\boldsymbol{\xi}$
$\displaystyle=\frac{1}{2}|\mathcal{S}^{d-1}|\frac{\Gamma(2-\alpha)}{\lambda^{2-\alpha}}\mathbb{E}_{r,\boldsymbol{\xi}}\left[\frac{2u(\boldsymbol{x})-u(\boldsymbol{x}-\boldsymbol{\xi}r)-u(\boldsymbol{x}+\boldsymbol{\xi}r)}{r^{2}}\right],$
(19)

where $r$ obeys the Gamma distribution
$\operatorname{Gamma}(2-\alpha,\lambda)$ and $\boldsymbol{\xi}$ is uniformly
distributed on the $d$-dimensional unit sphere.

Similar to the Monte Carlo fPINN case [12], we use Monte Carlo methods to
simulate the expectation over $r$ and $\boldsymbol{\xi}$. Also, we pick the
sampled $r\geq\epsilon>0$ to prevent the numerical instability due to
$r\approx 0$. When $\lambda>0$, thanks to the tempered term $\exp(-\lambda
r)$, we don’t need to decompose the sampling into two intervals like Monte
Carlo fPINN [12], reducing the hyperparameter number, simplifying the
algorithm and stabilizing the training. Besides, as $r\rightarrow\infty$, the
sampled Gamma distribution exhibits a short tail, reducing the variance.

We have introduced the Monte Carlo simulation for tempered fractional
Laplacian. Meanwhile, the approximation of the tempered time factional
derivative can be reduced to that of the conventional fractional case; see
equation (4). Hence, the methodology will be the same as MC-fPINN for the
tempered time factional derivative.

###  4.2 Improved MC-fPINN

In this subsection, we propose a new class of improved MC-fPINN to accelerate
training and achieve lower error based on quadrature. MC-fPINN suffers from
two issues. Firstly, MC-fPINN relies heavily on numerous hyperparameters,
particularly $r_{0}$ to divide the integral within the fractional Laplacian
into two parts, and $\epsilon$ to mitigate numerical instability, rendering
the model more sensitive to hyperparameter tuning. For instance, with
$\epsilon$, setting it too small leads to numerical instability and non-
convergence, while setting it too large results in bias and error in
estimating the fractional Laplacian.

To this end, we propose an improved MC-fPINN to reduce its hyperparameter
sensitivity. We choose $r_{0}$ as the diameter of the PDE’s support domain
$\Omega$ by default. Then, we can use the Gauss-Jacobi quadrature to compute
the integral concerning $r$ from 0 to $r_{0}$ with a singularity at $r=0$:

$\int_{-1}^{1}(1-x)^{\alpha}(1+x)^{\beta}f(x)dx\approx\sum_{i=0}^{N_{\text{quad}}}w_{i}f(x_{i}).$
(20)

We rescale the integral interval, choose $\alpha,\beta$ accordingly, and
rescale to perform the Gauss-Jacobi quadrature. Then, we compute the 1D
integral on $r$ via quadrature and the $(d-1)$-dimensional integral on the
unit sphere via Monte Carlo.

$\displaystyle\int_{\boldsymbol{y}\in
B_{r_{0}}(\boldsymbol{x})}\frac{u(\boldsymbol{x})-u(\boldsymbol{y})}{\|\boldsymbol{x}-\boldsymbol{y}\|^{d+\alpha}}d\boldsymbol{y}$
$\displaystyle=\frac{1}{2}\int_{\mathcal{S}^{d-1}}\int_{0}^{r_{0}}\frac{2u(\boldsymbol{x})-u(\boldsymbol{x}-\boldsymbol{\xi}r)-u(\boldsymbol{x}+\boldsymbol{\xi}r)}{r^{2}}r^{1-\alpha}drd\boldsymbol{\xi}$
(21)
$\displaystyle\approx\frac{|\mathcal{S}^{d-1}|}{2}\sum_{i=1}^{N}w_{i}\frac{2u(\boldsymbol{x})-u(\boldsymbol{x}-\boldsymbol{\xi}_{i}r_{i})-u(\boldsymbol{x}+\boldsymbol{\xi}_{i}r_{i})}{r_{i}^{2}},$
(22)

where $(w_{i},r_{i})_{i=1}^{N}$ are Gauss-Jacobi quadrature weight and points,
while $\boldsymbol{\xi}_{i}$ are sampled uniformly from the unit sphere. For
the second part of the fractional Laplacian, note that
$u(\boldsymbol{x}\pm\boldsymbol{\xi}r)$ will be 0 since $r\geq r_{0}$ and
$r_{0}$ is the diameter of the support $\Omega$:

$\displaystyle\int_{\boldsymbol{y}\notin
B_{r_{0}}(\boldsymbol{x})}\frac{u(\boldsymbol{x})-u(\boldsymbol{y})}{\|\boldsymbol{x}-\boldsymbol{y}\|^{d+\alpha}}d\boldsymbol{y}=\frac{|\mathcal{S}^{d-1}|r_{0}^{-\alpha}}{\alpha}u(\boldsymbol{x}),$
(23)

The Gauss-Jacobi quadrature for the integral about $r$ has multiple
advantages.

  * •

We do not require two hyperparameters, $r_{0}$ and $\epsilon$. For the
integral about $r$, we replace the Monte Carlo sampling for $r$ with a more
precise quadrature. This is far more accurate than Monte Carlo, especially
when dealing with integrals that have singularities.

  * •

No resampling over $r$ is performed as quadrature points are fixed. In the
regular MC-fPINN, each iteration requires sampling new $r$ values from the
Beta distribution. Thus, our improved method saves time when sampling.

  * •

The second part of the integral in our method is simplified and requires less
neural network inference, effectively reducing the overall computational load
by about 50%.

A similar quadrature can evaluate the time-fractional derivative. We use a
Gauss-Jacobi quadrature for the integral:

$\displaystyle\frac{\partial^{\gamma}f(t)}{\partial t^{\gamma}}$
$\displaystyle=\frac{1}{\Gamma(1-\gamma)}\int_{0}^{t}(t-\tau)^{-\gamma}\frac{\partial
f(\tau)}{\partial\tau}d\tau$
$\displaystyle=\gamma\int^{t}_{0}\tau^{-\gamma}\frac{f(t)-f(t-\tau)}{\tau}d\tau+\frac{f(t)-f(0)}{t^{\gamma}}$
$\displaystyle=\gamma\int^{1}_{0}(t\tau)^{-\gamma}\frac{f(t)-f(t-t\tau)}{t\tau}dt\tau+\frac{f(t)-f(0)}{t^{\gamma}}$
$\displaystyle=\gamma
t^{1-\gamma}\int^{1}_{0}(\tau)^{-\gamma}\frac{f(t)-f(t-t\tau)}{t\tau}d\tau+\frac{f(t)-f(0)}{t^{\gamma}}$
(24) $\displaystyle\approx\gamma
t^{1-\gamma}\left(\sum_{i=1}^{N}w_{i}\frac{f(t)-f(t-t\tau_{i})}{t\tau_{i}}\right)+\frac{f(t)-f(0)}{t^{\gamma}},$
(25)

where $(w_{i},\tau_{i})_{i=1}^{N}$ are Gauss-Jacobi quadrature weight and
points.

###  4.3 Improved tfPINN

For tfPINN, we estimate the integral on $r$ with generalized Gauss–Laguerre
quadrature, which is tailored for the following type of integral with a
singularity at $x=0$ for $\alpha>-1$ and $\lambda>0$:

$\int_{0}^{\infty}x^{\alpha}e^{-\lambda
x}f(x)dx\approx\sum_{i=0}^{N_{\text{quad}}}w_{i}f(x_{i}).$ (26)

We compute the 1D integral on $r$ via quadrature and the integral on the unit
sphere via Monte Carlo. By the equivalent form (4.1), we write

$\displaystyle\int_{\mathbb{R}^{d}}\frac{u(\boldsymbol{x})-u(\boldsymbol{y})}{\exp\left(\lambda\|\boldsymbol{x}-\boldsymbol{y}\|\right)\|\boldsymbol{x}-\boldsymbol{y}\|^{d+\alpha}}d\boldsymbol{y}$
$\displaystyle=\frac{1}{2}|\mathcal{S}^{d-1}|\mathbb{E}_{\boldsymbol{\xi}}\left[\int_{0}^{\infty}\frac{2u(\boldsymbol{x})-u(\boldsymbol{x}-\boldsymbol{\xi}r)-u(\boldsymbol{x}+\boldsymbol{\xi}r)}{r^{2}}\exp\left(-\lambda
r\right)r^{1-\alpha}dr\right]$ (27)
$\displaystyle\approx\frac{1}{2}|\mathcal{S}^{d-1}|\sum_{i=0}^{N}w_{i}\frac{2u(\boldsymbol{x})-u(\boldsymbol{x}-\boldsymbol{\xi}_{i}r_{i})-u(\boldsymbol{x}+\boldsymbol{\xi}_{i}r_{i})}{r_{i}^{2}}.$
(28)

where $(w_{i},r_{i})_{i=1}^{N}$ are Gauss-Laguerre quadrature weight and
points, while $\boldsymbol{\xi}_{i}$ are sampled uniformly from the unit
sphere.

###  4.4 Discussion: On Quasi-Monte Carlo

We have employed Gauss-quadrature rules to compute the 1D integral about $r$
more accurately and efficiently. For $\boldsymbol{\xi}$, we still utilize
standard Monte Carlo because $\boldsymbol{\xi}\in\mathbb{R}^{d}$ and
quadrature will suffer from the CoD. To further reduce the variance of Monte
Carlo over $\boldsymbol{\xi}\in\mathbb{R}^{d}$, Quasi-Monte Carlo (QMC) is an
appealing option. In Appendix A, we also compare the performance of using QMC
against using Monte Carlo methods. We do not observe any improvement in
accuracy using QMC.

##  5 Computational Experiments

###  5.1 Fractional Poisson Equation

Consider the fractional Poisson equation $(-\Delta)^{\alpha/2}u=f$ in the unit
ball. We will test the following exact solution:

$u_{\text{exact}}(\boldsymbol{x})=\left(1-\|\boldsymbol{x}\|^{2}\right)^{\alpha/2}\left(c_{1,0}+\sum_{i=1}^{d}c_{1,i}\boldsymbol{x}_{i}\right)+\left(1-\|\boldsymbol{x}\|^{2}\right)^{1+\alpha/2}\left(c_{2,0}+\sum_{i=1}^{d}c_{2,i}\boldsymbol{x}_{i}\right)$
(29)

where $c_{1,i}$ and $c_{2,i}$ are independently sampled from a unit Gaussian.
It is a linear combination of the functions in Table 1 to form a relatively
complicated and anisotropic exact solution whose forcing term can be computed
analytically.

$u(x)$ in the unit ball |  $(-\Delta)^{\alpha/2}u(x)$ in the unit ball
---|---
$\left(1-\|\boldsymbol{x}\|^{2}\right)^{\alpha/2}$ | $-2^{\alpha}\Gamma\left(\frac{\alpha}{2}+1\right)\Gamma\left(\frac{\alpha+d}{2}\right)\Gamma\left(\frac{d}{2}\right)^{-1}$
$\left(1-\|\boldsymbol{x}\|^{2}\right)^{1+\alpha/2}$ | $-2^{\alpha}\Gamma\left(\frac{\alpha}{2}+2\right)\Gamma\left(\frac{\alpha+d}{2}\right)\Gamma\left(\frac{d}{2}\right)^{-1}\left(1-\left(1+\frac{\alpha}{d}\right)\|\boldsymbol{x}\|^{2}\right)$
$\left(1-\|\boldsymbol{x}\|^{2}\right)^{\alpha/2}\boldsymbol{x}_{d}$ | $-2^{\alpha}\Gamma\left(\frac{\alpha}{2}+1\right)\Gamma\left(\frac{\alpha+d}{2}+1\right)\Gamma\left(\frac{d}{2}+1\right)^{-1}\boldsymbol{x}_{d}$
$\left(1-\|\boldsymbol{x}\|^{2}\right)^{1+\alpha/2}\boldsymbol{x}_{d}$ | $-2^{\alpha}\Gamma\left(\frac{\alpha}{2}+2\right)\Gamma\left(\frac{\alpha+d}{2}+1\right)\Gamma\left(\frac{d}{2}+1\right)^{-1}\left(1-\left(1+\frac{\alpha}{d+2}\right)\|\boldsymbol{x}\|^{2}\right)\boldsymbol{x}_{d}$
Table 1: Fractional Laplacian of some functions, from [7].

In MC-fPINN [12], a simple isotropic exact solution
$\left(1-\|\boldsymbol{x}\|^{2}\right)^{1+\alpha/2}$ is tested. We extend
their experiment to complicated anisotropic exact solutions and show that our
improved MC-fPINN outperforms.

For a fair comparison, we use the hyperparameters from MC-fPINN [12]. We
choose $\alpha=1.5$ for a test. We use $\epsilon=1e-6$ for preventing
numerical issues and truncating the sampled $r$ from Beta distribution as in
MC-fPINN. Here $r_{0}$ is set to 0.25 in MC-fPINN, according to the
hyperparameter study in Guo et al. [12].

We opt for 100 residual points and 64 points for the Monte Carlo/quadrature
batch size. We randomly sampled 20K test points within the unit ball. We
optimize both models for 1 million epochs with Adam [21], whose initial
learning rate is 1e-3 and linearly decays to zero. We adopt the following
model structure to satisfy the zero boundary condition with hard constraint
and to avoid the boundary loss [25]
$\text{ReLU}(1-\|\boldsymbol{x}\|_{2}^{2})u_{\theta}(\boldsymbol{x}),$ where
$u_{\theta}(\boldsymbol{x})$ is the neural neural network model. The neural
network model is a four-layer, fully-connected network with 128 hidden units
and a hyperbolic tangent activation function. Since we hard-constrained the
boundary condition, the only loss function will be the residual loss. Hence,
there is no weighting hyperparameter.

Method | Dimension | 100 | 1,000 | 10,000 | 100,000
---|---|---|---|---|---
MC-fPINN (Guo et al. [12]) | Speed | 261it/s | 223it/s | 77it/s | 9it/s
Rel. $L_{2}$ | 2.86E-02 | 4.99E-02 | 5.85E-02 | 8.02E-02
Improved MC-fPINN (Ours) | Speed | 1092it/s | 747it/s | 174it/s | 20it/s
Rel. $L_{2}$ | 2.84E-02 | 3.68E-02 | 4.94E-02 | 7.10E-02
Table 2: Results for the fractional Poisson equation. Our improved MC-fPINN
runs faster and performs better with slightly lower errors than the original
MC-fPINN by Guo et al. [12]. Faster speed and lower error are bold.

The results are presented in Table 2 with different methods’ speed (iteration
per second) and relative $L_{2}$ error. Our improved MC-fPINN consistently
outperforms MC-fPINN by Guo et al. [12] regarding speed and relative $L_{2}$
error performance. This is attributed to the Gauss quadrature for the integral
about $r$ in our method, rendering it more precise compared to Monte Carlo.
Additionally, there is no need to resample the integration points for $r$ at
each epoch, enhancing effectiveness and efficiency compared with vanilla MC-
fPINN [12]. Moreover, our improved MC-fPINN eliminates the need for additional
hyperparameters like $r_{0}$ and $\epsilon$, making it more robust in
hyperparameter selection during practical applications. Furthermore, in our
method, the second part of the integral in estimating the fractional Laplacian
is simplified and it requires less neural network inference, effectively
reducing the overall computational load by about 50%, since we choose $r_{0}$
as the diameter of the PDE support, which is $r_{0}=2$ in this test case.

###  5.2 Tempered Fractional Poisson Equation

We consider high-dimensional tempered fractional Poisson equations
$\left(-\Delta_{\lambda}\right)^{\alpha/2}u=f$ within the unit ball, zero
boundary conditions, and various anisotropic and inseparable solutions. The
following nontrivial exact solutions are considered in related work [16, 17,
15] to evaluate the high-dimensional PINN, specifically one with two-body
interaction and another with three-body interaction:

$u_{\text{two-
body}}(\boldsymbol{x})=\text{ReLU}\left(1-\|\boldsymbol{x}\|_{2}^{2}\right)\left(\sum_{i=1}^{d-1}c_{i}\sin(\boldsymbol{x}_{i}+\cos(\boldsymbol{x}_{i+1})+\boldsymbol{x}_{i+1}\cos(\boldsymbol{x}_{i}))\right),$
(30) $u_{\text{three-
body}}(\boldsymbol{x})=\text{ReLU}\left(1-\|\boldsymbol{x}\|_{2}^{2}\right)\left(\sum_{i=1}^{d-2}c_{i}\exp(\boldsymbol{x}_{i}\boldsymbol{x}_{i+1}\boldsymbol{x}_{i+2})\right),$
(31)

where $c_{i}\sim\mathcal{N}(0,1)$ and there are two-body interactions between
each $(\boldsymbol{x}_{i},\boldsymbol{x}_{i+1})$, three-body interactions
between each pair
$(\boldsymbol{x}_{i},\boldsymbol{x}_{i+1},\boldsymbol{x}_{i+2})$,
respectively. The forcing terms for training residual points are computed with
a Monte Carlo simulation with 1024 samples.

We choose $\alpha=0.5/1.5$ and $\lambda=1$ for the tempered fractional
Laplacian. We set $\epsilon=1e-6$ for MC-tfPINN to truncate $r$ to prevent
numerical instability. The Monte Carlo/quadrature sample size for the neural
network is 64 during training.

The model is a 4-layer fully connected network with 128 hidden units activated
by $\tanh$, which is trained via Adam [21] for 10K epochs, with an initial
learning rate 1e-3, which linearly decays to zero at the end of the
optimization. We select 100 random residual points at each Adam epoch and 20K
fixed testing points uniformly from the unit ball. We adopt the following
model structure to satisfy the zero boundary condition with hard constraint
and to avoid the boundary loss
[25]$\text{ReLU}(1-\|\boldsymbol{x}\|_{2}^{2})u_{\theta}(\boldsymbol{x}),$
where $u_{\theta}(\boldsymbol{x})$ is the neural network model. Since we hard-
constrained the boundary condition, the only loss function will be the
residual loss. Hence, there is no weighting hyperparameter.

| Metric/Dimension | $10^{1}$ | $10^{2}$ | $10^{3}$ | $10^{4}$ | $10^{5}$
---|---|---|---|---|---|---
Vanilla MC-tfPINN | Speed | 864it/s | 342it/s | 108it/s | 18it/s | 2it/s
Error_TwoBody $\alpha=0.5$ | 1.86E-3 | 7.27E-3 | 7.41E-3 | 3.56E-3 | 2.17E-3
Error_TwoBody $\alpha=1.5$ | 1.79E-3 | 7.85E-3 | 8.00E-3 | 4.57E-3 | 2.36E-3
Error_ThreeBody $\alpha=0.5$ | 4.01E-3 | 4.94E-4 | 4.32E-4 | 4.54E-4 | 5.69E-4
Error_ThreeBody $\alpha=1.5$ | 3.97E-3 | 5.02E-4 | 4.44E-4 | 5.02E-4 | 5.58E-4
Improved MC-tfPINN | Speed | 1036it/s | 468it/s | 117it/s | 21it/s | 3it/s
Error_TwoBody $\alpha=0.5$ | 2.03E-3 | 8.92E-3 | 6.57E-3 | 3.64E-3 | 2.09E-3
Error_TwoBody $\alpha=1.5$ | 1.97E-3 | 8.08E-3 | 6.42E-3 | 4.08E-3 | 2.11E-3
Error_ThreeBody $\alpha=0.5$ | 4.02E-3 | 5.62E-4 | 4.00E-4 | 4.81E-4 | 5.09E-4
Error_ThreeBody $\alpha=1.5$ | 3.81E-3 | 5.59E-4 | 4.25E-4 | 4.97E-4 | 5.28E-4
Table 3: Results for Monte Carlo tempered fractional PINN on various
inseparable exact solutions with a tempered fractional Poisson equation. The
error is stable, around 1e-3 to 1e-4, and the cost is sublinear. The two
models perform similarly so we do not bold the better result.

The results are shown in Table 3 with different methods’ speed (iteration per
second) and relative $L_{2}$ error. This test case is relatively simple. Thus,
both variants of tempered fractional PINN yield low errors ranging from 1e-3
to 1e-4. Additionally, we validated that the improved version exhibits faster
convergence. In this scenario, the forcing term of the tempered fractional PDE
is not analytical, hence we resort to Monte Carlo to unbiasedly estimate the
forcing term. In other words, the forcing term used during the optimization of
MC-tfPINN is noisy. Despite the presence of noise in the forcing term, both
variants of tfPINN still yield low errors, demonstrating the robustness of our
approach.

###  5.3 Time-Dependent Tempered Fractional Diffusion Equation

We demonstrate how our MC-tfPINN can deal with a time-dependent tempered
fractional diffusion equation:

$\displaystyle\frac{\partial^{\gamma}u(\boldsymbol{x},t)}{\partial
t^{\gamma}}+c\left(-\Delta_{\lambda}\right)^{\alpha/2}u(\boldsymbol{x},t)+\boldsymbol{v}\cdot\nabla_{\boldsymbol{x}}u(\boldsymbol{x},t)$
$\displaystyle=f(\boldsymbol{x},t),\quad(\boldsymbol{x},t)\in\Omega\times(0,T].$
(32) $\displaystyle u(\boldsymbol{x},t=0)$
$\displaystyle=0,\quad\boldsymbol{x}\in\Omega.$ (33) $\displaystyle
u(\boldsymbol{x},t)$
$\displaystyle=0,\quad(\boldsymbol{x},t)\in\Omega^{c}\times(0,T].$ (34)

The fabricated solution is highly coupled:

$u_{\text{two-
body}}(\boldsymbol{x},t)=\text{ReLU}\left(1-\|\boldsymbol{x}\|_{2}^{2}\right)\left(\sum_{i=1}^{d-1}c_{i}\sin\left[t\cdot(\boldsymbol{x}_{i}+\cos(\boldsymbol{x}_{i+1})+\boldsymbol{x}_{i+1}\cos(\boldsymbol{x}_{i}))\right]\right),$
(35)

where $c_{i}$ is randomly drawn from a unit Gaussian. The zero boundary and
initial conditions pose additional challenges in residual fitting on the given
complicated forcing term $f(\boldsymbol{x},t)$ since they will not leak any
information about the exact solution. The terminal time is set to $T=1$, the
time-fractional order is set to $\gamma=0.5$, and the velocity
$\boldsymbol{v}\in\mathbb{R}^{d}$ is uniformly drawn from $[0,1]$. For other
coefficients in the PDE operator, we choose the following combinations: (1)
$\alpha=0.5,\lambda=1$; (2) $\alpha=1.5,\lambda=1$; (3)
$\alpha=0.5,\lambda=2$; (4) $\alpha=0.5,\lambda=0.5$. We want to test if our
MC-tfPINN and its improved version can deal with diverse settings with
different PDE operator coefficients.

Here are the implementation details. The model is a 4-layer fully connected
network with 128 hidden units activated by $\tanh$, which is trained via Adam
[21] for 100K epochs, with an initial learning rate 1e-3, which linearly
decays to zero at the end of the optimization. We select 100 random residual
points at each Adam epoch and 20K fixed testing points uniformly from the unit
ball. We adopt the following model structure to satisfy the zero boundary
condition with hard constraint and to avoid the boundary and initial losses
[25] $t\text{ReLU}(1-\|\boldsymbol{x}\|_{2}^{2})u_{\theta}(\boldsymbol{x},t),$
where $u_{\theta}(\boldsymbol{x},t)$ is the neural network model. Since we
hard-constrained the boundary and initial conditions, the only loss function
will be the residual loss. Hence, there is no weighting hyperparameter.

Method | $\alpha$ | $\lambda$ | Dimension | 10 | 100 | 1,000 | 10,000
---|---|---|---|---|---|---|---
Vanilla MC-tfPINN | 0.5 | 1 | Speed | 182it/s | 161it/s | 20it/s | 3it/s
Rel. $L_{2}$ | 2.49E-02 | 4.94E-02 | 6.91E-02 | 7.13E-02
Improved MC-tfPINN | 0.5 | 1 | Speed | 202it/s | 189it/s | 21it/s | 4it/s
Rel. $L_{2}$ | 1.63E-02 | 2.69E-02 | 3.73E-02 | 4.43E-02
Method | $\alpha$ | $\lambda$ | Dimension | 10 | 100 | 1,000 | 10,000
Vanilla MC-tfPINN | 1.5 | 1 | Speed | 182it/s | 161it/s | 20it/s | 3it/s
Rel. $L_{2}$ | 2.97E-02 | 3.33E-02 | 4.79E-02 | 4.94E-02
Improved MC-tfPINN | 1.5 | 1 | Speed | 202it/s | 189it/s | 21it/s | 4it/s
Rel. $L_{2}$ | 2.10E-02 | 2.26E-02 | 4.79E-02 | 4.61E-02
Method | $\alpha$ | $\lambda$ | Dimension | 10 | 100 | 1,000 | 10,000
Vanilla MC-tfPINN | 0.5 | 2 | Speed | 182it/s | 161it/s | 20it/s | 3it/s
Rel. $L_{2}$ | 4.10E-03 | 3.67E-02 | 5.39E-02 | 5.53E-02
Improved MC-tfPINN | 0.5 | 2 | Speed | 202it/s | 189it/s | 21it/s | 4it/s
Rel. $L_{2}$ | 4.04E-03 | 4.89E-03 | 4.69E-02 | 4.59E-02
Method | $\alpha$ | $\lambda$ | Dimension | 10 | 100 | 1,000 | 10,000
Vanilla MC-tfPINN | 0.5 | 0.5 | Speed | 182it/s | 161it/s | 20it/s | 3it/s
Rel. $L_{2}$ | 9.31E-02 | 8.98E-02 | 6.91E-02 | 7.70E-02
Improved MC-tfPINN | 0.5 | 0.5 | Speed | 202it/s | 189it/s | 21it/s | 4it/s
Rel. $L_{2}$ | 1.61E-03 | 1.07E-02 | 4.56E-02 | 4.59E-02
Table 4: Results of MC-tfPINN on a time-dependent tempered fractional
diffusion equation. Improved MC-tfPINN is more efficient and has lower error
thanks to its accurate quadrature. Faster speed and lower error are bold.

The results are presented in Table 4. Firstly, this test case is notably more
challenging compared to the previous one, primarily due to the time-fractional
derivative. Additionally, the exact solution in this case is highly coupled in
both time and spatial $\boldsymbol{x}$, rendering it exceedingly complex and
inseparable. However, our two variants of tfPINN still manage to yield around
5% error across all dimensions, with speed growing sublinearly as the
dimension increases and the error increasing slowly with the dimension. The
speeds are the same if the PDE dimensionality and the model are the same since
the PDE operator coefficient does not affect the speed. Furthermore, we
reaffirm that improved MC-tfPINN with quadrature always trains faster and
achieves smaller errors. In conclusion, we extensively tested MC-tfPINN on
tempered fractional PDEs with different PDE operator parameters, demonstrating
that tfPINN is fairly stable with respect to the parameters alpha (fractional
order) and lambda (tempering factor) of the PDE operator. This showcases the
stability and scalability of our proposed MC-tfPINN across different
dimensions and various PDE operator parameter settings.

###  5.4 Inverse Tempered Fractional PDEs

We will test the two-body interaction exact solution in equation (30) with the
tempered fractional Laplacian equation
$\left(-\Delta_{\lambda}\right)^{\alpha/2}u=f$ in the unit ball. Unlike the
forward problem, we know the values of the exact solution on several points
and wish to infer the fractional order $\alpha$ or the tempering factor
$\lambda$. The implementation details on optimizing over $\alpha$ or $\lambda$
are presented in Appendix B. Specifically, we consider the following two
inverse problems.

Problem 1: Inverse problem of $\alpha$. We fix $\lambda=1$ and initial the
unknown fractional order as $\alpha=1.5$. The exact fractional order $\alpha$
is either 1.4 or 1.6.

Problem 2: Inverse problem of $\lambda$. We fix $\alpha=1.5$ and initial the
unknown tempering factor as $\lambda=1$. The exact $\lambda$ is set to 2.

The forcing terms for training residual points are computed with a Monte Carlo
simulation with 1024 samples, while the Monte Carlo/quadrature sample size for
the neural network is 128 during training. $\epsilon=1e-6$ for MC-tfPINN to
truncate $r$ to prevent numerical instability.

Here are more implementation details. The model is a 4-layer fully connected
network with 128 hidden units activated by $\tanh$, which is trained via Adam
[21] for 100K epochs, with an initial learning rate 1e-3, which linearly
decays to zero at the end of the optimization. We select 100 random residual
points and data points at each Adam epoch. We adopt the following model
structure to satisfy the zero boundary condition with hard constraint and to
avoid the boundary loss [25]
$\text{ReLU}(1-\|\boldsymbol{x}\|_{2}^{2})u_{\theta}(\boldsymbol{x}),$ where
$u_{\theta}(\boldsymbol{x})$ is a fully connected neural network. Since we
hard-constrained the boundary condition, the loss functions will be the
residual loss and the data loss in the inverse problem and the weights are
one.

Method | Problem | Dimension | Exact $\alpha$ | Identified $\alpha$ | Rel $L_{1}$ | Time
---|---|---|---|---|---|---
MC-tfPINN | 1 | 100 | 0.4 | 0.3902741 | 2.43% | 6min
Improved MC-tfPINN | 1 | 100 | 0.4 | 0.398084 | 0.48% | 5min
MC-tfPINN | 1 | 1000 | 0.4 | 0.3905588 | 2.36% | 17min
Improved MC-tfPINN | 1 | 1000 | 0.4 | 0.3989002 |  0.27% | 16min
MC-tfPINN | 1 | 10000 | 0.4 | 0.3907224 | 2.32% | 103min
Improved MC-tfPINN | 1 | 10000 | 0.4 | 0.3988998 | 0.28% | 102min
MC-tfPINN | 1 | 100 | 0.6 | 0.5896347 | 1.73% | 6min
Improved MC-tfPINN | 1 | 100 | 0.6 | 0.5983624 | 0.27% | 5min
MC-tfPINN | 1 | 1000 | 0.6 | 0.5901451 | 1.64% | 17min
Improved MC-tfPINN | 1 | 1000 | 0.6 | 0.5992336 | 0.13% | 16min
MC-tfPINN | 1 | 10000 | 0.6 | 0.5902596 | 1.62% | 103min
MC-tfPINN | 1 | 10000 | 0.6 | 0.5991701 | 0.14% | 102min
Method | Problem | Dimension | Exact $\lambda$ | Identified $\lambda$ | Rel $L_{1}$ | Time
MC-tfPINN | 2 | 100 | 2 | 2.005856 | 0.29% | 5min
Improved MC-fPINN | 2 | 100 | 2 | 2.00189 | 0.09% | 4min
MC-tfPINN | 2 | 1000 | 2 | 2.004311 | 0.22% | 18min
Improved MC-tfPINN | 2 | 1000 | 2 | 2.000428 | 0.02% | 17min
MC-tfPINN | 2 | 10000 | 2 | 2.00E+00 | 0.20% | 116min
Improved MC-tfPINN | 2 | 10000 | 2 | 2.000902 | 0.05% | 114min
Table 5: Results for Monte Carlo tempered fractional PINN on inverse problems,
where we report the problem setting and the relative $L_{1}$ error for the
identified PDE operator coefficient. We also present the total running time
for the algorithms. We bold the best results.

Table 5 presents the identified PDE operator coefficient value and its
corresponding relative $L_{1}$ error given the exact coefficient value. MC-
tfPINN can satisfactorily solve both forward and inverse PDE problems,
demonstrating its versatility. In addition, we further prove that the improved
MC-tfPINN using quadrature can be trained faster and yield better results
thanks to its improved accuracy and efficiency.

##  6 Summary

We extend the application of physics-informed neural networks (PINNs) to high-
dimensional fractional and tempered fractional partial differential equations
(PDEs) by developing the Monte Carlo tempered fractional PINN (MC-tfPINN)
based on the previous work Monte Carlo fPINN [12] (MC-fPINN). Furthermore,
this research enhances MC-fPINN and MC-tfPINN with Gaussian quadrature to
replace traditional Monte Carlo methods when computing the 1D radial singular
integral within their PDE operators. This new approach significantly reduces
computational errors and increases convergence speed by addressing the
limitations associated with high variance and hyperparameter sensitivity in
previous models MC-fPINN [12], which rely on a pure Monte Carlo simulation.
The proposed methods are validated across a range of forward and inverse
problems of (tempered) fractional PDEs, demonstrating their superior accuracy
and efficiency in handling up to 100,000 dimensions, a marked improvement over
existing methods. We also show that the improved methods based on quadrature
can consistently outperform the vanilla models based on Monte Carlo solely in
terms of speed and accuracy.

##  Appendix A Unit Test on Quasi-Monte Carlo (QMC)

We conduct a unit test on the accuracy of different fPINN simulations in
estimating the fractional Laplacian. Specifically, we consider the following
three models.

  1. 1.

The vanilla MC-fPINN [12]. It uses Monte Carlo for the integral of both $r$
and $\boldsymbol{\xi}$.

  2. 2.

Our improved MC-fPINN. It uses Gauss-Jacobi quadrature for the 1D integral of
$r$ and conventional Monte Carlo for the integral of $\boldsymbol{\xi}$.

  3. 3.

QMC-based improved fPINN. It uses Gauss-Jacobi quadrature for the 1D integral
of $r$ and Quasi-Monte Carlo for the integral of $\boldsymbol{\xi}$.

We will test the following exact solution with the fractional Laplacian
equation $(-\Delta)^{\alpha/2}u=f$ in the unit ball:

$u_{\text{exact}}(\boldsymbol{x})=\left(1-\|\boldsymbol{x}\|^{2}\right)^{\alpha/2}\left(c_{1,0}+\sum_{i=1}^{d}c_{1,i}\boldsymbol{x}_{i}\right)+\left(1-\|\boldsymbol{x}\|^{2}\right)^{1+\alpha/2}\left(c_{2,0}+\sum_{i=1}^{d}c_{2,i}\boldsymbol{x}_{i}\right),$
(36)

where $c_{1,i}$ and $c_{2,i}$ are independently sampled from a unit Gaussian.
The exact forcing term $f$ can be computed analytically based on Dyda [7]. We
conducted the 100D and 1000D test cases and repeated the experiment ten times
using ten different random seeds. We report the relative $L_{2}$ error for
fractional Laplacian forcing term estimation in Table 6. We randomly sampled
20K test points within the unit ball to test different methods’ accuracy in
estimating the forcing term. We choose 64 points for the Monte
Carlo/quadrature batch size. $\epsilon$ for preventing numerical issues and
truncating the sampled $r$ from Beta distribution in MC-fPINN is chosen as
1e-6. $\alpha$ is set to 1.5. First, our improved MC-fPINN outperforms MC-
fPINN [12] in lower error and variance thanks to the accurate quadrature.
Besides, improved MC-fPINN and QMC-fPINN perform similarly, i.e., the QMC in
the latter does not necessarily reduce variance and error in high dimensions.
These results validate the effectiveness of the key quadrature component in
our improved MC-fPINN. We observe from the table that QMC in the integral of
$\boldsymbol{\xi}$ does not improve the accuracy.

Method/Dimension | 100D | 1,000D
---|---|---
MC-fPINN | 1.32E-2$\pm$1.23E-2 | 4.05E-3$\pm$1.035E-3
Improved MC-fPINN | 1.03E-2$\pm$6.89E-3 | 4.02E-5$\pm$3.83E-5
QMC-fPINN | 1.24E-2$\pm$7.02E-3 | 6.21E-5$\pm$1.66E-5
Table 6: Results for the unit test on QMC.

##  Appendix B Implementation Detail on Inverse Problem

Application to inverse problems requires careful handling of the fractional
order. Unlike the forward problem, we need to optimize the PDE coefficient in
the inverse problem.

###  B.1 Time-fractional Derivative

By the representation of the Caputo time-fractional derivative (24) and (17),
we have

$\displaystyle\frac{\partial^{\gamma}f(t)}{\partial t^{\gamma}}$
$\displaystyle=\gamma
t^{1-\gamma}\int^{1}_{0}(\tau)^{-\gamma}\frac{f(t)-f(t-t\tau)}{t\tau}d\tau+\frac{f(t)-f(0)}{t^{\gamma}}$
$\displaystyle=\frac{\gamma}{1-\gamma}t^{1-\gamma}\mathbb{E}_{\tau\sim\text{Beta}(1-\gamma,1)}\left[\frac{f(t)-f(t-t\tau)}{t\tau}\right]+\frac{f(t)-f(0)}{t^{\gamma}}.$

We want to infer $\gamma$ while the Monte Carlo simulation is based on a
distribution with $\gamma$ as a parameter. Thus, resampling is required once
we have an update on $\gamma$. Hence, we assume the unknown order
$\gamma\in(\gamma_{L},\gamma_{H})$, then

$\displaystyle\frac{\partial^{\gamma}f(t)}{\partial t^{\gamma}}$
$\displaystyle=\gamma
t^{1-\gamma}\int^{1}_{0}(\tau)^{-\gamma}\frac{f(t)-f(t-t\tau)}{t\tau}d\tau+\frac{f(t)-f(0)}{t^{\gamma}}$
$\displaystyle=\gamma
t^{1-\gamma}\int^{1}_{0}(\tau)^{-\gamma_{H}}(\tau)^{\gamma_{H}-\gamma}\frac{f(t)-f(t-t\tau)}{t\tau}d\tau+\frac{f(t)-f(0)}{t^{\gamma}}$
$\displaystyle=\frac{\gamma}{1-\gamma_{H}}t^{1-\gamma}\mathbb{E}_{\tau\sim\text{Beta}(1-\gamma_{H},1)}\left[(\tau)^{\gamma_{H}-\gamma}\frac{f(t)-f(t-t\tau)}{t\tau}\right]+\frac{f(t)-f(0)}{t^{\gamma}}.$

Then we may optimize over $\gamma$ without resampling since $\gamma_{H}$ is
fixed. We parameterize the trainable $\gamma$ as follows:

$\displaystyle\gamma=\operatorname{Sigmoid}(\Gamma)\cdot(\gamma_{H}-\gamma_{L})+\gamma_{L},$
(37)

where $\operatorname{Sigmoid}$ controls its output value in $(0,1)$ and
$\Gamma$ is a unconstraint trainable parameter. We initial $\Gamma=0$ to
initialize $\gamma=(\gamma_{H}+\gamma_{L})/2$ as
$\operatorname{Sigmoid}(0)=0.5$.

###  B.2 Tempered Fractional Laplacian

For the tempered fractional Laplacian, we don’t need to assume prior knowledge
on $\lambda>0$ but require assuming $\alpha\in(\alpha_{L},\alpha_{H})$. By the
representation (4.1), we have

$\displaystyle\quad\int_{\mathbb{R}^{d}}\frac{u(\boldsymbol{x})-u(\boldsymbol{y})}{\exp\left(\lambda\|\boldsymbol{x}-\boldsymbol{y}\|\right)\|\boldsymbol{x}-\boldsymbol{y}\|^{d+\alpha}}d\boldsymbol{y}$
$\displaystyle=\frac{1}{2}\int_{\mathcal{S}^{d-1}}\int_{0}^{\infty}\frac{2u(\boldsymbol{x})-u(\boldsymbol{x}-\boldsymbol{\xi}r)-u(\boldsymbol{x}+\boldsymbol{\xi}r)}{r^{2}}\exp\left(-\lambda
r\right)r^{1-\alpha}drd\boldsymbol{\xi}$
$\displaystyle=\frac{1}{2}\lambda^{\alpha}\int_{\mathcal{S}^{d-1}}\int_{0}^{\infty}\frac{2u(\boldsymbol{x})-u(\boldsymbol{x}-\boldsymbol{\xi}r/\lambda)-u(\boldsymbol{x}+\boldsymbol{\xi}r/\lambda)}{r^{2}}\exp\left(-r\right)r^{1-\alpha}drd\boldsymbol{\xi}$
$\displaystyle=\frac{1}{2}\lambda^{\alpha}\int_{\mathcal{S}^{d-1}}\int_{0}^{\infty}\frac{2u(\boldsymbol{x})-u(\boldsymbol{x}-\boldsymbol{\xi}r/\lambda)-u(\boldsymbol{x}+\boldsymbol{\xi}r/\lambda)}{r^{2}}\exp\left(-r\right)r^{1-\alpha_{H}}r^{\alpha_{H}-\alpha}drd\boldsymbol{\xi}$
$\displaystyle=\frac{1}{2}\lambda^{\alpha}|\mathcal{S}^{d-1}|\Gamma(2-\alpha_{H})\mathbb{E}_{r,\boldsymbol{\xi}}\left[\frac{2u(\boldsymbol{x})-u(\boldsymbol{x}-\boldsymbol{\xi}r/\lambda)-u(\boldsymbol{x}+\boldsymbol{\xi}r/\lambda)}{r^{2}}r^{\alpha_{H}-\alpha}\right],$

where $r$ obeys the Gamma distribution $\operatorname{Gamma}(2-\alpha_{H},1)$
and $\boldsymbol{\xi}$ is uniformly distributed the $d$-dimensional unit
sphere. Hence, we can optimize over $\alpha$ or $\lambda$ without resampling.

## References

  * [1] Boris Baeumer and Mark M. Meerschaert.  Tempered stable lévy motion and transient super-diffusion.  Journal of Computational and Applied Mathematics, 233(10):2438–2448, 2010.
  * [2] Andrew R Barron.  Universal approximation bounds for superpositions of a sigmoidal function.  IEEE Transactions on Information theory, 39(3):930–945, 1993.
  * [3] Xiaoli Chen, Liu Yang, Jinqiao Duan, and George Em Karniadakis.  Solving inverse stochastic problems from discrete particle observations using the Fokker–Planck equation and physics-informed neural networks.  SIAM Journal on Scientific Computing (SISC), 43(3):B811–B830, 2021.
  * [4] Tim De Ryck, Ameya D Jagtap, and Siddhartha Mishra.  Error estimates for physics informed neural networks approximating the navier-stokes equations.  arXiv preprint arXiv:2203.09346, 2022.
  * [5] Weihua Deng, Buyang Li, Wenyi Tian, and Pingwen Zhang.  Boundary problems for the fractional and tempered fractional operators.  Multiscale Modeling & Simulation, 16(1):125–149, 2018.
  * [6] Siwei Duo and Yanzhi Zhang.  Numerical approximations for the tempered fractional Laplacian: Error analysis and applications.  Journal of Scientific Computing, 81:569–593, 2019.
  * [7] Bartłlomiej Dyda.  Fractional calculus for power functions and eigenvalues of the fractional Laplacian.  Fractional Calculus and Applied Analysis, 15(4):536–555, 2012.
  * [8] Marta D’Elia, Qiang Du, Christian Glusa, Max Gunzburger, Xiaochuan Tian, and Zhi Zhou.  Numerical methods for nonlocal and fractional models.  Acta Numerica, 29:1–124, 2020.
  * [9] Hafiz Muhammad Fahad, Arran Fernandez, Mujeeb Ur Rehman, and Maham Siddiqi.  Tempered and hadamard-type fractional calculus with respect to functions.  Mediterranean Journal of Mathematics, 18(4):143, 2021.
  * [10] Ali Nosrati Firoozsalari, Hassan Dana Mazraeh, Alireza Afzal Aghaei, and Kourosh Parand.  deepfdenet: A novel neural network architecture for solving fractional differential equations.  arXiv preprint arXiv:2309.07684, 2023.
  * [11] Ting Gao, Jinqiao Duan, and Xiaofan Li.  Fokker–Planck equations for stochastic dynamical systems with symmetric lévy motions.  Applied Mathematics and Computation, 278:1–20, 2016.
  * [12] Ling Guo, Hao Wu, Xiaochen Yu, and Tao Zhou.  Monte Carlo fPINNs: Deep learning method for forward and inverse problems involving high dimensional fractional partial differential equations.  Computer Methods in Applied Mechanics and Engineering, 400:115523, 2022.
  * [13] Emmanuel Hanert and Cécile Piret.  A Chebyshev pseudospectral method to solve the space-time tempered fractional diffusion equation.  SIAM Journal on Scientific Computing, 36(4):A1797–A1812, 2014.
  * [14] Zheyuan Hu, Ameya D Jagtap, George Em Karniadakis, and Kenji Kawaguchi.  When do extended physics-informed neural networks (XPINNs) improve generalization?  SIAM Journal on Scientific Computing (SISC), 44(5):A3158–A3182, 2022.
  * [15] Zheyuan Hu, Zekun Shi, George Em Karniadakis, and Kenji Kawaguchi.  Hutchinson trace estimation for high-dimensional and high-order physics-informed neural networks.  Computer Methods in Applied Mechanics and Engineering, 424:116883, 2024.
  * [16] Zheyuan Hu, Khemraj Shukla, George Em Karniadakis, and Kenji Kawaguchi.  Tackling the curse of dimensionality with physics-informed neural networks.  Neural Networks, 176:106369, 2024.
  * [17] Zheyuan Hu, Zhouhao Yang, Yezhen Wang, George Em Karniadakis, and Kenji Kawaguchi.  Bias-variance trade-off in physics-informed neural networks with randomized smoothing for high-dimensional PDEs.  arXiv preprint arXiv:2311.15283, 2023.
  * [18] Zheyuan Hu, Zhongqiang Zhang, George Em Karniadakis, and Kenji Kawaguchi.  Score-based physics-informed neural networks for high-dimensional Fokker-Planck equations.  arXiv preprint arXiv:2402.07465, 2024.
  * [19] George Em Karniadakis, editor.  Handbook of fractional calculus with applications, volume 3: Numerical Methods.  De Gruyter, Berlin, Boston, 2019.
  * [20] Kenji Kawaguchi, Leslie Pack Kaelbling, and Yoshua Bengio.  Generalization in deep learning.  arXiv preprint arXiv:1710.05468, 1(8), 2017.
  * [21] Diederik P Kingma and Jimmy Ba.  Adam: A method for stochastic optimization.  arXiv preprint arXiv:1412.6980, 2014.
  * [22] Nikolai Leonenko and Igor Podlubny.  Monte Carlo method for fractional-order differentiation extended to higher orders.  Fractional Calculus and Applied Analysis, 25(3):841–857, 2022.
  * [23] Can Li and Weihua Deng.  High order schemes for the tempered fractional diffusion equations.  Advances in computational mathematics, 42:543–572, 2016.
  * [24] Anna Lischke, Guofei Pang, Mamikon Gulian, Fangying Song, Christian Glusa, Xiaoning Zheng, Zhiping Mao, Wei Cai, Mark M. Meerschaert, Mark Ainsworth, and George Em Karniadakis.  What is the fractional Laplacian? a comparative review with new results.  Journal of Computational Physics, 404:109009, 2020.
  * [25] Lu Lu, Raphael Pestourie, Wenjie Yao, Zhicheng Wang, Francesc Verdugo, and Steven G Johnson.  Physics-informed neural networks with hard constraints for inverse design.  SIAM Journal on Scientific Computing, 43(6):B1105–B1132, 2021.
  * [26] Tao Luo and H. Yang.  Two-layer neural networks for partial differential equations: Optimization and generalization theory.  ArXiv, abs/2006.15733, 2020.
  * [27] Lei Ma, Fanhai Zeng, Ling Guo, George Em Karniadakis, et al.  Bi-orthogonal fPINN: A physics-informed neural network method for solving time-dependent stochastic fractional PDEs.  arXiv preprint arXiv:2303.10913, 2023.
  * [28] Siddhartha Mishra and Roberto Molinaro.  Estimates on the generalization error of physics informed neural networks (PINNs) for approximating PDEs.  arXiv preprint arXiv:2006.16144, 2020.
  * [29] Guofei Pang, Lu Lu, and George Em Karniadakis.  fPINNs: Fractional physics-informed neural networks.  SIAM Journal on Scientific Computing, 41(4):A2603–A2626, 2019.
  * [30] Maziar Raissi, Paris Perdikaris, and George E Karniadakis.  Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations.  Journal of Computational Physics, 378:686–707, 2019.
  * [31] Farzad Sabzikar, Mark M. Meerschaert, and Jinghua Chen.  Tempered fractional calculus.  Journal of Computational Physics, 293:14–28, 2015.  Fractional PDEs.
  * [32] Changtao Sheng, Bihao Su, and Chenglong Xu.  Efficient Monte Carlo method for integral fractional Laplacian in multiple dimensions.  SIAM Journal on Numerical Analysis, 61(5):2035–2061, 2023.
  * [33] Jing Sun, Daxin Nie, and Weihua Deng.  Algorithm implementation and numerical analysis for the two-dimensional tempered fractional Laplacian.  BIT Numerical Mathematics, 61(4):1421–1452, 2021.
  * [34] Shupeng Wang and George Em Karniadakis.  GMC-PINNs: A new general Monte Carlo PINNs method for solving fractional partial differential equations on irregular domains.  arXiv preprint arXiv:2405.00217, 2024.
  * [35] Lijing Zhao, Weihua Deng, and Jan S Hesthaven.  Spectral methods for tempered fractional differential equations.  Mathematics of Computation, 2016.


# Non-Markovian magnetization dynamics for uniaxial nanomagnets

Pascal Thibaudeau<EMAIL_ADDRESS>CEA DAM/Le Ripault, BP 16,
F-37260, Monts, FRANCE  Julien Tranchida<EMAIL_ADDRESS>CEA DAM/Le
Ripault, BP 16, F-37260, Monts, FRANCE  CNRS-Laboratoire de Mathématiques et
Physique Théorique (UMR 7350), Fédération de Recherche ”Denis Poisson”
(FR2964), Département de Physique, Université de Tours, Parc de Grandmont,
F-37200, Tours, FRANCE  Stam Nicolis<EMAIL_ADDRESS>CNRS-
Laboratoire de Mathématiques et Physique Théorique (UMR 7350), Fédération de
Recherche ”Denis Poisson” (FR2964), Département de Physique, Université de
Tours, Parc de Grandmont, F-37200, Tours, FRANCE

(August 27, 2024)

###### Abstract

A stochastic approach for the description of the time evolution of the
magnetization of nanomagnets is proposed, that interpolates between the
Landau-Lifshitz-Gilbert and the Landau-Lifshitz-Bloch approximations, by
varying the strength of the noise. Its finite autocorrelation time, i.e. when
it may be described as colored, rather than white, is, also, taken into
account and the consequences, on the scale of the response of the
magnetization are investigated. It is shown that the hierarchy for the moments
of the magnetization can be closed, by introducing a suitable truncation
scheme, whose validity is tested by direct numerical solution of the moment
equations and compared to the averages obtained from a numerical solution of
the corresponding colored stochastic Langevin equation. This comparison is
performed on magnetic systems subject to both an external uniform magnetic
field and an internal one-site uniaxial anisotropy.

##  I Introduction

Idealizations are often employed, when modeling real magnetization dynamics.
Far from equilibrium, magnetic fluctuations, modeled by white noise are
examples of such an idealization and the magnetization processes predicted in
this way have been observed in ultrafast reversal magnetization experiments.
However real fluctuations always have finite auto-correlation time and the
corresponding spin system dynamics is, rather, described by non-Markovian
stochastic processes atxitia2009ultrafast . Despite the practical relevance of
non-Markovian processes, however, these have received much less attention than
Markovian processes, chiefly, because of the simplicity of the theory and the
familiarity of the corresponding mathematical tools Gardiner:1998vf . The
study of the magnetization dynamics of super-paramagnetic systems, however, is
a suitable framework and provides motivation to describe magnetic
fluctuations, induced by contact to a heat bathGarcia-Palacios:1998uq , that
are sensitive to non–Markovian effects.

The effective magnetization dynamics of a spin system and a noise is described
by a stochastic differential equation of Langevin type, also known as the
stochastic Landau-Lifshitz-Gilbert (sLLG) equation :

$(1+\lambda^{2})\frac{\partial s_{i}}{\partial
t}=\epsilon_{ijk}(\omega_{j}+\Omega_{j})s_{k}+\lambda(\omega_{i}s_{j}s_{j}-\omega_{j}s_{j}s_{i}).$
(1)

The Einstein summation convention is adopted, with Latin indices standing for
vector components, $\epsilon_{ijk}$ is the anti-symmetric Levi-Civita pseudo-
tensor and $\bm{s}$ is a magnetization 3-vector, whose norm is fixed and can
be set to unity. Here $\bm{\omega}$ is the precession frequency which derives
from a spin Hamiltonian and is such connected to the spin vector. $\lambda$ is
the damping coefficient. The stochastic properties of ${\bf\Omega}$ describe
the fluctuations of external and/or internal degrees of freedom, random
perturbations and noise. Its correlation functions are assumed to satisfy the
relations

$\frac{\partial}{\partial
t}\langle\Omega_{i}(t)\dots\Omega_{i_{n}}(t_{n})\rangle=-\frac{1}{\tau}\langle\Omega_{i}(t)\dots\Omega_{i_{n}}(t_{n})\rangle$
(2)

with $\langle\rangle$ and $t\geq t_{1}\geq\dots\geq t_{n}$ is the time-ordered
statistical average product of stochastic variables Van-Kampen:1992it  and
$\tau$ is a single characteristic constant time. With proper conditions
Gaussian process, Poisson process and two-state (dichotomy) Markov process
enjoy this “time derivative property” Shapiro:1978qa . In general a detailed
description of the stochastic processes in the extended phase space $({\bf
s},{\bf\Omega})$ is not easy to establish Miyazaki:1998ek ; Thibaudeau:2012zr
.

Stochastic processes are completely described by an infinite set of
multivariate probability distributions. In practice, not all multivariate
distributions are not desirable but only a few, mainly the one-variable
distribution $P({\bf s},t)$ and the two-variable distribution, $P({\bf
s}_{1},t_{1};{\bf s}_{2},t_{2})$. The former is used to determine the average
$\langle f({\bf s},t)\rangle$ of any functional $f$ of the process ${\bf
s}(t)$, e.g. its statistical moments. The latter is used to determine two-time
statistics, e.g., the correlation function $\langle
s_{i}(t_{1})s_{j}(t_{2})\rangle$. In the white noise limit, a Fokker-Planck
equation on $P({\bf s},t)$ has been established only for a thermally
interacting spin system subject to a constant external field Brown-Jr:1963tp ;
Garcia-Palacios:1998uq . Beyond the mean-field approximation, the explicit
spin dependance of the effective field, including a magnetic anisotropy, on
the probability distribution $P$, may be addressed.

##  II Theory

In contact with a thermostat, described by the stochastic vector,
$\bm{\Omega}(t)$, the spin system precesses about the axis of magnetic
effective field, defined by both a uniform external magnetic field with
$\bm{\omega}^{(0)}=\gamma\mu_{0}\bm{H}$ and an anisotropic single-axis. This
can be observed in super-paramagnetic nano-magnets which behave as effective,
two-state, anisotropic spin systems Garcia-Palacios:1998uq . The deterministic
frequency vector is given by
$\bm{\omega}=\bm{\omega}^{(0)}+\kappa\bm{n}\left(\bm{n}\cdot\bm{s}\right)$,
where $\bm{n}$ is the anisotropy single-axis and $\kappa$ the anisotropy
frequency. The noise vector, $\bm{\Omega}(t)$, is assumed to be described by a
Gaussian stochastic process, that is completely defined by its 1-point
function $\langle\bm{\Omega}(t)\rangle$ and 2–point function
$\langle\Omega_{i}(t)\Omega_{j}(t^{\prime})\rangle$. Therefore the only
non–trivial information is contained in the these, that we assume to define an
Ornstein-Uhlenbeck process Uhlenbeck:1930kn , which means that
$\langle\Omega_{i}(t)\rangle=0$ and

$\langle\Omega_{i}(t)\Omega_{j}(t^{\prime})\rangle=\frac{D}{\tau}\exp(-\frac{|t-t^{\prime}|}{\tau})\delta_{ij}.$
(3)

It can be shown that this process does converge to the centered white noise
process, when $\tau\rightarrow 0$.

Moreover, it is assumed that the thermostat can be described by a fluctuation-
dissipation theorem that relates the strength of the noise, $D$, with the
temperature $T$ by the relationship $D=\lambda k_{b}T/\hbar$, where $\lambda$
is the parameter that appeared in Eq.(1).

These microscopic degrees of freedom give rise to average quantities, that
probe collective properties. Such quantities can be operationally defined by
taking statistical averages over the noise of Eq.(1) at equal times. In the
limit of low damping (i.e. $\lambda\ll 1$), the equation for the
magnetization, $\left\langle s_{i}\right\rangle$, is :

$\displaystyle\frac{\partial\langle s_{i}\rangle}{\partial t}$
$\displaystyle=$
$\displaystyle\epsilon_{ijk}(\langle\omega_{j}s_{k}\rangle+\langle\Omega_{j}s_{k}\rangle)$
(4)
$\displaystyle+\lambda(\langle\omega_{i}s_{j}s_{j}\rangle-\langle\omega_{j}s_{j}s_{i}\rangle),$

and implicates higher order correlation functions of noise and spin that have
to be determined in turn. By applying the Shapiro-Loginov formulae of
differentiation Shapiro:1978qa  to the components of the
$\langle{\Omega}_{i}s_{j}\rangle$ at equal times, we obtain the equation

$\frac{\partial\langle{\Omega}_{i}s_{j}\rangle}{\partial
t}=\left\langle{\Omega}_{i}\frac{\partial{s}_{j}}{\partial
t}\right\rangle-\frac{1}{\tau}\left\langle{\Omega}_{i}s_{j}\right\rangle.$ (5)

Introducing the right-hand side of Eq.(1) into Eq.(5), evolution equations for
the mixed moments $\langle{\Omega}_{i}s_{j}\rangle$ are obtained. Due to the
non-linearity of the sLLG equation, an infinite hierarchy arises
Nicolis:1998bd : the equations for the second-order moments depend on third-
order moments, and so on. In order to solve the system, closure relationships
need to be found. So the equations for the moments $\langle s_{i}s_{j}\rangle$
are now required and are provided by $\partial\langle
s_{i}s_{j}\rangle/\partial t=\langle
s_{i}\dot{s}_{j}\rangle+\langle\dot{s}_{i}{s}_{j}\rangle$. Assuming that this
identity holds, the time derivative of $s_{i}$ is replaced by the right-hand
side of Eq.(1) inside the statistical average. This leads to an additional set
of equations for the second-order moments of the microscopic spin degrees of
freedom, that are closely related to those written by Garanin et al (see
Eq.(6)Garanin:1990fk ) and Garcia-Palacios et al. (see Eq.(2.10) Garcia-
Palacios:1998uq ). As long as the stochastic variables follow individual
Gaussian processes, their connected moments (labelled as
$\langle\\!\langle.\rangle\\!\rangle$), beyond the second moment are zero Van-
Kampen:1992it . The case of the mixed connected moments of spin and noise is
more subtle. If a linear coupling between the spin and the noise is assumed,
the mixed cumulants, above the second will, also, vanish. If the Gaussian
hypothesis is dropped, the general relations, at equal times, between the
desired moments are :

$\displaystyle\langle\Omega_{i}s_{j}s_{k}\rangle$ $\displaystyle=$
$\displaystyle\langle\Omega_{i}s_{j}\rangle\langle
s_{k}\rangle+\langle\Omega_{i}s_{k}\rangle\langle
s_{j}\rangle+\langle\\!\langle\Omega_{i}s_{j}s_{k}\rangle\\!\rangle,$ (6)
$\displaystyle\langle\Omega_{i}s_{j}s_{k}s_{l}\rangle$ $\displaystyle=$
$\displaystyle\langle\\!\langle\Omega_{i}s_{j}s_{k}s_{l}\rangle\\!\rangle+\langle\Omega_{i}s_{j}\rangle\langle
s_{k}s_{l}\rangle+\langle\Omega_{i}s_{k}\rangle\langle s_{j}s_{l}\rangle$ (7)
$\displaystyle+\langle\Omega_{i}s_{l}\rangle\langle s_{j}s_{k}\rangle+\langle
s_{l}\rangle\langle\\!\langle\Omega_{i}s_{j}s_{k}\rangle\\!\rangle$
$\displaystyle+\langle
s_{k}\rangle\langle\\!\langle\Omega_{i}s_{j}s_{l}\rangle\\!\rangle+\langle
s_{j}\rangle\langle\\!\langle\Omega_{i}s_{k}s_{l}\rangle\\!\rangle,$
$\displaystyle\langle s_{i}s_{j}s_{k}\rangle$ $\displaystyle=$
$\displaystyle\langle s_{i}s_{j}\rangle\langle s_{k}\rangle+\langle
s_{j}s_{k}\rangle\langle s_{i}\rangle+\langle s_{i}s_{k}\rangle\langle
s_{j}\rangle$ (8) $\displaystyle-2\langle s_{i}\rangle\langle
s_{j}\rangle\langle s_{k}\rangle+\langle\\!\langle
s_{i}s_{j}s_{k}\rangle\\!\rangle,$ $\displaystyle\langle
s_{i}s_{j}s_{k}s_{l}\rangle$ $\displaystyle=$ $\displaystyle\langle
s_{i}s_{j}\rangle\langle s_{k}s_{l}\rangle+\langle s_{i}s_{k}\rangle\langle
s_{j}s_{l}\rangle+\langle s_{i}s_{l}\rangle\langle s_{j}s_{k}\rangle$ (9)
$\displaystyle-2\langle s_{i}\rangle\langle s_{j}\rangle\langle
s_{k}\rangle\langle s_{l}\rangle+\langle s_{k}\rangle\langle\\!\langle
s_{i}s_{j}s_{l}\rangle\\!\rangle$ $\displaystyle+\langle
s_{j}\rangle\langle\\!\langle s_{i}s_{k}s_{l}\rangle\\!\rangle+\langle
s_{l}\rangle\langle\\!\langle s_{i}s_{j}s_{k}\rangle\\!\rangle$
$\displaystyle+\langle s_{i}\rangle\langle\\!\langle
s_{j}s_{k}s_{l}\rangle\\!\rangle+\langle\\!\langle
s_{i}s_{j}s_{k}s_{l}\rangle\\!\rangle.$

These considerations lead to the following system of equations :

$\displaystyle\frac{\partial\langle s_{i}\rangle}{\partial t}$
$\displaystyle=$ $\displaystyle\epsilon_{ijk}(\omega^{(0)}_{j}\langle
s_{k}\rangle+\kappa n_{j}n_{l}\langle
s_{l}s_{k}\rangle+\langle\Omega_{j}s_{k}\rangle)$ (10)
$\displaystyle+\lambda\left(\omega_{i}^{(0)}\langle
s_{j}s_{j}\rangle-\omega_{j}^{(0)}\langle s_{j}s_{i}\rangle\right.$
$\displaystyle+\big{.}\kappa n_{l}(n_{i}\langle s_{l}s_{j}s_{j}\rangle-
n_{j}\langle s_{l}s_{j}s_{i}\rangle)\Big{)},$
$\displaystyle\frac{\partial\langle\Omega_{j}s_{k}\rangle}{\partial t}$
$\displaystyle=$
$\displaystyle-\frac{1}{\tau}\langle\Omega_{j}s_{k}\rangle+\epsilon_{klm}\left(\langle\Omega_{j}\omega_{l}s_{m}\rangle+\langle\Omega_{j}\Omega_{l}s_{m}\rangle\right),$
(11)
$\displaystyle+\lambda\left(\langle\Omega_{j}\omega_{k}s_{l}s_{l}\rangle-\langle\Omega_{j}\omega_{l}s_{l}s_{k}\rangle\right),$
$\displaystyle\frac{\partial\langle s_{i}s_{j}\rangle}{\partial t}$
$\displaystyle=$
$\displaystyle\epsilon_{ikl}\left(\langle\omega_{k}s_{l}s_{j}\rangle+\langle\Omega_{k}s_{l}s_{j}\rangle\right)$
(12)
$\displaystyle+\lambda\left(\langle\omega_{i}s_{k}s_{k}s_{j}\rangle-\langle\omega_{k}s_{k}s_{i}s_{j}\rangle\right)$
$\displaystyle+(i\leftrightarrow j).$

These expressions are closed by considering first the independence of Gaussian
processes, the form of the effective frequency vector and the expression of
the autocorrelation function of the noise. One has explicitly :

$\displaystyle\langle\Omega_{j}\omega_{l}s_{m}\rangle$ $\displaystyle=$
$\displaystyle\omega^{(0)}_{l}\langle\Omega_{j}s_{m}\rangle$
$\displaystyle+\kappa n_{l}n_{u}\left(\langle\Omega_{j}s_{u}\rangle\langle
s_{m}\rangle+\langle\Omega_{j}s_{m}\rangle\langle s_{u}\rangle\right)$
$\displaystyle\langle\Omega_{j}\Omega_{l}s_{m}\rangle$ $\displaystyle=$
$\displaystyle\frac{D}{\tau}\delta_{jl}\langle s_{m}\rangle$
$\displaystyle\langle\Omega_{j}\omega_{k}s_{l}s_{l}\rangle$ $\displaystyle=$
$\displaystyle 2\omega^{(0)}_{k}\langle\Omega_{j}s_{l}\rangle\langle
s_{l}\rangle$ $\displaystyle+\kappa
n_{k}n_{u}\left(\langle\Omega_{j}s_{u}\rangle\langle
s_{l}s_{l}\rangle+2\langle\Omega_{j}s_{l}\rangle\langle
s_{u}s_{l}\rangle\right)$
$\displaystyle\langle\Omega_{j}\omega_{l}s_{l}s_{k}\rangle$ $\displaystyle=$
$\displaystyle\omega^{(0)}_{l}\left(\langle\Omega_{j}s_{l}\rangle\langle
s_{k}\rangle+\langle\Omega_{j}s_{k}\rangle\langle s_{l}\rangle\right)$
$\displaystyle+\kappa n_{l}n_{u}\left(\langle\Omega_{j}s_{u}\rangle\langle
s_{l}s_{k}\rangle+\langle\Omega_{j}s_{l}\rangle\langle
s_{u}s_{k}\rangle\right.$
$\displaystyle+\left.\langle\Omega_{j}s_{k}\rangle\langle
s_{u}s_{l}\rangle\right)$ $\displaystyle\langle\omega_{k}s_{l}s_{j}\rangle$
$\displaystyle=$ $\displaystyle\omega^{(0)}_{k}\langle
s_{l}s_{j}\rangle+\kappa n_{k}n_{u}\left(\langle s_{u}\rangle\langle
s_{l}s_{j}\rangle\right.$ $\displaystyle\left.+\langle s_{l}\rangle\langle
s_{u}s_{j}\rangle+\langle s_{j}\rangle\langle s_{u}s_{l}\rangle-2\langle
s_{u}\rangle\langle s_{l}\rangle\langle s_{j}\rangle\right)$
$\displaystyle\langle\Omega_{k}s_{l}s_{j}\rangle$ $\displaystyle=$
$\displaystyle\langle\Omega_{k}s_{l}\rangle\langle
s_{j}\rangle+\langle\Omega_{k}s_{j}\rangle\langle s_{l}\rangle$
$\displaystyle\langle\omega_{i}s_{k}s_{k}s_{j}\rangle$ $\displaystyle=$
$\displaystyle\omega^{(0)}_{i}\left(\langle s_{j}\rangle\langle
s_{k}s_{k}\rangle+2\langle s_{k}\rangle\langle s_{k}s_{j}\rangle\right.$
$\displaystyle\left.-2\langle s_{k}\rangle\langle s_{k}\rangle\langle
s_{j}\rangle\right)+\kappa n_{i}n_{u}\left(\langle s_{k}s_{k}\rangle\right.$
$\displaystyle\left.\langle s_{u}s_{j}\rangle+2\langle
s_{u}s_{k}\rangle\langle s_{k}s_{j}\rangle-2\langle s_{u}\rangle\langle
s_{k}\rangle^{2}\langle s_{j}\rangle\right)$
$\displaystyle\langle\omega_{k}s_{k}s_{i}s_{j}\rangle$ $\displaystyle=$
$\displaystyle\omega^{(0)}_{k}\left(\langle s_{j}\rangle\langle
s_{k}s_{i}\rangle+\langle s_{k}\rangle\langle s_{i}s_{j}\rangle\right.$
$\displaystyle\left.+\langle s_{k}s_{j}\rangle\langle s_{i}\rangle-2\langle
s_{k}\rangle\langle s_{i}\rangle\langle s_{j}\rangle\right)$
$\displaystyle+\kappa n_{k}n_{u}\left(\langle s_{u}s_{k}\rangle\langle
s_{i}s_{j}\rangle+\langle s_{u}s_{i}\rangle\langle s_{k}s_{j}\rangle\right.$
$\displaystyle\left.+\langle s_{u}s_{j}\rangle\langle
s_{k}s_{i}\rangle-2\langle s_{u}\rangle\langle s_{k}\rangle\langle
s_{i}\rangle\langle s_{j}\rangle\right)$

Combined with Eqs.(10,11,12), the derived system of equations is called the
dynamical Landau-Lifshitz-Bloch model (dLLB). In the limit of $\tau\rightarrow
0$, Eq.(11) gives a direct expression for the mixed moment of the noise and
the magnetization as $\langle\Omega_{i}s_{j}\rangle=-D\epsilon_{ijk}\langle
s_{k}\rangle$, which makes the $\langle\Omega_{i}s_{j}\rangle$ matrix
antisymmetric. Once introduced into Eq.(10), a longitudinal contribution for
the spin dynamics appears in which the amplitude of the noise is proportional
to the square root of the temperature. For finite values of $\tau$, Eq.(10)
exhibits a ”memory effect” which shows that the value of $\langle
s_{i}(t)\rangle$ for the time $t$ depends on all the previous values of
$\langle s_{i}(t^{\prime})\rangle$ before the time $t$.

##  III Numerical experiments

Numerical experiments were performed to test the consistency of the closure
hypothesis. For “low” values of $\tau$, the dynamics is close to Markovian and
Eqs.(10,11,12) are expected to match the results obtained from an averaged
Markovian stochastic dynamics. In order to realize these comparisons, a
reference model is constructed from Eq.(1) with a white noise random vector
$\bm{\Omega}$. Multiple spin trajectories with the same initial conditions are
generated and followed in time according to a previously detailed symplectic
integration scheme Lewis:2003lq ; Thibaudeau:2012zr ; Beaujouan:2012zk . A
statistical average is taken over all the realizations. Eqs.(10,11,12) are
integrated with an eight-order Runge-Kutta integration scheme, with adaptive
steps, constructed from the Jacobian of the whole system.

The figure (1-left) (resp. 1-right) exhibits the average magnetization (resp.
spin-spin autocorrelation function) from the stochastic LLG equation and our
dLLB model for the same value of the amplitude $D$ in the white noise limit.

Figure 1: Dynamics of the averaged magnetization (left) and auto-correlation
spin function (right) around a single-axis in the $z$-direction in the case of
small correlation time of the noise. The dot lines plot the sLLG case
($10^{3}$ repetitions), solid lines stand for our dLLG model with the Gaussian
hypothesis of independance. $D=1\times
10^{-3}{\rm{GHz}};\lambda=0.1;\tau=10^{-6}{\rm{ns}};{\bf{\omega}}^{(0)}=(0,0,0);\kappa=1$GHz;
${\bf n}=(0,0,1)$ are taken with the initial conditions ${\bf
s}(0)=(\frac{1}{\sqrt{3}},\frac{1}{\sqrt{3}},\frac{1}{\sqrt{3}})$, $\forall
i,j$ $\langle s_{i}(0)s_{j}(0)\rangle=\frac{1}{3}$ and
$\langle\Omega_{i}(0)s_{j}(0)\rangle=0$.

For low value of $D$ (i.e. in the normal phase), the average stochastic
dynamics follows our dLLB model. A similar situation has been already
investigated for an external constant field only Tranchida:hl . When $\tau$
increases, the noise is no longer white and does not follow a Markov’s
process. A complete stochastic system of equations for both the spin and the
noise is built by supplementing Eq.(1) with

$\frac{d\Omega_{i}}{dt}=-\frac{1}{\tau}\left(\Omega_{i}-\xi_{i}(t)\right)$
(13)

where $\xi_{i}$ is a Gaussian centered stochastic process Fox:1988tg , but
with $\langle\xi_{i}(t)\rangle=0$ and
$\langle\xi_{i}(t)\rangle\xi_{j}(t^{\prime})\rangle=2D\delta_{ij}\delta(t-t^{\prime})$.
Formally, this system is integrated as $\Omega_{i}(t)\equiv
e^{tL_{\Omega}}\Omega_{i}(0)=\int_{0}^{t}e^{-(t-u)/\tau}\xi_{i}(u)du$ with
$L_{\Omega}$ is the Liouville operator corresponding to Eq.(13). Because the
Liouville operators of the spin and the noise commute, an integration of the
colored system is performed by the combination of the two operators as
$e^{tL}(s_{i}(0),\Omega_{i}(0))=e^{tL_{s}}e^{tL_{\Omega}}(s_{i}(0),\Omega_{i}(0))=(s_{i}(t),\Omega_{i}(t))$.
This is performed for larger values of $\tau$ for which the figure (2)
exhibits $\langle s_{i}(t)\rangle$ and the diagonal part of $\langle
s_{i}(t)s_{j}(t)\rangle$.

Figure 2: Stochastic average magnetization (left) and auto-correlation spin
function (right) around a single-axis in the $z$-direction and a field in the
$x$-direction. The dot lines are for the sLLG for $10^{4}$ repetitions and the
solid lines stand for our dLLB model. $D=1\times
10^{-3}{\rm{GHz}};\lambda=0.05;{\bf{\omega}}^{(0)}=(\pi/5,0,0)$GHz ;
$\kappa=1/10$GHz; ${\bf n}=(0,0,1)$; $\tau=10^{-1}$ ns are taken with the same
initial conditions as previously.

It is observed that, even for “large” values of $\tau$, our dLLB model
continues to track the dynamics of the colored stochastic LLG system. The
influence of the higher than second-order cumulants on the dynamics of the
average magnetization and spin autocorrelation function is mitigated by a
large damping factor $\lambda$ which prevents detailed descriptions of the
magnetization dynamics. Moreover, because of the link between the amplitude
$D$ and the damping factor $\lambda$, through a to-be-derived fluctuation-
dissipation theorem in the special case of an explicit dependance of the spin
in the effective field, the domain of validity of our dLLB model in
temperature could be significantly controlled.

When the intensity of the noise increases, the spin dynamics is no longer
Gaussian. It is observed numerically, in the stochastic dynamics, by
calculating the cumulant function $\langle\\!\langle
s_{i}(t)s_{j}(t)s_{k}(t)\rangle\\!\rangle$. The ten independent spin cumulants
are depicted, for several values of the amplitude $D$ in the figure (3).

Figure 3: Stochastic average over for $10^{5}$ repetitions of the independent
$\langle\\!\langle s_{i}s_{j}s_{k}\rangle\\!\rangle$ terms for various values
of $D$. $\lambda=0.1;{\bf{\omega}}^{(0)}=(\pi/5,0,0)$GHz ; $\kappa=0.5$GHz;
${\bf n}=(0,0,1)$; $\tau=10^{-3}$ns are taken with the same initial conditions
as previously.

As shown, for this particular situation, some terms go to zero in time,
whereas some remain finite, depending on the underlying symmetries. At the
expense of higher order expressions, the Eqs. (10,11,12) can be supplemented
by an expansion of $\frac{d\langle s_{i}s_{j}s_{k}\rangle}{dt}$ and integrated
together. This represents an option where the deviation off the Gaussian
character of the spin distribution is explicitly treated at the third-order.

The third-order correlation functions for the dynamical variables ${\bf s}$
satisfying a stochastic Langevin-like equation are generally determined by
first solving this equation in the presence of the noise and then averaging
the product of these variables over the ensemble of the random noise. A path
integral formalism can be established for deducing identities between the
correlation functions. This strategy has been tested in various approximation
schemes, including the quenched approximation (i.e. $D=0$). A shortcut might
be to look for an alternate description of an effective system of equations in
a closed time path formalism of Schwinger allowing to describe a general non-
equilibrium field theory at finite temperature Kleinert:2006ye . This is why a
closure approximation was proposed by Martin, Siggia and Rose Martin:1973lr
which represents the classical limit of the closed time path formalism and
will be the subject of further investigation.

JT acknowledges financial support through a joint doctoral fellowship “Région
Centre-CEA” under grant number 00086667.

## References

  * (1) U. Atxitia, O. Chubykalo-Fesenko, R. W. Chantrell, U. Nowak, and A. Rebei.  Ultrafast spin dynamics: the effect of colored noise.  Phys. Rev. Lett., 102(5):057203, 2009.
  * (2) C. W. Gardiner.  Handbook of Stochastic Methods in Physics, Chemistry and Natural Sciences.  Springer-Verlag, Berlin–Heidelberg–New York–Tokyo, 2nd edition, 1998\.
  * (3) J. L. García-Palacios and F. J. Lázaro.  Langevin-dynamics study of the dynamical properties of small magnetic particles.  Phys. Rev. B, 58(22):14937, 1998.
  * (4) N. G. Van Kampen.  Stochastic processes in physics and chemistry, volume 1.  Elsevier, 1992.
  * (5) V. E. Shapiro and V. M. Loginov.  Formulae of differentiation and their use for solving stochastic equations.  Physica A: Statistical Mechanics and its Applications, 91(3):563–574, 1978.
  * (6) K. Miyazaki and K. Seki.  Brownian motion of spins revisited.  J. Chem. Phys., 108(17):7052–7059, 1998.
  * (7) P. Thibaudeau and D. Beaujouan.  Thermostatting the atomic spin dynamics from controlled demons.  Physica A: Statistical Mechanics and its Applications, 391:1963–1971, 2012.
  * (8) W. F. Brown Jr.  Thermal fluctuations of a single-domain particle.  J. Appl. Phys., 34(4):1319–1320, 1963.
  * (9) G. E. Uhlenbeck and L. S. Ornstein.  On the theory of the brownian motion.  Phys. Rev., 36:823–841, 1930.
  * (10) C. Nicolis and G. Nicolis.  Closing the hierarchy of moment equations in nonlinear dynamical systems.  Phys. Rev. E, 58(4):4391, 1998.
  * (11) D. A. Garanin, V. V. Ishchenko, and L. V. Panina.  Dynamics of an ensemble of single-domain magnetic particles.  Theor. Math. Phys., 82(2):169–179, 1990.
  * (12) D. Lewis and N. Nigam.  Geometric integration on spheres and some interesting applications.  J. Comp. Appl. Math., 151:141–170, 2003.
  * (13) D. Beaujouan, P. Thibaudeau, and C. Barreteau.  Anisotropic magnetic molecular dynamics of cobalt nanowires.  Phys. Rev. B, 86(17):174409, 2012.
  * (14) J. Tranchida, P. Thibaudeau, and S. Nicolis.  Closing the hierarchy for non-markovian magnetization dynamics.  Physica B, 2015.  arXiv:1506.00544 (in press).
  * (15) R. F. Fox, I. R. Gatland, R. Roy, and G. Vemuri.  Fast, accurate algorithm for numerical simulation of exponentially correlated colored noise.  Phys. Rev. A, 38(11):5938–5940, 1988.
  * (16) H. Kleinert.  Path integrals in quantum mechanics, statistics, polymer physics, and financial Markets.  World Scientific, Singapore, 5th edition, 2006.
  * (17) P. C. Martin, E. D. Siggia, and H. A. Rose.  Statistical dynamics of classical systems.  Phys. Rev. A, 8(1):423–437, 1973.