# Fractional Brownian motion in confining potentials: non-equilibrium
distribution tails and optimal fluctuations

Baruch Meerson<EMAIL_ADDRESS>Racah Institute of Physics, Hebrew
University of Jerusalem, Jerusalem 91904, Israel  Pavel V. Sasorov
<EMAIL_ADDRESS>ELI Beamlines Facility, Extreme Light infrastructure
ERIC, 252 41 Dolni Brezany, Czech Republic

###### Abstract

At long times, a fractional Brownian particle in a confining external
potential reaches a non-equilibrium (non-Boltzmann) steady state. Here we
consider scale-invariant power-law potentials $V(x)\sim|x|^{m}$, where $m>0$,
and employ the optimal fluctuation method (OFM) to determine the large-$|x|$
tails of the steady-state probability distribution $\mathcal{P}(x)$ of the
particle position. The calculations involve finding the optimal (that is, the
most likely) path of the particle, which determines these tails, via a
minimization of the exact action functional for this system, which has
recently become available. Exploiting dynamical scale invariance of the model
in conjunction with the OFM ansatz, we establish the large-$|x|$ tails of
$\ln\mathcal{P}(x)$ up to a dimensionless factor $\alpha(H,m)$, where $0<H<1$
is the Hurst exponent. We determine $\alpha(H,m)$ analytically (i) in the
limits of $H\to 0$ and $H\to 1$, and (ii) for $m=2$ and arbitrary $H$,
corresponding to the fractional Ornstein-Uhlenbeck (fOU) process. Our results
for the fOU process are in agreement with the previously known exact
$\mathcal{P}(x)$ and autocovariance. The form of the tails of $\mathcal{P}(x)$
yields exact conditions, in terms of $H$ and $m$, for the particle confinement
in the potential. For $H\neq 1/2$, the tails encode the non-equilibrium
character of the steady state distribution, and we observe violation of time
reversibility of the system except for $m=2$. To compute the optimal paths and
the factor $\alpha(H,m)$ for arbitrary permissible $H$ and $m$, one needs to
solve an (in general nonlinear) integro-differential equation. To this end we
develop a specialized numerical iteration algorithm which accounts
analytically for an intrinsic cusp singularity of the optimal paths for
$H<1/2$.

##  I Introduction

The Mandelbrot -– van Ness fractional Brownian motion (fBm) Kolmogorov ;
Mandelbrot ; Qian2003 , a Gaussian non-Markovian process with stationary
increments, provides a prototypical example of a system with long-range
temporal correlations and anomalous diffusion. Here we focus on the following
question: how do the properties of the fBm modify in the presence of
confinement? A natural model system for studying this type of question
involves overdamped stochastic motion of a fractional Brownian particle in an
external potential Metzler2021 ; Metzler2022 ; MBO ; Metzler2023 . Here we
consider a family of scale-invariant power-law potentials $V(x)=k|x|^{m}/m$,
where $m>0$. This brings about a non-Markovian Langevin equation

$\dot{x}=f[x(t)]+\xi(t)$ (1)

with the deterministic force

$f(x)=-\frac{dV(x)}{dx}=-k|x|^{m-1}\,\text{sign}\,x$ (2)

and the long-range correlated fractional Gaussian noise (fGn) $\xi(t)$. In the
absence of external potential, the particle position $x(t)$ in Eq. (1) follows
the fractional Brownian motion (fBm), introduced by Kolmogorov Kolmogorov  and
by Mandelbrot and van Ness (MvN) Mandelbrot . Like the standard Brownian
motion, the MvN fBm is a scale-invariant Gaussian process. It is defined, for
all $|t|<\infty$, by its mean (which can be set to zero) and the
autocovariance Mandelbrot

$\\!\kappa(t,t^{\prime})=\langle
x(t)x(t^{\prime})\rangle=D\left(|t|^{2H}\\!+\\!|t^{\prime}|^{2H}\\!-\\!|t-t^{\prime}|^{2H}\right),$
(3)

where $0<H<1$ is the Hurst exponent, and $D$ is the generalized diffusion
constant with units $\text{length}^{2}/\text{time}^{2H}$. The standard
Brownian motion corresponds to $H=1/2$. For $H\neq 1/2$, the MvN process is
non-Markovian and long-range correlated. It describes subdiffusion and
superdiffusion for $H<1/2$ and $H>1/2$, respectively. For $H=1/2$ the standard
Brownian motion is recovered.

The fGn – the time derivative of the MvN fBm – is a stationary Gaussian
process, and it is also scale invariant. Its autocovariance
$c(t_{1}-t_{2})=\langle\xi(t_{1})\xi(t_{2})\rangle$ is the following
Mandelbrot ; Qian2003 :

$c(\tau)=2DH\frac{d}{d\tau}\left(|\tau|^{2H-1}\text{sgn}\,\tau\right)\equiv
D\frac{d^{2}}{d\tau^{2}}|\tau|^{2H}\,.$ (4)

In the particular case of $H=1/2$, Eq. (4) yields

$c(\tau)=D\frac{d}{d\tau}\,\text{sgn}\,\tau=2D\delta(\tau)\,,$ (5)

the familiar expression for the white noise.

In the presence of the confining potential $V(x)$, the fractional Brownian
particle, as described by Eq. (1), reaches a steady state. Unless $H=1/2$,
this is a _non-equilibrium_ steady state (NESS). In particular, the steady-
state probability distribution $\mathcal{P}(X)$ of observing the particle at
the position $x=X$ is _not_ given by the Boltzmann formula Metzler2021 . This
and other nonequilibrium features, that we discuss here, make this model
interesting not also in particular applications (see e.g. Ref. Metzler2021
for an extensive list of such), but also in the more general context of
statistical mechanics far from equilibrium.

For $H\neq 1/2$, the steady-state distribution $\mathcal{P}(X)$ is known
exactly only in the special case of the quadratic potential, $m=2$ Metzler2021
. Here we focus on the large-$|X|$ tails of $\mathcal{P}(X)$ for all $H$ and
$m$ that allow for the particle confinement. These (identical) tails describe
large deviations of the particle position from its mean value $X=0$.
Exploiting the corresponding large parameter, we apply the optimal fluctuation
method, which relies on finding the optimal (that is, the most likely) path of
the particle which determines the tails. The calculations involve a
minimization of the proper action functional for Eq. (1), recently obtained in
Ref. MBO .

Using scale invariance properties of the model in conjunction with the OFM,
one can immediately establish the large-$|X|$ asymptotic of the logarithm of
$\mathcal{P}(X)$ up to an unknown dimensionless factor $\alpha(H,m)$. We find
$\alpha(H,m)$ analytically (i) in the limits of $H\to 0$ and $H\to 1$, and
(ii) for $m=2$ and arbitrary $H$. The latter result agrees with the previously
known exact expression for $\mathcal{P}(x)$ for $m=2$ Metzler2021 . The
knowledge of the tail asymptotic provides exact conditions for the confinement
of the particle in terms of $H$ and $m$, which agree with previous numerical
results Metzler2022 . As an additional manifestation of the nonequilibrium
nature of the steady state, we also observe violation of the time
reversibility of the systems’ paths, except for $H=1/2$ or $m=2$.

The minimization of the action functional MBO  can be reduced to solving an
(in general nonlinear) integro-differential equation for the optimal path
conditioned on reaching a specified point $x=X$. We developed a specialized
numerical iteration algorithm for this purpose. The algorithm accounts
analytically for the intrinsic cusp singularities of the optimal paths and, in
principle, makes it possible to compute the dimensionless factor $\alpha(H,m)$
for any permissible $H$ and $m$.

Here is a plan of the rest of the paper. In Sec. II we revisit and slightly
extend the path integral formulation MBO  for Eq. (1), introduce the OFM,
establish the scaling behavior of the distribution tails, obtain the
confinement condition, exactly solve the OFM problem for the quadratic
potential $m=2$, determine the factor $\alpha(H,M)$ in the limits of $H\to 0$
and $H\to 1$ for any permissible $m$ and demonstrate the violation of the time
reversibility of the system’s paths for $H\neq 1/2$. In Sec. III we present
our numerical algorithm for computing the optimal paths and evaluating the
function $\alpha(H,m)$. A brief summary of main results is presented in Sec.
IV. Some technical details of the numerical algorithm are relegated to the two
appendices.

##  II Optimal Fluctuation Method

###  II.1 Path integral

Here we will exploit exact path-integral representations for Eq. (1) which
have been recently derived separately for the subdiffusion, $0<H<1/2$, and the
superdiffusion, $1/2<H<1$ MBO . For the subdiffusion, the action functional of
the path integral reads

$\mathcal{S}[x(t)]=\frac{1}{2}\int_{-\infty}^{\infty}dt_{1}\int_{-\infty}^{\infty}dt_{2}\,C(t_{1}-t_{2})[\dot{x}-f(x)]_{t_{1}}[\dot{x}-f(x)]_{t_{2}}\,,$
(6)

where the subscripts $t_{1}$ and $t_{2}$ denote the arguments of the functions
inside the square brackets. In its turn,

$C(\tau)=\frac{\text{cot}\,(\pi H)}{4\pi DH}\frac{1}{|\tau|^{2H}}$ (7)

is the inverse kernel,

$\int_{-\infty}^{\infty}dt^{\prime\prime}\,C(t-t^{\prime\prime})c(t^{\prime}-t^{\prime\prime})=\delta(t-t^{\prime})\,,$
(8)

for the autocovariance $c(\tau)$ of the fGn, see Eq. (4).

For the superdiffusion Ref. MBO  gives the action functional

$\mathcal{S}[x(t)]=\frac{1}{2}\int_{-\infty}^{\infty}dt_{1}\int_{-\infty}^{\infty}dt_{2}\,Q(t_{1}-t_{2})\left[\ddot{x}-f^{\prime}[x(t)]\,\dot{x}(t)\right]_{t_{1}}\left[\ddot{x}(t)-f^{\prime}[x(t)]\,\dot{x}(t)\right]_{t_{2}},$
(9)

where

$Q(\tau)=\frac{{\rm cot}(\pi H)}{8\pi DH(1-H)(2H-1)}|\tau|^{2-2H}\,,$ (10)

and $f^{\prime}(x)\equiv df/dx$.

In fact, both expressions (6) and (9) can be represented as Eq. (6) with the
universal kernel

$C(t_{1}-t_{2})=\frac{{\rm cot}(\pi
H)\,\partial_{t_{1}}\\!\partial_{t_{2}}|t_{1}-t_{2}|^{2-2H}}{8\pi
DH(1-H)(2H-1)}\,.$ (11)

For $H<1/2$ the differentiations in Eq. (11) involve classical functions. For
$H>1/2$, the differentiations should be interpreted in terms of generalized
functions, or Schwartz distributions.

It is useful for the following to also present the governing equations in the
frequency domain. The action for a given realization of the fGn $\xi(t)$ can
be written as

$S[\xi(t)]=\frac{1}{4\pi}\int\limits_{-\infty}^{\infty}C_{\omega}\xi_{\omega}\xi_{\omega}^{*}\,d\omega\,,$
(12)

where $\xi_{\omega}=\int_{-\infty}^{\infty}dt\,\xi(t)e^{i\omega t}$ is the
Fourier transform of the fGn,

$C_{\omega}=\int_{-\infty}^{\infty}dt\,C(t)e^{i\omega
t}=\frac{|\omega|^{2H-1}}{2D\sin(\pi H)\Gamma(2H+1)}$ (13)

is the Fourier transform of the inverse kernel $C(t)$, and $\Gamma(\dots)$ is
the gamma function. The expression (12) for the action with this positive
$C_{\omega}$ is well defined for any $0<H<1$.

The Fourier transform of the autocovariance itself is

$c_{\omega}=\int_{-\infty}^{\infty}dt\,c(t)e^{i\omega t}=DH\sin(\pi
H)\Gamma(2H)|\omega|^{1-2H}\,,$ (14)

as Eq. (8) becomes simply $c_{\omega}\,C_{\omega}=1$. As a result, we can also
express in the frequency domain the action (6):

$\mathcal{S}[x(t)]=\frac{1}{4\pi}\int_{-\infty}^{\infty}d\omega\,C_{\omega}[\dot{x}-f(x)]_{\omega}\,[\dot{x}-f(x)]_{\omega}^{*}\,.$
(15)

When $H\to 1/2$, the action functional (6) reduces to the familiar local
action for the white-noise driving:

$\mathcal{S}[x(t)]=\frac{1}{4D}\int_{-\infty}^{\infty}dt\left[\dot{x}(t)-f(x(t))\right]^{2}\,.$
(16)

###  II.2 OFM formulation

Sufficiently far in the tail, $\mathcal{P}(X)$ is dominated by a single
_optimal path_ , for which we will keep the notation $x(t)$. The optimal path
minimizes the action (6) subject to the conditions

$x(t=0)=X$ (17)

and

$x(t\to-\infty)=x(t\to\infty)=0.$ (18)

It is convenient to incorporate the condition (17) by introducing the modified
action

$\mathcal{S}_{\lambda}[x(t)]=\mathcal{S}[x(t)]-\lambda\int_{-\infty}^{\infty}dt\,x(t)\delta(t)\,.$
(19)

where $\mathcal{S}[x(t)]$ is given by Eq. (6). The Lagrange multiplier
$\lambda$ should be ultimately expressed through $X$ and the rest of the
parameters of the problem.

Demanding that the linear variation of the modified action vanish, and using
the symmetry condition $C(-\tau)=C(\tau)$ and two integrations by part, we
arrive at the following equation:

$\int_{-\infty}^{\infty}dt^{\prime}C(t-t^{\prime})\left[f^{\prime}(x(t))+\partial_{t^{\prime}}\right]\left[\dot{x}-f(x)\right]_{t^{\prime}}=-\lambda\delta(t)$
(20)

which, alongside with the conditions (17) and (18), describe the optimal path
$x(t)$. In general, the integro-differential equation (20) is nonlinear. It is
linear, however, for the quadratic potential $m=2$, which will be considered
in Sec. II.5.

Once the optimal path $x(t)$ is determined, one can evaluate, up to an unknown
pre-exponential factor, the position distribution $\mathcal{P}(X)$ by using
the approximate relation $-\ln\mathcal{P}(X)\simeq S(X)$, where $S(X)$ is the
action (6) evaluated along the optimal path. A useful shortcut in the
calculation of the action is provided by the relation

$\frac{dS}{dX}=\lambda(X)\,,$ (21)

which follows from the fact that $X$ and $\lambda$ are conjugate variables,
see, e.g. Ref. Cunden . Equation (21) enables one to immediately calculate the
action once the Lagrange multiplier $\lambda$ is expressed through $X$, or
vice versa.

###  II.3 Scaling behavior

The scale invariance properties of the fGn and of the power-law potential
$V(x)=k|x|^{m}/m$ make it possible to immediately establish some important
scaling properties of the action $S(X)$ by using dimensional argument
Barenblatt . Indeed, the dimensionless action $S(X)\simeq-\ln\mathcal{P}(X)$
depends only on the dimensional parameters $X$, $k$ and $D$ and on the
dimensionless parameters $H$ and $m$. There is only one independent
dimensionless combination of $X$, $k$ and $D$. This fact, combined with the
$1/D$ scaling of $S$, predicted by Eqs. (6) and (11), leads to the scaling
form

$S(X)=\alpha(H,m)\,\frac{k^{2H}|X|^{2(1-2H+Hm)}}{D}\,.$ (22)

Equation (22) establishes the form of the distribution tails
$\mathcal{P}(|X|\to\infty)$ as a function of $|X|$ which is manifestly non-
Boltzmann for all $H\neq 1/2$. At this stage, the important dimensionless
factor $\alpha(H,m)$ is yet unknown except for $H=1/2$, where $S(X)$ describes
the Boltzmann distribution

$S_{H=1/2}(X)=\frac{k|X|^{m}}{mD}\,,$ (23)

so that $\alpha(1/2,m)=1/m$.

Equation (22) also enables us to establish the confinement condition which
makes the steady state possible. Indeed, for the subdiffusion $0<H<1/2$ the
power of $X$ in Eq. (22) is positive for any $m>0$, automatically implying
confinement. For the superdiffusion $1/2<H<1$ the power of $X$ becomes
negative if

$m<m_{0}(H)\equiv\frac{2H-1}{H}\,.$ (24)

In this case the potential $V(x)$ is non-confining, and a steady state is
impossible. Therefore, the potential is confining, and the steady state can be
reached, for $H<1/2$ (for all $m>0$) or for $H>1/2$ and $m>m_{0}(H)$. These
conditions coincide with those presented in Ref. Metzler2022 , where they were
conjectured by using an analogy with the Lévy flights, and then verified in
numerical simulations. The criterion (24) is most stringent when $H\to 1$.

In the rest of this section we determine $\alpha(H,m)$ exactly for $m=2$, and
find the leading-order asymptotics of $\alpha(H,m)$ in the limit of $H\to 0$
and $H\to 1$. Before that, however, let us briefly review the simple case of
white-noise driving, $H=1/2$.

###  II.4 $H=1/2$

This case corresponds, for any $m>0$, to thermal equilibrium, and there is no
need to address the dynamics in order to find the exact (Boltzmann)
distribution. It is still interesting, however, to determine the large-$X$
optimal paths of the particle for different $m$. This is easily done, because
the activation path in equilibrium coincides with the time-reversed relaxation
path. The latter can be found, for any $m>0$, by solving the deterministic
equation $\dot{x}=-k|x|^{m-1}\,\text{sgn}(x)$ with the initial condition
$x(t=0)=X$. As a result, the whole optimal path is

$x(t)=\left[X^{2-m}+k(m-2)|t|\right]^{-\frac{1}{m-2}}\,,$ (25)

where we assumed that $X>0$. For $m\geq 2$ the optimal paths live on the whole
line $|t|<\infty$. Note that as $m>2$ increases the large-$|t|$ power-law
tails become very long.

For $m<2$ the optimal paths have a compact support, $|t|<T<\infty$. For
example, for $m=1$ one obtains

$x(t)=\begin{cases}X-k|t|\,,&|t|<\frac{X}{k}\,,\\\
0\,,&|t|>\frac{X}{k}\,.\end{cases}$ (26)

For a quadratic potential $m=2$, and only in this case, the tails are
exponential. In the next subsection we will consider this important particular
case for all permissible $m$.

###  II.5 Quadratic potential

####  II.5.1 $m=2$: OFM solution

For $m=2$ the Langevin equation (1) becomes linear,

$\dot{x}=-kx+\xi(t)\,,$ (27)

and it describes the fractional Ornstein-Uhlenbeck (fOU) process Cheredito2003
; Kaarakka2015 , a useful generalization of the celebrated Ornstein-Uhlenbeck
process OU  (the latter arises in the particular case of $H=1/2$).

As this linear equation is driven by a Gaussian (albeit non-white) noise, the
steady-state distribution $\mathcal{P}$ is also Gaussian, and it has been
known for some time. Here, once the factor $\alpha(H,2)$ is determined, the
Gaussian tail

$S_{m=2}(X)=\alpha(H,2)\,\frac{k^{2H}X^{2}}{D}\,,$ (28)

predicted by the OFM scaling relation (22) actually gives the whole exact
distribution (upon normalization to unity). Although this factor is known from
the exact solution, see Ref. Metzler2021 , it is instructive to rederive it
within the framework of the OFM. The OFM also provides the (observable and
quite informative) optimal path of the system, not easily accessible
otherwise. In addition, as we see below, it provides the exact autocovariance
of the stationary fOU process.

For $m=2$, Eq. (20) for the optimal path $x(t)$ simplifies to

$\int_{-\infty}^{\infty}dt^{\prime}C(t-t^{\prime})\left[\ddot{x}(t^{\prime})-k^{2}x(t^{\prime})\right]=-\lambda\delta(t)\,.$
(29)

This linear equation can be easily solved by Fourier transform. One obtains an
algebraic equation,

$(\omega^{2}+k^{2})C_{\omega}x_{\omega}=\lambda\,,$ (30)

where $x_{\omega}=\int_{-\infty}^{\infty}dt\,x(t)e^{i\omega t}$. Equation (30)
yields

$x_{\omega}=\frac{\lambda}{(\omega^{2}+k^{2})C_{\omega}}=\frac{\lambda
c_{\omega}}{(\omega^{2}+k^{2})}=\frac{2D\lambda\sin(\pi
H)\Gamma(2H+1)|\omega|^{1-2H}}{\omega^{2}+k^{2}}\,.$ (31)

Now we perform the inverse transform,
$x(t)=(2\pi)^{-1}\int_{-\infty}^{\infty}d\omega\,x_{\omega}e^{-i\omega t}$.
Demanding that $x(t=0)=X$, we obtain

$\lambda=\lambda(X)=\frac{k^{2H}X}{D\Gamma(2H+1)}\,,$ (32)

and the optimal path is given, for all $0<H<1$, by

$\frac{x(t)}{X}=\cosh(kt)-\frac{(k|t|)^{2H}}{\Gamma(2H+1)}\,_{1}F_{2}\left(1;H+\frac{1}{2},H+1;\frac{k^{2}t^{2}}{4}\right)\,,$
(33)

where ${}_{1}F_{2}$ is the hypergeometric function Wolfram . Remarkably, Eq.
(33) coincides, up to a constant factor, with the known autocovariance
$\kappa(t)$ of the _stationary_ fOU process Cheredito2003 . We explain this
coincidence in subsubsection II.5.3.

Using Eqs. (21) and (32), we determine the action:

$S_{m=2}(H)=\int_{0}^{X}dX\,\lambda(X)=\frac{k^{2H}X^{2}}{2D\Gamma(2H+1)}\,.$
(34)

Comparing Eqs. (28) and (34), we obtain

$\alpha(H,2)=\frac{1}{2\Gamma(2H+1)}\,.$ (35)

Interestingly, $\alpha(H,2)$ is a non-monotonic function of $H$, see Fig. 1:
it has a maximum $0.5645\dots$ at $H=0.2308\dots$, the only root of the
digamma function $\psi^{(0)}(2H+1)$ for positive $H$.

Equation (34) yields the variance of $X$,

$\text{Var}_{X,m=2}=\frac{D\Gamma(2H+1)}{k^{2H}}\,,$ (36)

which perfectly agrees with the exact expression presented in Ref. Metzler2021
, see also the footnote half .

It is useful for the following to rederive these results in the frequency
domain. Using Eq. (15), we obtain

$\mathcal{S}[x(t)]=\frac{1}{4\pi}\int_{-\infty}^{\infty}d\omega\,C_{\omega}(\omega^{2}+k^{2})x_{\omega}\,x_{\omega}^{*}\,.$
(37)

The modified action (19) becomes

$\mathcal{S}_{\lambda}[x(t)]=\frac{1}{4\pi}\int_{-\infty}^{\infty}d\omega\,\left[C_{\omega}(\omega^{2}+k^{2})x_{\omega}\,x_{\omega}^{*}-2\lambda
x_{\omega}\right]\,.$ (38)

Its minimization yields Eq. (30) and then Eq. (31). Substituting Eq. (31) into
Eq. (37), we obtain the action

$\mathcal{S}[x(t)]=\frac{\lambda}{4\pi}\int_{-\infty}^{\infty}d\omega\,x_{\omega}\,.$
(39)

Hence

$\mathcal{S}=\frac{1}{2}\lambda X(\lambda)=\frac{1}{2}\lambda(X)X\,,$ (40)

where
$X(\lambda)=(2\pi)^{-1}\int_{-\infty}^{\infty}d\omega\,x_{\omega}=D\Gamma(2H+1)k^{-2H}$.
Using these expressions, we reproduce Eqs. (32)-(35). The most important for a
further use are Eqs. (31) and (39).

Figure 1: Shown is the coefficient $\alpha(H,m)$ in Eq. (22) for $m=2$ (orange
curve) and $m=5/2$ (blue points). The orange curve corresponds to Eq. (35).
The blue points for $H=0$,1/2 and 1 are plotted in accordance to Eq. (54),
whereas three points in the interval $0<H<1/2$ are plotted using numerical
solutions for the optimal paths, described in Sec. III.2, Eq. (74).

####  II.5.2 $m=2$: asymptotics of the optimal paths

Now let us return to the optimal paths (33). Their asymptotic at $t\to 0$,

$\frac{x(t)}{X}\Bigr{|}_{k|t|\ll 1}\simeq
1-\frac{(k|t|)^{2H}}{\Gamma(2H+1)}\,,$ (41)

exhibits a cusp singularity for $H<1/2$. At $H=1/2$ this singularity becomes a
corner singularity, whereas for $H>1/2$ the optimal path is already once
differentiable at $t=0$.

No less instructive is the power-law asymptotic behavior at large $k|t|$,

$\frac{x(t)}{X}\Bigr{|}_{k|t|\gg
1}\simeq\frac{(k|t|)^{2H-2}}{\Gamma(2H-1)}\,,$ (42)

which holds at any $H\neq 1/2$. As one can see, the $k|t|\to\infty$ tails
become very fat as $H$ approaches $1$. Remarkably, in order to reach the point
$x=X>0$ along the least-action path, a subdiffusive particle starting at $x=0$
first moves to the left, and only then to the right.

At $H=1/2$ the optimal path (33) has a simple form $x(t)=X\,e^{-k|t|}$. The
same result follows, in the limit of $m\to 2$, from Eq. (25). As we have
already mentioned, in this case the activation trajectory $0\to X$ coincides
with the time-reversed relaxation trajectory $X\to 0$.

Figure 2: Left panel: The optimal paths $x(t)$ for the quadratic potential
$V(x)=x^{2}/2$ for several values of $H$, as described by Eq. (33). Right
panel: the optimal realization of the fGn, $p(t)\equiv\dot{x}(t)+x(t)$ for the
same set of $H$.

Several examples of the optimal paths $x(t)$ for and $k=1$ and different $H$
are shown in the upper panel of Fig. 2. The lower panel shows, for the same
values of $H$, the optimal realizations of the fGn $\xi(t)$ itself, which are
described by the expression $p(t)\equiv\dot{x}(t)+kx(t)$, see Eq. (1). Notice
that the path $X\to 0$, which brings the particle back to $x=0$, involves a
nonzero noise. For $H$ not close to $1/2$, this noise is quite significant.
(To remind the reader, for $H=1/2$ the “relaxation path” is zero-noise, that
is deterministic.) Overall, like in other systems, the knowledge of the
optimal paths provides a valuable intuition about the physical mechanism
behind the large deviations in question.

The reader may have noticed that the non-equilibrium optimal paths (33) still
exhibit a perfect mirror symmetry in time, $x(-t)=x(t)$, for all values of
$H$. This is just a consequence of the time-reversibility of the fGn and the
linearity of the Langevin equation (1) for $m=2$.

####  II.5.3 $m=2$: autocovariance of the stationary fOU process

The optimal paths, dominating the single-point statistics for permissible $m$
and all $0<H<1$, are defined on the whole time axis, $|x|<\infty$, and they
can be also found by assuming that the process $x(t)$ is _stationary_. In the
particular case of $m=2$, this stationary process – the fOU process – is also
Gaussian. As a result, its one-point statistics can be found by directly
minimizing the Gaussian action of the process conditional on the value $X$.
The corresponding Gaussian action is

$S[x(t)]=\frac{1}{2}\int_{-\infty}^{\infty}dt_{1}\int_{-\infty}^{\infty}dt_{2}K(t_{1}-t_{2})x(t_{1})x(t_{2})\,,$
(43)

where $K(\tau)$ is the inverse kernel for the autocovariance $\kappa(\tau)$ of
the stationary fOU process. Recasting the condition $x(t=0)=X$ as an integral
constraint and introducing a Lagrange multiplier $\lambda$ in the same way as
in Eq. (19), we can reduce the minimization problem to solving a simple
integral equation:

$\int_{-\infty}^{\infty}dt_{1}K(t-t_{1})\kappa(t_{1})=\lambda\delta(t)\,.$
(44)

Its solution is

$x(t)=\lambda\kappa(t)\,.$ (45)

That is, the optimal path of the one-point statistics of the stationary fOU
process coincides, up to a constant factor, with its autocovariance. An
explicit formula for the latter is obtained by dividing $x(t)$ from Eq. (33)
by $\lambda$ given by Eq. (32), and the result coincides with the known
expression Cheredito2003 .

In fact, this important connection between the optimal path of the one-point
statistics and the autocovariance is more general, as it holds for _all_
stationary Gaussian processes Meerson2022 ; MO2022 .

###  II.6 Asymptotic limits of $H\to 0$ and $H\to 1$

####  II.6.1 $H\to 0$

Sending $H$ to $0$ in Eq. (22), we obtain $S(X,H=0)=\alpha(0,m)X^{2}/D$.
Remarkably, the deterministic relaxation rate $k$ drops out. This implies
that, in this limit, the action is entirely determined by the noise and is
independent of the external potential $V(x)$. In particular, the coefficient
$\alpha(0,m)$ must be independent of $m$. Plugging $H=0$ into Eq. (35), we
obtain $\alpha(0,2)=1/2$. Therefore, $\alpha(0,m)=1/2$ for all permissible
$m$, and we arrive at a universal (potential-independent) $H\to 0$ asymptotic
of the action:

$S_{H\to 0}(X)=\frac{X^{2}}{2D}\,.$ (46)

In the frequency domain, this asymptotic is dominated by the high frequencies
of the fGn, $|\omega|\gg k$. For $m=2$ this fact becomes evident. Indeed, the
calculation of the action by using Eq. (40) is reduced, at any $H$, to
evaluating the integral

$\int_{-\infty}^{\infty}d\omega\,\frac{|\omega|^{1-2H}}{\omega^{2}+k^{2}}\,.$
(47)

At $H\to 0$, this integral is dominated by large $\omega$, up to $|\omega|\sim
k\exp[1/(2H)]$.

####  II.6.2 $H\to 1$

When $H$ approaches $1$ in Eq. (22), we obtain $S(X,H\to
1)=\alpha(1,m)f(X)^{2}/D$. As we will see, this will lead us to the conclusion
that $\alpha(1,m)$ does not depend on $m$ as well. Let us consider the
particular case of $m=2$. When $H\to 1$, Eq. (35) becomes

$S_{H\to 1}(X)=\frac{f(X)^{2}}{4D}=\frac{(kX)^{2}}{4D}\,.$ (48)

In this limit the integral in Eq. (47) is mostly contributed to by the low
frequencies, $|\omega|\ll k$. Here one can neglect the term $\dot{x}$ in the
Langevin equation (1). Therefore, to zero order, we have a balance
$-kx(t)\simeq\bar{\xi}_{T}(t)$ at $T\gg k^{-1}$, where

$\bar{\xi}_{T}(t)=\frac{1}{2T}\int_{-T}^{T}\xi(t+\tau)\,d\tau\,,$ (49)

is the low-frequency part of the fGn $\xi$. Using Eq. (4) for
$\langle\xi(t)\xi(t^{\prime})\rangle$, we obtain at a given $T$ and $H$
tending to 1 that $\langle\bar{\xi}_{T}^{2}\rangle\bigr{|}_{H\to 1}=2D$.
Therefore, Eq. (48) describes the Gaussian distribution
$\mathcal{P}(\bar{\xi}_{T})$ of $\bar{\xi}_{T}$, that is,

$\mathcal{P}_{H\to 1}(\bar{\xi}_{T})\simeq\frac{1}{\sqrt{4\pi
D}}\exp\left(-\frac{\bar{\xi}_{T}^{2}}{4D}\right)\,.$ (50)

Similarly, the approximate force balance $-f(x)\simeq\xi_{T}(t)$ holds, in the
limit of $H\to 1$, for all monotonically growing forces $f(x)$. Using the
force balance and Eq. (50), we obtain the PDF of $X$ we are after:

$\mathcal{P}_{H\to 1}[f(X)]\simeq\frac{1}{\sqrt{4\pi
D}}\exp\left[-\frac{f(X)^{2}}{4D}\right]\,.$ (51)

In particular, for the power-law potentials (2) with $m>1$, Eq. (51) becomes

$\mathcal{P}_{H\to 1}(X)\simeq\frac{k(m-1)|X|^{m-2}}{\sqrt{4\pi
D}}\exp\left(-\frac{k^{2}|X|^{2m-2}}{4D}\right)\,,$ (52)

As one can see, in this limit we have been able to account for the pre-
exponential factor of the distribution as well, and therefore to describe the
whole distribution (52). Strikingly, for $m>2$ the distribution (52) is
bimodal: it vanishes at $X=0$. The bimodality of $\mathcal{P}(X)$ at $m>2$ in
the particular limit of $H\to 1$ agrees with the findings of Ref. Metzler2021
who observed bimodality numerically for all $H>1/2$.

The action, corresponding to Eq. (52), is

$S_{H\to 1}(X)\simeq\frac{k^{2}|X|^{2m-2}}{4D}\,.$ (53)

Comparing this expression with Eq. (22), we conclude that $\alpha(1,m)=1/4$
independently of $m$.

Altogether, we know the values of the factor $\alpha(H,m)$ in Eq. (22) in
three limits:

$\alpha(H,m)=\begin{cases}\frac{1}{2}\,,&H\to 0\,,\\\
\frac{1}{m}\,,&H=1/2\,,\\\ \frac{1}{4}\,,&H\to 1\,,\end{cases}$ (54)

where the first and third lines do not depend on $m$.

###  II.7 Path irreversibility

One important signature of nonequilibrium is path irreversibility: a disparity
between the back and forth probabilities of the system moving along a closed
path in the parameter space. We can use the expression (6) to calculate the
difference $S[x(t)]-S[x(-t)]$:

$S[x(t)]-S[x(-t)]=-2\int\limits_{-\infty}^{\infty}\int\limits_{-\infty}^{\infty}C(t-t^{\prime})\dot{x}(t)f(x(t^{\prime}))\,.$
(55)

This expression vanishes for $H=1/2$, where closed paths, obeying (18), are
reversible for any $m>0$. For $H\neq 1/2$ they are irreversible except for
$m=2$ when the Langevin equation (1) is linear.

As an illustrative example of irreversibility, let us consider an asymmetric
closed path $0\to X\to 0$:

$x(t)=\begin{cases}X\,e^{t}\,,&t<0\,,\\\ X\,e^{-2t}\,,&t>0\\\ \end{cases}$
(56)

in a quartic potential, $m=4$. We evaluated the action (6) for several values
of $0<H<1/2$ by performing the integrations over time from $t=-\infty$ to
$\infty$ and backward. Let us denote the resulting forward and backward
actions by $S_{+}(H)$ and $S_{-}(H)$, respectively. Figure 3 shows the ratio
$S_{+}/S_{-}$ as a function of $H$, and irreversibility is evident. As $H$
approaches $1/2$, the ratio $S_{+}/S_{-}$ approaches $1$ as to be expected.

Figure 3: The ratio of the forward and backward actions (6) vs. $H$ for the
closed path (56). The parameters are $D=1/2$ and $k=X=1$.

##  III Regularization of Eq. (20) and numerical optimal paths

###  III.1 Regularization

Throughout this section we confine ourselves to the subdiffusion case,
$0<H<1/2$. Inspecting the $t\to 0$ asymptotic behavior of the exact optimal
path $x(t)$ for $m=2$, Eq. (41), and the relationship (32) between $\lambda$
and $X$, one can notice that they are independent of $m$. Therefore, the
character of the cusp singularity of $x(t)$ at $t\to 0$,

$x(t)\bigr{|}_{t\to 0}=X-\lambda D|t|^{2H}+\dots\,,$ (57)

remains the same for all $m>0$. This reflects the fact that, for $H<1/2$, the
second term in the asymptotic (57) is universal and independent of $m$
(actually, of the confining force $f(x)$ in general). This observation plays a
critical role in our approach to a numerical solution of Eq. (20).

The most singular term in the left hand side of Eq. (20) comes from the second
derivative of the second term in Eq. (57). A further observation is that the
first two terms of the asymptotic (57) provide an exact solution of a
_truncated_ Eq. (20) which includes only the most singular term on the left
hand side – the term originating from $\ddot{x}$) – and the delta-function
term on the right hand side. All the other contributions, which come from the
rest of the left hand side terms and from the higher-order terms in the
expansion (57), are much less singular at $t\to 0$.

As one check, the second term in the expansion (57) is so singular that it
makes it impossible to treat the integral in Eq. (20) as an improper integral
at $t^{\prime}=0$ of a function that is regular elsewhere. As a result, Eq.
(20) cannot be used for numerical analysis in its present form. The main goal
of this section is to obtain a regularized form of Eq. (20), and we achieve
this goal via an analytical treatment of the most singular part of $x(t)$
coming from the 2nd term in the expansion (57). To this end we present $x(t)$,
for all $t$, in the following form:

$x(t)=\tilde{x}(t)-\lambda D|t|^{2H}\,,$ (58)

where $\tilde{x}$ is a much smoother function at $t=0$ than $x(t)$, and try to
derive a regular equation for $\tilde{x}$.

For brevity, we rewrite Eq. (20) in the following form:

$C*\ddot{x}+f^{\prime}(x)\cdot
C*\dot{x}-C*f^{\prime}(x)\cdot\dot{x}-f^{\prime}(x)\cdot
C*f(x)+\lambda\delta(t)=0\,.$ (59)

Here all functions are functions of $t$, $C(t)$ is defined by Eq. (7), the
symbol $F\cdot$ denotes ordinary multiplication by the function $F$, and the
asterisk denotes the convolution of two functions of $t$. Rearranging the
terms in this equation and using the identity $c(t)*C(t)=\delta(t)$, we can
rewrite Eq. (59) in an equivalent form,

$C*\left[\ddot{x}-f(x)f^{\prime}(x)+\lambda
c(t)\right]=\left[f^{\prime}(x)\cdot,\,C\,*\right]\left(f-\dot{x}\right)\,,$
(60)

where the commutator $\left[f^{\prime}(x)\cdot,\,\hat{C}\,*\right]$ is defined
as

$\left[f^{\prime}(x)\cdot,\,C\,*\right]\,F=f^{\prime}(x)\cdot
C*F-C*f^{\prime}(x)\cdot F\,.$ (61)

[Notice that the commutator in Eq. (60) vanishes at $m=2$ for any $H$, or at
$H=1/2$ for any $f(x)$.] Making the convolution of the both sides of Eq. (60)
with the kernel $c(t)$, we obtain

$\ddot{x}-f(x)f^{\prime}(x)+\lambda
c(t)=c*\left[f^{\prime}(x)\cdot,\,C\,*\right]\left(f(x)-\dot{x}\right)\,.$
(62)

Using Eqs. (4) and (7) for $c(t)$ and $C(t)$, respectively, we can rewrite Eq.
(62) as

$\frac{d^{2}}{dt^{2}}\left\\{x(t)+D\lambda|t|^{2H}+Z[x](t)\right\\}=f(x(t))f^{\prime}(x(t))\,,$
(63)

where

$Z[x](t)=C_{0}\,|t|^{2H}\left[f^{\prime}(x(t))\cdot,\,|t|^{-2H}\,*\right]\left[\dot{x}(t)-f(x(t))\right]$
(64)

and $C_{0}=\cot(\pi H)/(4\pi H)$. Going over from the dependent variable
$x(t)$ to the new dependent variable

$\bar{x}(t)=x(t)+D\lambda|t|^{2H}+Z[x](t)\,,$ (65)

we arrive at the following equation for $\bar{x}(t)$:

$\frac{d^{2}}{dt^{2}}\bar{x}(t)=f(x)f^{\prime}(x)\,,$ (66)

where the old variable $x(t)$ in the right-hand-side should be expressed
through $\bar{x}(t)$:

$x(t)=\bar{x}(t)-\lambda D|t|^{2H}-Z[x](t)=\tilde{x}(t)-Z[x](t)\,.$ (67)

Notice that $Z[x](t)$ vanishes identically at $m=2$ for any $H$ and at $H=1/2$
for any $m$. In these cases Eq. (66) turns into a second-order ODE (linear or
nonlinear, respectively) which can be exactly solved for $\bar{x}(t)$. As to
be expected, the resulting $x(t)=\bar{x}(t)-\lambda D|t|^{2H}$ perfectly
agrees with the corresponding solutions in these cases that we presented
above.

As we argue below, in the general case $Z[x](t)$ is a continuous function of
$t$ at least on the solutions of Eq. (63) [or (66)]. Let us assume for a
moment that this continuous function $Z[x](t)$ is known. Then Eq. (66) becomes
an ordinary differential equation (ODE), and it has solutions $\bar{x}(t)$
which are twice differentiable. Some examples of iterative numerical solutions
of the problem, defined by Eqs. (66), (67) and (18), will be presented in Sec.
III.2.

Rigorously speaking, the integral corresponding to the second convolution,
$|t|^{2H}*\dots$, in Eq. (64) for $Z[x](t)$, diverges at
$t^{\prime}\to\pm\infty$. However, using the freedom of adding to $Z[x](t)$
any linear function of $t$ without affecting the ODE, we can “improve” the
definition of $Z$ so that it becomes bounded for all finite times, with the
same solution $x(t)$. Such an improved definition is considered in Appendix A.
This redefinition, however, becomes unnecessary in the numerical solution of
Eq. (63) due to the finiteness of the numerical time interval, see Appendix B.

Now we present our argument that $Z[x](t)$ is a continuous function of $t$ at
least on the solutions of Eq. (63). Let us recast the definition of $Z[x](t)$
in Eq. (64) in a shorter form,

$Z[x](t)\propto|t|^{2H}*Y[x](t)\,,$ (68)

where

$Y[x](t)=\left[f^{\prime}(x(t))\cdot,\,|t|^{-2H}\,*\right]\left[\dot{x}(t)-f(x(t))\right]\,.$
(69)

The most singular term at $t=0$ in Eq. (69) comes from the derivative of the
second term in the asymptotic (57). Keeping only this singular term in Eq.
(69), we see that the most singular part of the expression for $Y$ at $t\to 0$
is the following:

$Y[x](t)_{t\to
0}\propto\mbox{p.v.}\int\frac{|t^{\prime}|^{2H}\left(|t^{\prime}|^{2H}-|t|^{2H}+\dots\right)\,dt^{\prime}}{t^{\prime}\,|t-t^{\prime}|^{2H}}+\dots\,,$
(70)

where we have kept only the most singular terms in the vicinity of
$t^{\prime}=0$ when $t\to 0$. Now we can see that all the singularities of the
integrand are integrable at $t^{\prime}=0$ and $t^{\prime}=t$. As a result, we
have:

$Y[x](t)_{t\to 0}=A+B|t|^{2H}+\dots\,,$ (71)

with some constant $A$ and $B$. It is evident now that the leading singularity
of $Y$ at $t\to 0$ is the same as the singularity of $x(t)$, and that the
function $Y[x](t)$ is continuous at $t=0$. Therefore, by virtue of Eq. (68),
$Z[x](t)$ is also a continuous function of $t$ . Our numerical solutions of
this problem, presented in Sec. III.2, confirm these (not entirely rigorous)
arguments.

###  III.2 Numerical solution

In the general case of a nonzero $Z[x](t)$, Eq. (63) is a nonlinear integro-
differential equation which cannot be solved analytically. Here we propose a
numerical iteration procedure which allows one to get sufficiently accurate
numerical solutions of Eq. (63), and hence Eq. (20), for $H<1/2$ in spite of
the expected singularity of the solution at $t=0$.

The numerical iteration procedure treats Eq. (66) as a second-order ODE for
$\bar{x}(t)$, defined by Eq. (65), with the boundary conditions

$x(t=-T)=x(t=T)=0\,,$ (72)

where $T$ is sufficiently large but finite. To obey the boundary conditions we
used a shooting method which involves a numerical solution of the Cauchy
problem with the fixed condition $x(t=-T)=0$ but an adjustable condition for
$\dot{\bar{x}}(t=-T)$. The function $Z[x](t)$ is calculated by using
$\tilde{x}(t)$ from the previous iteration. We found out that this iteration
procedure converges quite fast to its limiting solution. Some additional
details of the iteration procedure are presented in Appendix B.

We used this iteration procedure to compute the optimal paths and the action
for $m=5/2$ and three values of $H$: $1/8$, $1/4$ and $3/8$. We set $k=1$ and
$\lambda D=1/2$ without limiting generality. Having found numerically
$x(t=0)\equiv X$, we computed the dimensionless coefficient $\alpha(H,m)$ for
the chosen values of $H$ and $m$ by substituting Eq. (22) for the action $S$
into the shortcut relation (21):

$\frac{dS}{dX}\biggr{|}_{X=x(0)}=\frac{1}{2D}\,.$ (73)

This gives

$\alpha(H,m)=\frac{X^{4H-1-2Hm}}{4(1-2H+Hm)}\,.$ (74)

The coefficients $\alpha(H,m)$ in Eq. (22), computed in this way, are
presented in Fig. 1.

Figure 4: Numerical optimal path $x(t)$ for $m=5/2$ and $H=1/4$. The
parameters $k$ amd $\lambda D$ are $k=1$ and $\lambda D=1/2$.

As an example, Fig. 4 shows the numerically found optimal path $x(t)$ for
$H=1/4$. It has a cusp singularity at $t=0$ which is of the same type as the
one for $m=2$, see Fig. 2. This supports our argument presented in Sec. III.1.
As explained in Sec. II.7, the optimal path is expected to show
irreversibility, that is to be asymmetric with respect to the reflection
$t\leftrightarrow-t$. The asymmetry is not so evident in Fig. 4, but it
becomes evident in the plots of the functions $(x(t)-x(-t))/2$ and
$(x(t)+x(-t))/2$, see Fig. 5.

Figure 5: Irreversibility of the optimal path shown in Fig. 4. Plotted are the
functions $(x(t)-x(-t))/2$ and $(x(t)+x(-t))/2$.

Figure 6 shows $Z(t)$ computed by using the numerically found $x(t)$ shown in
Fig. 4. The function $Z(t)$ is manifestly continuous, in line with our
argumentation in Sec. III.1. We also checked that the continuity of $Z[x](t)$
is preserved at each step of the iteration procedure.

Figure 6: The function $Z[x](t)$ evaluated on the numerical optimal path
$x(t)$ for $H=1/4$, shown in Fig. 4.

##  IV Summary and Discussion

As we have demonstrated, the OFM provides a valuable insight into the (in
general, nonequilibrium) steady state of fBM in scale-invariant confining
potentials. In the particular case of a quadratic potential, which corresponds
to the fOU process, the OFM supplies exact results not only for the steady-
state probability distribution of the particle position, but also for the
autocovariance of the stationary process. Additional analytical results can be
obtained in the limits of $H\to 0$ and $H\to 1$. Importantly, these asymptotic
results apply to a broad class of confining potentials, not necessarily scale-
invariant. For non-quadratic potentials, the OFM problem for the position
distribution can be solved numerically. In all cases, the OFM predicts the
optimal paths of the conditional process, which highlight the physics of the
process.

It would be very interesting to test our predictions for the tails of
$\mathcal{P}(X)$, and for the corresponding optimal paths, in large-deviation
simulations of the fractional Langevin equation (1). Very recently, large-
deviation simulations of a free (unconfined) fBm, have been successfully
accomplished in the context of the first-passage area distribution and its
optimal paths HM2024 .

A natural, and very interesting, next step would be applying the OFM to the
problem of activated escape of fBM from a potential well due to external fGn
Metzler2010 ; Wang2017 . Here the OFM should provide the optimal path to
escape and an asymptotically exact estimate for the logarithm of the mean
escape time.

## Acknowledgments

We are very grateful to Alexander K. Hartmann and Alexander Valov for useful
discussions. B. M. was supported by the Israel Science Foundation (Grant No.
1499/20).

##  Appendix A Redefinition of $Z[x](t)$

Let us recast the definition (64) of $Z[x](t)$ in the following form:

$Z[x](t)=C_{0}\,\int\limits_{-\infty}^{\infty}dt^{\prime}\,|t-t^{\prime}|^{2H}Y[x](t^{\prime})\,,$
(75)

where

$Y[x](t)=\left[f^{\prime}(x(t))\cdot,\,|t|^{-2H}\,*\right]\left[\dot{x}(t)-f(x(t))\right]\,.$
(76)

The integral in Eq. (75) over $t^{\prime}$ diverges at
$t^{\prime}\to\pm\infty$ for any $t$. To cure this divergence, we can improve
the definition (75) by introducing

$Z_{T}[x](t)=C_{0}\,\int\limits_{-T}^{T}dt^{\prime}\,|t-t^{\prime}|^{2H}Y[x](t^{\prime})$
(77)

and by defining $Z[x](t)$ as follows:

$Z[x](t)=\lim\limits_{T\to\infty}\biggl{\\{}Z_{T}[x](t)-Z_{T}[x](0)-t\frac{Z_{T}[x](T_{1})-Z_{T}[x](-T_{1})}{2T_{1}}\biggr{\\}}\,,$
(78)

where $T_{1}$ is an arbitrary positive constant of order $1$.

##  Appendix B Details of the iteration procedure

We implemented the iteration procedure, outlined in Sec. III.2 for the
function $\tilde{x}(t)$, which is defined as

$\tilde{x}(t)=x(t)+\lambda D|t|^{2H}$ (79)

and presumably has a continuous first derivative. We consider this function
only on the interval $|t|<T$, where $0<T<\infty$. As a result, all the
integrals, entering the definition of the convolutions in Eq. (64), are
computed on the same time interval. For this reason there is no need for a
redefinition of $Z[x](t)$, caused by the formal divergence of the second
convolution in the expression (64) at $T=\infty$.

During the iterations of $\tilde{x}_{n}$,

$\tilde{x}_{n}\to\tilde{x}_{n+1}\quad(n=0,1,2,\dots)\,,$ (80)

the terms $x_{n}(t)$, which appear below, are evaluated by using Eq. (79) for
the known $\tilde{x}_{n}(t)$. The derivative of $x_{n}(t)$ is evaluated by
computing the numerical derivative of $\tilde{x}_{n}$ and the explicit
analytical derivative of the second term in Eq. (79). Each step of the
iterations consists of several sub-steps which are listed below. The
expression

$\displaystyle Y[\tilde{x}_{n}](t)$ $\displaystyle=$
$\displaystyle\int\limits_{-T}^{T}dt^{\prime}\,|t-t^{\prime}|^{-2H}\left[f^{\prime}(x_{n}(t))-f^{\prime}(x_{n}(t^{\prime}))\right]\left[\dot{x}_{n}(t^{\prime})-f(x_{n}(t^{\prime}))\right]$
(81) $\displaystyle=$
$\displaystyle\int\limits_{-T}^{T}dt^{\prime}\,|t-t^{\prime}|^{-2H}\left[f^{\prime}(x_{n}(t))-f^{\prime}(x_{n}(t^{\prime}))\right]\left[\dot{\tilde{x}}_{n}(t^{\prime})-2\lambda
DH|t^{\prime}|^{2H-1}\,\mbox{sign}\,t^{\prime}-f(x_{n}(t^{\prime}))\right]$

evaluates
$\left[f^{\prime}(x(t))\cdot,\,|t|^{-2H}\,*\right]\left[\dot{x}(t)-f(x(t))\right]$
in Eq. (64). Further,

$Z[\tilde{x}_{n}](t)=C_{0}\,\int\limits_{-T}^{T}dt^{\prime}\,|t-t^{\prime}|^{2H}Y[\tilde{x}_{n}](t^{\prime})\,.$
(82)

We introduce

$X[\tilde{x}_{n}](t)=\tilde{X}[\tilde{x}_{n}](t)-\lambda
D|t|^{2H}\,;\quad\tilde{X}[\tilde{x}_{n}](t)=\bar{X}[\tilde{x}_{n}](t)-Z[\tilde{x}_{n}](t)\,.$
(83)

Now we solve numerically the ODE

$\ddot{\bar{X}}-f(X)f^{\prime}(X)=0\qquad\left(X(t)_{|t|=T}=0\right)\,.$ (84)

Then we calculate $\tilde{x}_{n}$ of the next step of the iterations as
follows:

$\tilde{x}_{n+1}(t)=(1-\alpha_{n})\,\tilde{x}_{n}(t)+\alpha_{n}\tilde{X}[\tilde{x}_{n}](t)\,.$
(85)

$\alpha_{n}=1$ would correspond to the straightforward iterations briefly
considered in Sec. III.1: $\tilde{x}_{n+1}(t)=\tilde{X}[\tilde{x}_{n}](t)$.
However, to avoid a possible instability of the iteration process we used the
more conservative weights

$\alpha_{0}=1\quad\mbox{and}\quad\alpha_{n\geq 1}=\alpha=0.5\,.$

We started the iteration process by settting $x_{0}(t)=0$, so that
$Z[x_{0}](t)=0$, and $\bar{x}_{1}=\tilde{x}_{1}=\tilde{X}=\bar{T}$ is the
solution of the ODE
$\ddot{\bar{X}}-f(X)f^{\prime}(X)\bigr{|}_{X=\bar{X}-\lambda D|t|^{2H}}=0$.

We observed, for $H=1/4$ and $m=5/2$, that replacing $T=\infty$ by $T=15$ does
not lead to an error in the value of $x(0)$ that would exceed $\pm 0.5\%$. The
iterations converge relatively fast: about 10-15 iterations were sufficient to
find the limiting optimal path for $T=15$ within the above-mentioned
inaccuracy of less $\pm 0.5\%$ caused by the finite $T=15$.

## References

  * (1) A. N. Kolmogorov, C. R. Dokl. Acad. Sci. URSS 26, 115 (1940).
  * (2) B. B. Mandelbrot and J. W. van Ness, SIAM Review 10, 422 (1968).
  * (3) H. Qian, in “Processes with Long-Range Correlations. Theory and Applications”, edited by G. Rangarajan and M. Ding (Springer, Berlin, 2003), p. 22.
  * (4) T. Guggenberger, A. Chechkin, and R. Metzler, J. Phys. A: Math. Theor. 54, 29 (2021).
  * (5) T. Guggenberger, A. Chechkin, and R. Metzler, New J. Phys. 24, 073006 (2022).
  * (6) B. Meerson, O. Bénichou, and G. Oshanin, Phys. Rev. E 106, L062102 (2022).
  * (7) Y. Liang, W. Wang, R. Metzler, and A. G. Cherstvy Phys. Rev. E 108, L052101 (2023).
  * (8) F. D. Cunden, P. Facchi, and P. Vivo, J. Phys. A: Math. Theor. 49, 135202 (2016).
  * (9) G. I. Barenblatt, Scaling, Self-similarity, and Intermediate Asymptotics: Dimensional Analysis and Intermediate Asymptotics (Cambridge Texts in Applied Mathematics, Series Number 14) (Cambridge University Press, Cambridge, UK, 1996).
  * (10) P. Cheredito, H. Kawaguchi, and M. Maejima, Electron. J. Probab. 8, 1 (2003).
  * (11) T. Kaarakka, Fractional Ornstein-Uhlenbeck Processes (Tampere Univ. Technol., Tampere, Finland, 2015).
  * (12) G. E. Uhlenbeck and L. S. Ornstein. Phys. Rev. 36, 823 (1930).
  * (13) Wolfram Research, Inc., Mathematica, Version 14.0.0.0, Champaign, IL (2023).
  * (14) When comparing our formulas with those of Refs. Metzler2021 ; Metzler2022 ; MBO , one should notice that in those papers $D$ was set to be $1/2$.
  * (15) B. Meerson, Phys. Rev. E 105, 034106 (2022); 107, 039902 (2023).
  * (16) B. Meerson and G. Oshanin, Phys. Rev. E 105, 064137 (2022).
  * (17) A. K. Hartmann and B. Meerson, Phys. Rev. E 109, 014146 (2024).
  * (18) O. Yu. Sliusarenko, V. Yu. Gonchar, A. V. Chechkin, I. M. Sokolov, and R. Metzler, Phys. Rev. E 81, 041119 (2010).
  * (19) W. Wang, Z. Yan, and X. Liu, Phys. Lett. A 381, 2324 (2017).


Having discussed the form of the AGPs wavefunction, we now discuss how it can be implemented using FermiNet.
In the original FermiNet architecture, the outputs of the one-electron stream are used to build FermiNet orbitals $\phi^{k\alpha}_{i}(\mathbf{r}^\alpha_j;\{\mathbf{r}^\alpha_{/j}\}; \{\mathbf{r}^{\bar\alpha}\})$.
The full many-fermion wavefunction is a weighted sum of terms, each of which is the product of one up-spin and one down-spin determinant of the FermiNet orbital matrices, as shown in Eq. (<ref>).

To build a pairing orbital, one can make use of the outputs $\mathbf{h}^{L\alpha}_i$ from the last layer $L$ of the one-electron stream.
Instead of using these outputs to build FermiNet orbitals, as in Eq. (<ref>), they can be used to build FermiNet pairing orbitals as follows:

\begin{multline}
\varphi^{k}(\mathbf{r}^\alpha_i, \mathbf{r}^{\bar\alpha}_j;
\{\mathbf{r}^\alpha_{/i}\};\{\mathbf{r}^{\bar\alpha}_{/j}\}) \\
= [\mathbf{w}^{k} \cdot (\mathbf{h}^{L\alpha}_i \odot
\mathbf{h}^{L\bar\alpha}_j)] \chi^{k}(\mathbf{r}^\alpha_i)
\chi^{k}(\mathbf{r}^{\bar\alpha}_j),
\label{eq:ggno2e}
\end{multline}

where $\chi(\mathbf{r})$ are the envelope functions, $\mathbf{w}^{k}$ are vectors and $\odot$ denotes the element-wise product.
This procedure generates a pairing function between particles $\mathbf{r}^\alpha_i$ and $\mathbf{r}^{\bar\alpha}_j$, retaining the permutation invariant property possessed by FermiNet orbitals.
Depending on the number of FermiNet geminals generated, the AGPs can be written as one or a weighted sum of multiple determinants of the FermiNet geminals

\begin{align}
&\Psi^\text{AGPs}_\text{FermiNet}\left(\mathbf{r}^\uparrow_1, \dots, \mathbf{r}^\downarrow_{N^\downarrow} \right) \notag \\
&= \sum_k \omega_k \det[\varphi^{k}(\mathbf{r}^\alpha_i, \mathbf{r}^{\bar\alpha}_j; \{\mathbf{r}^\alpha_{/i}\};\{\mathbf{r}^{\bar\alpha}_{/j}\})] .
\label{eq:agps_ferminet}
\end{align}

This is analogous to a weighted sum of single-determinant AGPs wavefunctions of the type defined in Eq. (<ref>), but the replacement of the two-particle geminals by FermiNet geminals that depend on the positions of all of the particles makes it much more general.

Although using the outputs from the one-electron stream is sufficient to build an AGPs, one can also include the outputs from the two-electron stream:

\begin{align}
\varphi^{k}&(\mathbf{r}^\alpha_i, \mathbf{r}^{\bar\alpha}_j;
\{\mathbf{r}^\alpha_{/i}\};\{\mathbf{r}^{\bar\alpha}_{/j}\}) & \notag \\
= [&\mathbf{w}^{k}_1 \cdot (\mathbf{h}^{L\alpha}_i \odot
\mathbf{h}^{L\bar\alpha}_j) \notag \\
&+ \mathbf{w}^{k}_2 \cdot
(\mathbf{h}^{L\alpha\bar\alpha}_{ij} \odot
\mathbf{h}^{L\bar\alpha\alpha}_{ji})]
\chi^{k}(\mathbf{r}^\alpha_i)\chi^{k}(\mathbf{r}^{\bar\alpha}_j) ,
\label{eq:ggw2e}
\end{align}

where $\mathbf{w}^{k}_1$ and $\mathbf{w}^{k}_2$ are vectors.

Note that Eqs. (<ref>) and (<ref>) are two possible ways of building a pairing function.
There are many others ways and they are all valid as long as the appropriate symmetries are preserved.
An alternative method is given by Xie et al. <cit.>.

The benefit of building AGPs-like wavefunctions using the FermiNet is that the pairing function $\varphi(\mathbf{r}^\alpha_i, \mathbf{r}^{\bar\alpha}_j;\{\mathbf{r}^\alpha_{/i}\};\{\mathbf{r}^{\bar\alpha}_{/j}\})$ now depends not only on $\mathbf{r}_i^{\alpha}$ and $\mathbf{r}_j^{\bar{\alpha}}$ but also on the positions of the other particles in the system.
Correlations between the singlet pair and the other particles can thus be captured.
In a similar way, the original unpaired FermiNet replaced Hartree-Fock-like single-particle orbitals $\phi_i^{\alpha}(\mathbf{r}_j^{\alpha})$ by FermiNet orbitals $ \phi^{\alpha}_{i}(\mathbf{r}^\alpha_j;\{\mathbf{r}^\alpha_{/j}\}; \{\mathbf{r}^{\bar\alpha}\})$, helping to capture correlations between the particle at $\mathbf{r}_j^{\alpha}$ and all other particles.

§.§ Relations between the original FermiNet and the AGPs FermiNet

Next, we clarify the relation between the original FermiNet and the AGPs FermiNet, showing that the AGPs FermiNet is the more general of the two.
A FermiNet geminal with a two-particle stream term is even more general than a FermiNet geminal without, so it is sufficient for this purpose to omit the two-particle stream term.
We also neglect the envelope functions.
Including them circumvents numerical difficulties in finite systems and can speed up the network optimization, but does not affect the generality of the Ansatz.

Let us first define a pairing function in the following way:

\begin{equation}
\varphi(\mathbf{r}^\uparrow_i, \mathbf{r}^\downarrow_j;
\{\mathbf{r}^\uparrow_{/i}\};\{\mathbf{r}^{\downarrow}_{/j}\})
= \sum_{kl} W_{kl} h^{L\uparrow (k)}_i h^{L\downarrow(l)}_j ,
\label{eq:ggij}
\end{equation}

where $h^{L\alpha (k)}_i = [\mathbf{h}^{L\alpha}_i]_k$ are the outputs from the final layer of the one-electron stream for particle $i$ of spin $\alpha$.
As we explain below, Eq. (<ref>) is equivalent to the simpler FermiNet geminal described above:

\begin{equation}
\varphi(\mathbf{r}^\uparrow_i, \mathbf{r}^{\downarrow}_j;
\{\mathbf{r}^\uparrow_{/i}\};\{\mathbf{r}^{\downarrow}_{/j}\})
= \mathbf{w} \cdot (\mathbf{h}^{l\uparrow}_i \odot \mathbf{h}^{l\downarrow}_j).
\label{eq:ggwob}
\end{equation}

We choose to write the pairing function in the form of Eq. (<ref>) only because this makes it easier to relate to FermiNet orbitals.
Since Eqs. (<ref>) and (<ref>) are equivalent, the choice does not affect the conclusions of the argument.
In the rest of this section, for the sake of simplicity, we omit the sets $\{\mathbf{r}_{/i}^{\uparrow}\}$ and $\{\mathbf{r}_{/j}^{\downarrow}\}$ from the arguments of FermiNet pairing functions and orbitals.

To explain the equivalence of Eqs. (<ref>) and (<ref>), it is helpful to represent the matrix $W_{kl}$ as its singular-value decomposition (SVD):

\begin{equation}
W_{kl} = \sum^{n_L}_{\gamma=1} \sigma_\gamma U_{\gamma k} V_{\gamma l},
\label{eq:svd}
\end{equation}

where $U \in \mathbb{R}^{n_L \times n_L}$ and $V  \in \mathbb{R}^{n_L \times n_L}$ are orthogonal matrices and $n_L$ is the size of the vectors $\mathbf{h}_i^{L\alpha}$ output by the final layer $L$ of the one-electron stream.
This is also known as the number of hidden units in layer $L$.
The pairing function in Eq. (<ref>) becomes

\begin{align}
\varphi(\mathbf{r}^\uparrow_i, \mathbf{r}^\downarrow_j)
&= \sum_{kl} W_{kl} h^{L\uparrow (k)}_i h^{L\downarrow(l)}_j \\
&= \sum_\gamma \sigma_\gamma \left(\sum_k U_{\gamma k} h^{L\uparrow (k)}_i\right)
\left(\sum_k V_{\gamma l} h^{L\uparrow (l)}_j\right) \\
&= \sum_\gamma \sigma_\gamma (U \mathbf{h}^{L\uparrow}_i)_\gamma
(V \mathbf{h}^{L\downarrow}_j)_\gamma.
\end{align}

Given the universal approximation theorem <cit.>, and the fact that every layer of the network contains an arbitrary linear transformation, it is reasonable to assume that the functions $\mathbf{h}^{L\alpha}_i$ and $O\mathbf{h}^{L\alpha}_i$, where $O$ is $U$ or $V$, have the same variational freedom and information content.
In other words, we assume that any network capable of representing $\mathbf{h}^{L\alpha}_i$ can also represent $O\mathbf{h}^{L\alpha}_i$, since this is merely a rotation of the vectors $\mathbf{h}^{L\alpha}_i$ in the last layer.
We thus define $\tilde{\mathbf{h}}^{L\uparrow}_i = U \mathbf{h}^{L\uparrow}_i$ and $\tilde{\mathbf{h}}^{L\downarrow}_i = V \mathbf{h}^{L\downarrow}_i$, such that the pairing function becomes

\begin{align}
\varphi(\mathbf{r}^\uparrow_i, \mathbf{r}^\downarrow_j)
&= \sum_\gamma \sigma_\gamma (U \mathbf{h}^{L\uparrow}_i)_\gamma
(V \mathbf{h}^{L\downarrow}_j)_\gamma \notag \\
&= \sum_\gamma \sigma_\gamma (\tilde{\mathbf{h}}^{L\uparrow}_i
\odot \tilde{\mathbf{h}}^{L\downarrow}_j)_{\gamma},
\end{align}

which is equivalent to Eq. (<ref>).

To relate the AGPs FermiNet and the original FermiNet, we expand an AGPs determinant constructed using the pairing function from Eq. (<ref>) as a sum of block-diagonal determinants of FermiNet orbitals.
It will be sufficient to consider matrices $W \in \mathbb{R}^{n_L \times n_L}$ of rank $M$, with $p\le M\le n_L$.
We can decompose any such matrix using rank factorization,

\begin{equation}
W_{kl} = \sum^M_{\gamma=1} F_{\gamma k} G_{\gamma l},
\label{eq:rf}
\end{equation}

where $F$ and $G$ are matrices in $\mathbb{R}^{M\times n_L}$.
Equation (<ref>) then becomes

\begin{align}
\varphi(\mathbf{r}^\uparrow_i, \mathbf{r}^\downarrow_j)
&= \sum_{\gamma=1}^{M} \left(\sum_{k=1}^{n_L} F_{\gamma k} h^{L\uparrow (k)}_i\right) \left(\sum_{l=1}^{n_L} G_{\gamma l}  h^{L\downarrow (l)}_j\right) \notag \\
&= \sum_{\gamma=1}^{M} \phi^\uparrow_\gamma(\mathbf{r}^\uparrow_i) \phi^\downarrow_\gamma(\mathbf{r}^\downarrow_j),
\end{align}

where the last line defines the functions $\phi_{\gamma}^{\uparrow}$ and $\phi_{\gamma}^{\downarrow}$.
In the case when $M = p$, where $p = N/2$ is the number of pairs in the system, the determinant of the pairing function can be written as a block-diagonal determinant of FermiNet orbitals:

\begin{align}
\det[\varphi(\mathbf{r}^\uparrow_i, \mathbf{r}^\downarrow_j)]
&= \det[\sum_{\gamma=1}^p \phi^\uparrow_\gamma(\mathbf{r}^\uparrow_i) \phi^\downarrow_\gamma(\mathbf{r}^\downarrow_j)] \notag \\
&= \det[\left(\phi^{\uparrow T} \phi^{\downarrow}\right)_{ij}] \notag \\
&= \det[\phi^\uparrow_i(\mathbf{r}^\uparrow_j)] \det[\phi^\downarrow_i(\mathbf{r}^\downarrow_j)]
\label{eq:M_equals_p}
\end{align}

where $[\phi^{\alpha}]_{ij} = \phi^\alpha_i(\mathbf{r}^\alpha_j)$ are matrices in $\mathbb{R}^{M\times p}$ with $M=p$.
The product of two $p\times p$ determinants can be written as the determinant of a single $2p \times 2p$ matrix, with the $p\times p$ spin-up and spin-down blocks on the diagonal.
Therefore, any AGPs wavefunction constructed using the pairing function from Eq. (<ref>) with a rank-$p$ matrix $W_{kl}$ is equivalent to a $2p \times 2p$ block-diagonal determinant of conventional FermiNet orbitals.
The equivalence is already well known <cit.> for APGs wavefunctions constructed using conventional two-particle orbitals.

Now consider the more general case where $p\le M\le n_L$.
The Cauchy-Binet formula states that

\begin{align}
&\det[\varphi(\mathbf{r}^\uparrow_i, \mathbf{r}^\downarrow_j)]
= \det[\sum_{\gamma=1}^{M} \phi^\uparrow_\gamma(\mathbf{r}^\uparrow_i) \phi^\downarrow_\gamma(\mathbf{r}^\downarrow_j)] \notag \\
&= \sum_{1 \le j_1 < j_2 < \cdots < j_p \le M}
\det
\begin{pmatrix}
\phi^\uparrow_{j_1}(\mathbf{r}^\uparrow_1)
& \phi^\uparrow_{j_1}(\mathbf{r}^\uparrow_2)
& \cdots & \phi^\uparrow_{j_1}(\mathbf{r}^\uparrow_p)
\\
\phi^\uparrow_{j_2}(\mathbf{r}^\uparrow_1)
& \phi^\uparrow_{j_2}(\mathbf{r}^\uparrow_2)
& \cdots & \phi^\uparrow_{j_2}(\mathbf{r}^\uparrow_p)
\\
\vdots & \vdots & \ddots & \vdots
\\
\phi^\uparrow_{j_p}(\mathbf{r}^\uparrow_1)
& \phi^\uparrow_{j_p}(\mathbf{r}^\uparrow_2)
& \cdots & \phi^\uparrow_{j_p}(\mathbf{r}^\uparrow_p)
\end{pmatrix}
\det
\begin{pmatrix}
\phi^\downarrow_{j_1}(\mathbf{r}^\downarrow_1)
& \phi^\downarrow_{j_1}(\mathbf{r}^\downarrow_2)
& \cdots & \phi^\downarrow_{j_1}(\mathbf{r}^\downarrow_p)
\\
\phi^\downarrow_{j_2}(\mathbf{r}^\downarrow_1)
& \phi^\downarrow_{j_2}(\mathbf{r}^\downarrow_2)
& \cdots & \phi^\downarrow_{j_2}(\mathbf{r}^\downarrow_p)
\\
\vdots & \vdots & \ddots & \vdots
\\
\phi^\downarrow_{j_p}(\mathbf{r}^\downarrow_1)
& \phi^\downarrow_{j_p}(\mathbf{r}^\downarrow_2)
& \cdots & \phi^\downarrow_{j_p}(\mathbf{r}^\downarrow_p)
\end{pmatrix},
\end{align}

where the sum is over all $\binom{M}{p}$ distinct choices of $p$ rows from the two $M \times p$ matrices $\phi^\uparrow_\gamma(\mathbf{r}^\uparrow_i)$ and $\phi^\downarrow_\gamma(\mathbf{r}^\downarrow_i)$.
The products of the determinants of the two $p \times p$ matrices associated with each such choice are summed to reproduce the APGs.
This is similar to the block-diagonal multiple determinant expansion of FermiNet orbitals without weights given by Eq. (<ref>) and in the original FermiNet paper <cit.>
To get a rough estimate of the number of terms in the sum, we can take the number of hidden units in the one-electron stream to be $M = n_L \sim 200$ and the number of pairs to be $p \sim 20$, this number is approximately $10^{27}$.

Note that the intermediate layers, i.e., the one and two-electron streams, are identical in the original FermiNet and the AGPs FermiNet.
The only modifications are made at the orbital shaping layer, or, equivalently, the method of antisymmetrization has changed.
Thus, the representational power of the intermediate layers of the AGPs FermiNet remain the same as for the original FermiNet. Thus, it must be the method of antisymmetrization that limits the performance of the original FermiNet when applied to the UFG.

We have shown that a single APGs determinant constructed using the pairing function from Eq. (<ref>) with a matrix $W_{kl}$ of rank greater than $p$ contains multiple block-diagonal determinants of FermiNet orbitals.
If the rank of $W_{kl}$ is equal to $p$, the AGPs is equivalent to a single block-diagonal FermiNet determinant.
Conversely, any single-determinant FermiNet wavefunction can be written as an AGPs of rank $p$.
Therefore, the AGPs FermiNet provide a more powerful Ansatz with fewer variational parameters than the original FermiNet, since the former contains the latter.

It is worth mentioning another advantage of using the FermiNet to build geminals.
By generating more sets of independent parameters $\mathbf{w}^{k}_i$ in Eq. (<ref>), one can easily construct an arbitrary number $N_{\det}$ of FermiNet geminals $\varphi^{k}(\mathbf{r}^\alpha_i, \mathbf{r}^{\bar\alpha}_j;
\{\mathbf{r}^\alpha_{/i}\};\{\mathbf{r}^{\bar\alpha}_{/j}\})$ with $k \in \{1, 2, \dots, N_{\det}\}$, all without the use of a basis set.
It allows one to use weighted sums of AGPs determinants as trial wavefunctions, similar to the weighted sum of conventional FermiNet determinants seen in Eq. (<ref>).

§.§ AGPs FermiNet with Unpaired States

To extend the AGPs FermiNet to systems with unpaired states, such as an odd-number of particle system, we use FermiNet geminals and orbitals to replace both the pairing orbitals and the single-particle orbitals in Eq. (<ref>).
In this work we consider systems with equal numbers of up-spin and down-spin particles, which are assumed to be fully paired, and systems containing one additional unpaired particle, which may have spin up or spin down.
For example, the AGPs FermiNet with an extra spin-up particle is given by

\begin{align}
&\Psi^\text{AGPs}_\text{FermiNet}\left(\mathbf{r}^\uparrow_1, \dots, \mathbf{r}^\uparrow_{p+1}, \mathbf{r}^\downarrow_{1}, \dots,  \mathbf{r}^\downarrow_{p}\right) \notag \\
&= \sum_k \omega_k \det \left(
\begin{array}{cccc}
\varphi^k(\mathbf{r}^\uparrow_1, \mathbf{r}^\downarrow_1;
\{\mathbf{r}^\uparrow_{/1}\};\{\mathbf{r}^{\downarrow}_{/1}\})
& \cdots
& \varphi^k(\mathbf{r}^\uparrow_1, \mathbf{r}^\downarrow_{p};
\{\mathbf{r}^\uparrow_{/1}\};\{\mathbf{r}^{\downarrow}_{/p}\})
& \phi^{k\uparrow}_1(\mathbf{r}^\uparrow_1;\{\mathbf{r}^\uparrow_{/1}\}; \{\mathbf{r}^{\downarrow}\})
\\
\varphi^k(\mathbf{r}^\uparrow_2, \mathbf{r}^\downarrow_1;
\{\mathbf{r}^\uparrow_{/2}\};\{\mathbf{r}^{\downarrow}_{/1}\})
& \cdots
& \varphi^k(\mathbf{r}^\uparrow_2, \mathbf{r}^\downarrow_{p};
\{\mathbf{r}^\uparrow_{/2}\};\{\mathbf{r}^{\downarrow}_{/p}\})
& \phi^{k\uparrow}_1(\mathbf{r}^\uparrow_2;\{\mathbf{r}^\uparrow_{/2}\}; \{\mathbf{r}^{\downarrow}\})
\\
\vdots
& \ddots
& \vdots
& \vdots
\\
\varphi^k(\mathbf{r}^\uparrow_{p+1}, \mathbf{r}^\downarrow_1;
\{\mathbf{r}^\uparrow_{/p+1}\};\{\mathbf{r}^{\downarrow}_{/1}\})
& \cdots
& \varphi^k(\mathbf{r}^\uparrow_{p+1}, \mathbf{r}^\downarrow_{p};
\{\mathbf{r}^\uparrow_{/p+1}\};\{\mathbf{r}^{\downarrow}_{/p}\})
& \phi^{k\uparrow}_1(\mathbf{r}^\uparrow_{p+1};\{\mathbf{r}^\uparrow_{/p+1}\}; \{\mathbf{r}^{\downarrow}\})
\end{array}
\right)
\label{eq:agpsferminet_uppaired}
\end{align}

Having explained how to construct AGPs wavefunctions for fully paired systems and for systems with a single unpaired particle, we are ready to used them to study the UFG.