# Oscillatory force autocorrelations in equilibrium odd-diffusive systems

Erik Kalz  University of Potsdam, Institute of Physics and Astronomy, D-14476
Potsdam, Germany  Hidde Derk Vuijk  University of Augsburg, Institute of
Physics, D-86159 Augsburg, Germany  Jens-Uwe Sommer  Leibniz-Institute for
Polymer Research, Institute Theory of Polymers, D-01069 Dresden, Germany
Technical University of Dresden, Institute for Theoretical Physics, D-01069
Dresden, Germany  Technical University of Dresden, Cluster of Excellence
Physics of Life, D-01069 Dresden, Germany  Ralf Metzler  University of
Potsdam, Institute of Physics and Astronomy, D-14476 Potsdam, Germany  Asia
Pacific Centre for Theoretical Physics, KR-37673 Pohang, Republic of Korea
Abhinav Sharma  University of Augsburg, Institute of Physics, D-86159
Augsburg, Germany  Leibniz-Institute for Polymer Research, Institute Theory of
Polymers, D-01069 Dresden, Germany

###### Abstract

The force autocorrelation function (FACF), a concept of fundamental interest
in statistical mechanics, encodes the effect of interactions on the dynamics
of a tagged particle. In equilibrium, the FACF is believed to decay
monotonically in time which is a signature of slowing down of the dynamics of
the tagged particle due to interactions. Here we analytically show that in
odd-diffusive systems, which are characterized by a diffusion tensor with
antisymmetric elements, the FACF can become negative and even exhibit temporal
oscillations. We also demonstrate that, despite the isotropy, the knowledge of
FACF alone is not sufficient to describe the dynamics: the full
autocorrelation tensor is required and contains an antisymmetric part. These
unusual properties translate into enhanced dynamics of the tagged particle
quantified via the self-diffusion coefficient that, remarkably, increases due
to particle interactions.

Introduction. Time integrals of appropriate correlation functions are related
to transport coefficients via the Green-Kubo relations [1, 2]. The self-
diffusion coefficient, for instance, is determined by the integral of the
force autocorrelation function (FACF), which quantifies the effect of
interactions on the diffusive dynamics of a particle. Previous studies have
shown that in overdamped equilibrium systems, the FACF decays monotonically
for all densities independently of the nature of interaction between particles
[3, 4, 5, 6]. As a consequence, the self-diffusion coefficient is always
reduced in interacting systems. A nonmonotonic decay is known to exist in
active [7], driven [8], and harmonically trapped systems [9], as well as in
fluid systems with inertia [10], but it has never been shown in equilibrium
systems. In fact, such a behavior was even shown to be prohibited in
overdamped equilibrium systems [11, 12, 7]. The lack of any nonmonotonic
features is intuitively expected, for there exists neither inertia nor any
internal or external driving which could introduce additional time scales to
the system.

Here we show for the first time that the FACF can be nonmonotonic and even
oscillatory in overdamped equilibrium systems. Systems showing this behavior
are characterized by probability fluxes, which are perpendicular to density
gradients and are referred to as _odd-diffusive_ systems [13]. We further
demonstrate that the unusual behavior of the self-diffusion coefficient in
these systems, it increases with increasing density [14], is a natural
consequence of the nonmonotonicity of the FACF. The transverse response to the
perturbation is the fundamental property of _odd_ systems which have received
much interest lately [15]. In addition to odd-diffusive systems, there are odd
systems characterized by odd viscosity [16, 17, 18, 19, 20, 21], odd
elasticity [22, 23] and odd viscoelasticity [24, 25]. With the advent of
experimental odd systems such as spinning biological organisms [26], chiral
fluids [27, 28] and colloidal spinners [29], the interest in odd systems has
increased rapidly.

The odd-diffusion tensor for a two-dimensional isotropic system can be written
as

$\mathbf{D}=D_{0}\left(\mathbf{1}+\kappa\boldsymbol{\varepsilon}\right),$ (1)

where $\mathsf{1}$ is the identity matrix, $\boldsymbol{\varepsilon}$ is the
antisymmetric Levi-Civita symbol in two dimensions
($\boldsymbol{\varepsilon}_{xy}=-\boldsymbol{\varepsilon}_{yx}=1$ and
$\boldsymbol{\varepsilon}_{xx}=\boldsymbol{\varepsilon}_{yy}=0$), $D_{0}$ is
the diffusivity and $\kappa$ is the odd-diffusion parameter. A nonzero
$\kappa$ results in probability fluxes perpendicular to density gradients.
Examples of odd-diffusive systems are Brownian particles diffusing under the
effect of Lorentz force [30, 31, 32, 33, 34, 35], and diffusing skyrmions [36,
37, 38, 39, 40, 41], see also the Supplementary Material (SM) [42]. Although
these are equilibrium odd-diffusive systems, there exist also driven odd-
diffusive systems such as active chiral particles (also called circle
swimmers) [43, 44, 45, 46] and strongly damped particles subjected to Magnus
[47] or Coriolis force [48]. In contrast to equilibrium systems which are
invariant under time-reversal, the odd-diffusive behavior in nonequilibrium
systems is a consequence of broken time-reversal and parity symmetries [21].

While an exact calculation of the FACF is a formidable task, near-exact
analytical results can be obtained in the dilute limit in which the dynamics
are dominated by two-body effects. To this end, we generalize the first-
principles approach developed by Hanna, Hess, and Klein [3, 4] to calculate
the FACF in a dilute odd-diffusive system of hard-core interacting particles.
We show analytically that odd diffusion qualitatively alters the time
correlations: the correlation function becomes negative for finite $\kappa$
indicating the anticorrelated nature of the force experienced by an odd-
diffusive particle due to collisions with other particles. Moreover, the
correlation function exhibits temporal oscillations for certain values of
$\kappa$ it crosses zero twice. We further show that for sufficiently large
$\kappa$, the integral of the correlation function becomes negative which
gives rise to the increase in the self-diffusion coefficient. Using the Green-
Kubo relation, we derive exactly the same expression for the self-diffusion
coefficient as in Ref. [14] which was obtained using an alternative approach.

Theoretical background. We consider a two-dimensional system of two
interacting, odd-diffusive hard disks with coordinates
$\vec{\mathbf{x}}=(\mathbf{x}_{1},\mathbf{x}_{2})$. The two-particle
conditional probability density function for the particles to evolve from
$\vec{\mathbf{x}}^{\prime}$ at time $t^{\prime}\leq t$ to $\vec{\mathbf{x}}$
at time $t$, $P=P(\vec{\mathbf{x}},t|\vec{\mathbf{x}}^{\prime},t^{\prime})$,
satisfies the Fokker-Planck equation

$\displaystyle\frac{\partial}{\partial t}P$
$\displaystyle=\nabla_{1}\cdot\mathbf{D}\left[\nabla_{1}+\beta\nabla_{1}U(r)\right]P$
$\displaystyle\quad+\nabla_{2}\cdot\mathbf{D}\left[\nabla_{2}+\beta\nabla_{2}U(r)\right]P,$
(2)

with the odd-diffusion tensor (1) and $\nabla_{1},\nabla_{2}$ as the partial
differential operator with respect to the coordinates of particle one and two,
respectively. $U(r)$ is the potential energy with
$r=|\mathbf{x}_{1}-\mathbf{x}_{2}|$ as the relative distance between the
particles and $\beta=1/k_{\mathrm{B}}T$, where $k_{\mathrm{B}}$ is the
Boltzmann constant and $T$ is the temperature. We assume hard-core
interactions between the two disks of diameter $\sigma$, which can be written
as $U(r)=\small\begin{cases}\textstyle\infty,&r\leq\sigma\\\
0,&r>\sigma\end{cases}$. The analytical solution to the two-particle Fokker-
Planck equation was obtained for normal-diffusing particles, i.e.
$\mathbf{D}=D_{0}\mathbf{1}$ [3, 4]. While the hard-core interactions are
modeled via Neumann boundary conditions in normal-diffusing systems, they are
modeled as oblique boundary conditions in odd-diffusive systems due to the
transverse fluxes [49, 14]. This has profound consequences for the solution
and therefore for the application of our theory. We solve the two-particle
problem (Oscillatory force autocorrelations in equilibrium odd-diffusive
systems) for odd-diffusive hard disks exactly in the SM [42].

Force autocorrelation tensor. The force autocorrelation tensor (FACT), which
is defined as
$\mathsf{C}_{F}(\tau)=\left<\mathbf{F}(\tau)\otimes\mathbf{F}(0)\right>$, can
be written as [50]

$\displaystyle\mathsf{C}_{F}(\tau)$
$\displaystyle=\int\mathrm{d}\vec{\mathbf{x}}\int\mathrm{d}\vec{\mathbf{x}}_{0}~{}\mathbf{F}\left(\vec{\mathbf{x}}\right)\otimes\mathbf{F}\left(\vec{\mathbf{x}}_{0}\right)$
$\displaystyle\quad\times
P\left(\vec{\mathbf{x}},\tau|\vec{\mathbf{x}}_{0},0\right)P_{\mathrm{eq}}\left(\vec{\mathbf{x}}_{0}\right),$
(3)

for $\tau>0$. Here $\mathbf{F}$ is the interaction force acting on a tagged
particle due to other particles, $\langle\cdot\rangle$ denotes an ensemble
average with the equilibrium distribution
$P_{\mathrm{eq}}\left(\vec{\mathbf{x}}_{0}\right)$ and the outer product is
defined as $[\mathbf{A}\otimes\mathbf{B}]_{\alpha\beta}=A_{\alpha}B_{\beta}$.
Throughout the paper, time is measured in units of
$\tau_{0}=\sigma^{2}/(2D_{0})$ which is the characteristic timescale of a
particle diffusing over a distance of diameter $\sigma$, i.e.
$\tau=t/\tau_{0}.$ The FACT can be calculated from Eq. (3) to first order in
the density, details of which are shown in SM [42]. Similar to the diffusion
tensor, the FACT can be split in a diagonal and an antisymmetric off-diagonal
part:

$\mathsf{C}_{F}(\tau)=C_{F}^{\mathrm{diag}}(\tau)\mathsf{1}+C_{F}^{\mathrm{off}}(\tau)\boldsymbol{\varepsilon},$
(4)

for $\tau>0$, where $C_{F}^{\mathrm{diag}}(\tau)$ and
$C_{F}^{\mathrm{off}}(\tau)$ are the diagonal and antisymmetric off-diagonal
elements of the FACT. In Laplace domain they read

$\displaystyle\tilde{C}_{F}^{\mathrm{diag}}(s)$ $\displaystyle=\frac{2\
\phi}{\beta^{2}D_{0}}\
\frac{K_{1}\left[\sqrt{s}K_{0}+K_{1}\right]}{\left[\sqrt{s}K_{0}+K_{1}\right]^{2}+\left[\kappa
K_{1}\right]^{2}},$ (5) $\displaystyle\tilde{C}_{F}^{\mathrm{off}}(s)$
$\displaystyle=\frac{2\ \phi}{\beta^{2}D_{0}}\ \frac{\kappa\
[K_{1}]^{2}}{\left[\sqrt{s}K_{0}+K_{1}\right]^{2}+\left[\kappa
K_{1}\right]^{2}},$ (6)

where $K_{n}=K_{n}(\sqrt{s})$ is the modified Bessel function of the second
kind of order $n$, $\phi=\pi\ (N/V)\ (\sigma/2)^{2}$ is the area fraction for
$N$ particles of diameter $\sigma$ in an area $V$, and $\tilde{(\cdot)}$
denotes the Laplace transform with $s$ as the Laplace variable conjugate to
$\tau$. Note that the off-diagonal elements $C^{\mathrm{off}}_{F}$ are
proportional to the odd-diffusion parameter $\kappa$ and therefore vanish in
the case of normal diffusion ($\kappa=0$). In this case the FACT reduces to
$\mathsf{C}_{F}(\tau)=C_{F}^{\mathrm{diag}}(\tau)\
\mathsf{1}=\frac{\mathsf{1}}{2}\langle\mathbf{F}(\tau)\cdot\mathbf{F}(0)\rangle$,
which is the usual FACF in normal systems.

Figure 1: Double-logarithmic plot of the diagonal and off-diagonal elements of
the force autocorrelation tensor (FACT) of interacting hard disks as a
function of reduced time $\tau=t/\tau_{0}$, where
$\tau_{0}=\sigma^{2}/(2D_{0})$. (a) The diagonal elements of the FACT
$C_{F}^{\mathrm{diag}}(\tau)$, corresponding to the force autocorrelation
function (FACF), can turn negative. The FACF diverges in the limit $\tau\to 0$
as $C_{F}^{\mathrm{diag}}(\tau)\simeq\tau^{-1/2}$. At long times the FACF
scales as $C_{F}^{\mathrm{diag}}(\tau)\simeq\tau^{-2}$. For $\kappa=1$ we find
an exceptional long-time behavior, where
$C_{F}^{\mathrm{diag}}(\tau)\simeq\tau^{-3}$. The inset shows the zero-
crossing time $\tau_{c}$ of $C_{F}^{\mathrm{diag}}(\tau)$ as a function of
$\kappa$, which in the main figure is marked by red circles. The onset of the
anticorrelation corresponds to $\kappa>\kappa_{th}\approx 0.88$ (b) The off-
diagonal elements of the FACT $C_{F}^{\mathrm{off}}(\tau)$ are independent of
time in the short-time limit $C_{F}^{\mathrm{off}}(\tau)\simeq\kappa\tau^{0}$
and are directly proportional to $\kappa$. In the long-time limit, they scale
similarly to the diagonal elements as
$C_{F}^{\mathrm{off}}(\tau)\simeq\tau^{-2}$ for all $\kappa$. The inset in (b)
shows typical configurations after a collision of particles, where the
orientational change of the force (orange arrow)
$\mathbf{F}_{i}=\mathbf{F}(\tau_{i}),\ i\in\\{0,1,2,3\\}$ of the tagged
particle (red) is indicated.

The diagonal and off-diagonal elements of the FACT are plotted in Fig.1 as a
function of time. We first consider the behavior of the diagonal elements of
the tensor in Fig.1(a), which correspond to the usual FACF for odd-diffusive
systems. For small values of $\kappa$, the FACF is a positive, monotonically
decaying function of time, qualitatively similar to a normal diffusive system.
For larger values of $\kappa$, however, a new feature appears in the FACF: it
crosses through zero and hence becomes negative, indicating an anticorrelation
of the force. The time scale of the force reversal on a tracer particle, i.e.,
when the FACF becomes negative, depends strongly on $\kappa$, as can be seen
in the inset of Fig.1(a). There exists a numerically obtained threshold
$\kappa_{th}\approx 0.88$ below which the FACF is strictly positive. The off-
diagonal elements of the FACT are shown in Fig.1(b). Unlike the diagonal
elements, which diverge as $t\to 0$, the off-diagonal elements remain finite.
Specifically they remain positive for all $\kappa$ and decay monotonically in
time.

It is interesting to investigate the short- and long-time behavior of the
elements of the FACT. Using the asymptotic behavior of the modified Bessel
functions $K_{0}$ and $K_{1}$, see SM [42] for details, from Eq. (5) and Eq.
(6) we have analytical access to the behavior on time scales $t\ll\tau_{0}$
and $t\gg\tau_{0}$, i.e. $s\gg 1$ and $s\ll 1$ in the Laplace domain,
respectively. At short times, the FACF behaves like
$C^{\mathrm{diag}}_{F}(\tau)\simeq\tau^{-1/2}$, as shown in Fig.1(a), and is
independent of $\kappa$. Here $\simeq$ is used to denote asymptotic
proportionality. The long-time behavior of the FACF can be obtained from the
$s\ll 1$ expansion and behaves asymptotically as

$\tilde{C}^{\mathrm{diag}}_{F}(s)\sim\frac{2\phi}{\beta^{2}D_{0}}\frac{1}{1+\kappa^{2}}\Bigg{(}1+\frac{1-\kappa^{2}}{1+\kappa^{2}}\bigg{(}\gamma-\ln(2)+\frac{\ln(s)}{2}\bigg{)}s\\\
+\frac{1-6\kappa^{2}+\kappa^{4}}{8(\kappa^{2}+1)^{3}}\
s^{2}\ln^{2}(s)\Bigg{)},$ (7)

for $s\to 0$ and where $\gamma=0.5772$ is the Euler-Mascheroni constant. For
$\kappa=0$, the asymptotic behavior of $\tilde{C}^{\mathrm{diag}}_{F}$
coincides with the form reported for related 2d Lorentz gas systems [51].
Furthermore, from Eq. (7) it can be seen that the long-time behavior of
$C_{F}^{\mathrm{diag}}(\tau)$ strongly depends on $\kappa$. The FACF decays as
$\tau^{-2}$ for all $\kappa$ except for $\kappa=1$, at which the leading order
contribution vanishes in Eq. (7) and
$C_{F}^{\mathrm{diag}}(\tau)\simeq\tau^{-3}\ln(\tau)$, as shown in Fig.1(a)
[52, 53]. The ordinary algebraic long-time decay $\simeq\tau^{-2}$
($\kappa\neq 1$) is consistent with the general prediction of a decay
$\simeq\tau^{-(d/2+1)},\ d=1,2,3,$ for correlation functions in systems, which
do not conserve momentum [54, 55]. This universal behavior was theoretically
and numerically exhaustively demonstrated specifically for the 2d Lorentz gas
model [56, 57, 58, 59, 51]. In three dimensions, the decay of the correlation
functions $\simeq\tau^{-5/2}$ [60, 3, 61, 62, 63] could recently be
demonstrated computationally [6]. In contrast, the short-time behavior
$\simeq\tau^{-1/2}$ is independent of dimensionality and attributed to the
hard interactions between the particles [3, 61].

The asymptotic short-time behavior of $C_{F}^{\mathrm{off}}(\tau)$ turns out
to be independent of time but depends linearly on $\kappa$,
$C_{F}^{\mathrm{off}}(\tau)\simeq\kappa\tau^{0}$, as as can be seen in
Fig.1(b). Such a scaling of the off-diagonal elements with $\kappa$ at short
times has been recently derived by Yasuda et. al in Ref. [64] for odd Langevin
systems. The authors also pointed out that this could be useful for estimating
the odd-diffusion parameter in experiments. The asymptotic long-time behavior
of $C_{F}^{\mathrm{off}}(\tau)$ shows a monotonic decay in time and also
depends on $\kappa$,
$C_{F}^{\mathrm{off}}(\tau)\simeq\kappa\tau^{-2}/(\kappa^{2}+1)^{2}$, as can
be seen in Fig.1(b).

In a low-density system, in which only two-body correlations are important, it
is quite surprising that the FACF can turn negative, as shown in Fig.1. It is
even more surprising that there exists a range of $\kappa\in(\kappa_{th},1)$
for which the FACF exhibits not one but two zero crossings, as shown in Fig.
2. It appears that for $\kappa$ slightly larger than $\kappa_{th}\approx
0.88$, which is obtained from numerical inversion of Eq. (5), the FACF first
becomes anticorrelated (first zero crossing) in time before it crosses the
time axis again (second and last zero crossing). Here, at long times, the FACF
decays to zero from above. We have numerically inverted the Laplace transform
over much longer times than shown here and did not find more than two zero
crossings. This ”temporal oscillation” in the FACF ceases to exist for
$\kappa\geq 1$. For $\kappa>1$, the asymptotic expansion in Eq. (7),
transformed back into time domain, is strictly negative and therefore the FACF
decays to zero from below, i.e. the second zero-crossing vanishes (see also
inset in Fig.2).

Green-Kubo relation for the self-diffusion coefficient. The self-diffusion
coefficient $D_{s}$ can be obtained from the velocity autocorrelation function
(VACF) $C_{v}(\tau)=\langle\mathbf{v}(\tau)\cdot\mathbf{v}(0)\rangle/2$, where
$\mathbf{v}(\tau)=\mathrm{d}\mathbf{x}/\mathrm{d}\tau$ and $\mathbf{x}$ is the
position of the tagged particle as the time integral

$D_{s}=\int_{0}^{\infty}\mathrm{d}t\ C_{v}(\tau),$ (8)

a Green-Kubo relation between an equilibrium autocorrelation function
($C_{v}(\tau)$) and a transport coefficient ($D_{s}$) [1].

Figure 2: Double-logarithmic plot of the absolute value of the diagonal
elements of the force autocorrelation tensor $C_{F}^{\mathrm{diag}}(\tau)$ of
interacting hard disks as a function of reduced time $\tau=t/\tau_{0}$, where
$\tau_{0}=\sigma^{2}/(2D_{0})$. Investigating the regime
$\kappa\in[0.88,1.0]$, we find oscillatory behavior of
$C_{F}^{\mathrm{diag}}(\tau)$. At short times $C_{F}^{\mathrm{diag}}(\tau)$
starts as a positive function, turns negative and after a second zero-crossing
becomes positive again. The inset shows the zero-crossing times $\tau_{c}$ of
$C_{F}^{\mathrm{diag}}(\tau)$ as a function of $\kappa$ in a linear-
logarithmic plot which in the main figure are marked as red circles. The
oscillatory behavior starts at $\kappa\geq\kappa_{th}=0.88$, whereas the
second zero-crossing drifts to infinity as $\kappa\to 1$. For $\kappa>1$,
$C_{F}^{\mathrm{diag}}(\tau)$ only shows one zero-crossing and remains anti-
correlated for the remaining $\tau\to\infty$ (see also inset in Fig.1 (a)).

In normal diffusive systems, the VACF is related to the FACF. In contrast, in
an odd-diffusive system, the knowledge of the FACF alone is not sufficient to
calculate the VACF. This is despite the fact that the system is isotropic. In
fact, one requires the entire FACT to calculate the velocity correlation
function. We show in SM [42] that in odd-diffusive systems, the VACF can be
written as

$C_{v}(\tau)=D_{0}\left(\delta_{+}(\tau)-D_{0}\beta^{2}C_{F}(\tau)\right),$
(9)

where

$C_{F}(\tau)=\frac{1}{2}\frac{1}{D_{0}^{2}}(\mathsf{D}^{2})^{\mathrm{T}}:\mathsf{C}_{F}(\tau),$
(10)

and where the double contraction is defined as
$\mathsf{A}:\mathsf{B}=\sum_{\alpha,\beta=1}^{2}A_{\alpha\beta}B_{\beta\alpha}$.
$\delta_{+}$ is the one-sided delta distribution, see also the SM [42]. We
refer to $C_{F}(\tau)$ as the generalized force autocorrelation function
(gFACF) which reads

$C_{F}(\tau)=\left(1-\kappa^{2}\right)C_{F}^{\mathrm{diag}}(\tau)-2\kappa
C_{F}^{\mathrm{off}}(\tau).$ (11)

For normal diffusive systems (i.e. $\kappa=0$), $C_{F}$ reduces to the
ordinary FACF. Note that even though the gFACF is diverging for all
$\kappa\neq 1$ in $\tau\rightarrow 0$ in the hard-disk system, the function
remains integrable. This is of physical significance since the integral of the
gFACF captures the effect of collisions on the self-diffusion as we see from
the Green-Kubo relation Eq. (8) together with Eq. (9).

The self-diffusion coefficient $D_{s}$ can be obtained from the time integral
of Eq. (9) or by using the limit theorem $\int_{0}^{\infty}f(t)\
\mathrm{d}t=\lim_{s\to 0}\tilde{f}(s)$ in Eq. (7) for
$\tilde{C}^{\mathrm{diag}}_{F}$ and similarly for
$\tilde{C}^{\mathrm{off}}_{F}$ which yields

$\displaystyle\lim_{s\to 0}\tilde{C}_{F}^{\mathrm{diag}}(s)$
$\displaystyle=\frac{1}{\kappa}\lim_{s\to
0}\tilde{C}_{F}^{\mathrm{off}}(s)=\frac{2\ \phi}{\beta^{2}D_{0}}\
\frac{1}{1+\kappa^{2}}.$ (12)

Together with Eq. (9) and Eq. (11), this gives the self-diffusion coefficient
in an odd-diffusive system,

$D_{s}=D_{0}\left(1-2\phi\ \frac{1-3\kappa^{2}}{1+\kappa^{2}}\right),$ (13)

valid up to first order in area concentration $\phi$ for a system of hard
disks. This result was previously derived by us in Refs. [14, 65] by a
different method.

For $\kappa=0$ the expression for $D_{s}$ reproduces the known result of
normal diffusive systems of hard disks in two dimensions
$D_{s}=D_{0}(1-2\phi)$ [4, 61]. The surprising result of $D_{s}$ in Eq. (13)
is that the prefactor of $\phi$ can change sign. This shows that odd
diffusivity ($\kappa>0$) results in a cancellation of the ordinary collision-
induced reduction of the self-diffusion. For $\kappa=\kappa_{c}=1/\sqrt{3}$,
up to first order in the area fraction, the effect of the collisions on the
self-diffusion vanishes ($D_{s}=D_{0}$), meaning that the on long time and
length scales hard disks appear to diffuse as non-interacting particles. For
$\kappa>\kappa_{c}$, collisions surprisingly increase the self-diffusion
coefficient: the system mixes more efficiently.

It is natural to ask whether our findings can be extended to three dimensions.
However, in three dimensions, odd systems cannot be isotropic because the
plane in which the rotation takes place breaks isotropy [15, 66, 13]. We
investigated the self-diffusion in such a system via Brownian dynamics
simulations and found that the in-plane odd diffusivity has no effect on the
diffusion along the axes of rotation, which turns out to be exactly the same
as that of a normal-diffusive system of hard spheres. The in-plane
diffusivity, however, shows the same $\kappa$-dependent behavior as in a two-
dimensional odd-diffusive system.

Discussion. We analytically demonstrated that equilibrium correlation
functions can be non-monotonic and even oscillatory in overdamped systems.
This finding is at odds with the statement that in an equilibrium system the
correlation function and all its derivatives decay monotonically [12, 7].
While the latter holds in systems where the time-evolution is described by a
Hermitian Fokker-Planck operator, for odd systems this is not applicable due
to their intrinsic antisymmetric off-diagonal elements in the diffusion tensor
(1).

Our work shows that rich physics is to be explored in equilibrium, odd-
diffusive systems. In normal-diffusive systems, for instance, there exists a
crossover between two diffusive regimes: short-time diffusion with diffusivity
$D_{0}$ and long-time diffusion with $D_{s}<D_{0}$ [50]. That the long-time
self-diffusion coefficient is smaller than the short-time is indicative of the
slowing down of the dynamics of the tracer particle in the crossover. In odd-
diffusive systems, in contrast, the dynamics can be enhanced, which is
reflected in the anticorrelated force autocorrelations. The anticorrelation
can be physically interpreted in terms of reversal of the force experienced by
a tagged particle such that rather than impeding, collisions with other odd-
diffusive particles enhance the motion of the tagged particle, see also the
inset in Fig. 1(b). Even though qualitatively this mutual rolling of particles
explains the enhancement of self-diffusion with collisions in an odd-diffusive
system through the reversal of force [14], a detailed mechanism is still
elusive. To this end, we believe it will be interesting to investigate the
structural rearrangements that occur in an odd-diffusive system and contrast
them with those in a normal diffusive system. We further expect that the
unusual behavior could also have implications for the rheological properties
of odd fluids, such as viscosity.

With increasing experimental interest in systems such as spinning biological
organisms [26], chiral fluids [27, 28], and colloidal spinners [29], our work
will contribute to the broadening interest of the physics community in these
systems, especially in the novel and interesting way how interactions modify
the particle dynamics here. Furthermore, we believe that our work will
stimulate fundamental research on extending statistical physics to the novel
case of odd-diffusive systems. Lastly, since exact analytical results are
rather rare in interacting systems, our work may serve as a reference to
validate approximate theories for dense systems or computer simulations.

###### Acknowledgements.

Acknowledgements. We would like to thank one of the anonymous Reviewers for
several suggestions and fruitful comments regarding the analytical results of
the auto-correlation functions discussed in this work. We further thank Felix
Büttner for illuminating discussions on skyrmionic systems. E. K., R. M., and
A. S. acknowledge support by the Deutsche Forschungsgemeinschaft (grants No.
ME 1535/16-1 and SH 1275/3-1). J.-U.S. thanks the cluster of excellence
“Physics of Life” at TU Dresden for support.

## References

  * Green [1954] M. S. Green, Markoff random processes and the statistical mechanics of time-dependent phenomena. ii. irreversible processes in fluids, J. Chem. Phys. 22, 398 (1954).
  * Kubo  _et al._ [2012] R. Kubo, M. Toda, and N. Hashitsume,  _Statistical physics II: nonequilibrium statistical mechanics_ , Vol. 31 (Springer Science & Business Media, 2012).
  * Hanna  _et al._ [1981] S. Hanna, W. Hess, and R. Klein, The velocity autocorrelation function of an overdamped brownian system with hard-core interaction, J. Phys. A 14, L493 (1981).
  * Hanna  _et al._ [1982] S. Hanna, W. Hess, and R. Klein, Self-diffusion of spherical brownian particles with hard-core interaction, Physica A 111, 181 (1982).
  * Sharma and Brader [2016] A. Sharma and J. M. Brader, Communication: Green-kubo approach to the average swim speed in active brownian systems, J. Chem. Phys. 145, 161101 (2016).
  * Mandal  _et al._ [2019] S. Mandal, L. Schrack, H. Löwen, M. Sperl, and T. Franosch, Persistent anti-correlations in brownian dynamics simulations of dense colloidal suspensions revealed by noise suppression, Phys. Rev. Lett. 123, 168001 (2019).
  * Caraglio and Franosch [2022] M. Caraglio and T. Franosch, Analytic solution of an active brownian particle in a harmonic well, Phys. Rev. Lett. 129, 158001 (2022).
  * Leitmann  _et al._ [2018] S. Leitmann, S. Mandal, M. Fuchs, A. M. Puertas, and T. Franosch, Time-dependent active microrheology in dilute colloidal suspensions, Phys. Rev. Fluids 3, 103301 (2018).
  * Franosch  _et al._ [2011] T. Franosch, M. Grimm, M. Belushkin, F. M. Mor, G. Foffi, L. Forró, and S. Jeney, Resonances arising from hydrodynamic memory in brownian motion, Nature 478, 85 (2011).
  * Jeney  _et al._ [2008] S. Jeney, B. Lukić, J. A. Kraus, T. Franosch, and L. Forró, Anisotropic memory effects in confined colloidal diffusion, Phys. Rev. Lett. 100, 240604 (2008).
  * Feller [1970] W. Feller,  _An introduction to probability theory and its applications_ (John Wiley & Sons, New York, 1970).
  * Leitmann and Franosch [2017] S. Leitmann and T. Franosch, Time-dependent fluctuations and superdiffusivity in the driven lattice lorentz gas, Phys. Rev. Lett. 118, 018001 (2017).
  * Hargus  _et al._ [2021] C. Hargus, J. M. Epstein, and K. K. Mandadapu, Odd diffusivity of chiral random motion, Phys. Rev. Lett. 127, 178001 (2021).
  * Kalz  _et al._ [2022] E. Kalz, H. D. Vuijk, I. Abdoli, J.-U. Sommer, H. Löwen, and A. Sharma, Collisions enhance self-diffusion in odd-diffusive systems, Phys. Rev. Lett. 129, 090601 (2022).
  * Fruchart  _et al._ [2023] M. Fruchart, C. Scheibner, and V. Vitelli, Odd viscosity and odd elasticity, Annu. Rev. Condens. Matter Phys. 14, 471 (2023).
  * Banerjee  _et al._ [2017] D. Banerjee, A. Souslov, A. G. Abanov, and V. Vitelli, Odd viscosity in chiral active fluids, Nat. Commun. 8, 1 (2017).
  * Han  _et al._ [2021] M. Han, M. Fruchart, C. Scheibner, S. Vaikuntanathan, J. J. De Pablo, and V. Vitelli, Fluctuating hydrodynamics of chiral active fluids, Nat. Phys. 17, 1260 (2021).
  * Markovich and Lubensky [2021] T. Markovich and T. C. Lubensky, Odd viscosity in active matter: microscopic origin and 3d effects, Phys. Rev. Lett. 127, 048001 (2021).
  * Zhao  _et al._ [2022] Z. Zhao, M. Yang, S. Komura, and R. Seto, Odd viscosity in chiral passive suspensions, Front. Phys. , 815 (2022).
  * Lier  _et al._ [2022a] R. Lier, C. Duclut, S. Bo, J. Armas, F. Jülicher, and P. Surówka, Lift force in odd compressible fluids, arXiv preprint arXiv:2205.12704 (2022a).
  * Hargus  _et al._ [2020] C. Hargus, K. Klymko, J. M. Epstein, and K. K. Mandadapu, Time reversal symmetry breaking and odd viscosity in active fluids: Green–kubo and nemd results, J. Chem. Phys. 152, 201102 (2020).
  * Scheibner  _et al._ [2020] C. Scheibner, A. Souslov, D. Banerjee, P. Surówka, W. T. M. Irvine, and V. Vitelli, Odd elasticity, Nat. Phys. 16, 475 (2020).
  * Braverman  _et al._ [2021] L. Braverman, C. Scheibner, B. VanSaders, and V. Vitelli, Topological defects in solids with odd elasticity, Phys. Rev. Lett. 127, 268001 (2021).
  * Banerjee  _et al._ [2021] D. Banerjee, V. Vitelli, F. Jülicher, and P. Surówka, Active viscoelasticity of odd materials, Phys. Rev. Lett. 126, 138001 (2021).
  * Lier  _et al._ [2022b] R. Lier, J. Armas, S. Bo, C. Duclut, F. Jülicher, and P. Surówka, Passive odd viscoelasticity, Phys. Rev. E 105, 054607 (2022b).
  * Tan  _et al._ [2022] T. H. Tan, A. Mietke, J. Li, Y. Chen, H. Higinbotham, P. J. Foster, S. Gokhale, J. Dunkel, and N. Fakhri, Odd dynamics of living chiral crystals, Nature 607, 287 (2022).
  * Soni  _et al._ [2019] V. Soni, E. S. Bililign, S. Magkiriadou, S. Sacanna, D. Bartolo, M. J. Shelley, and W. T. M. Irvine, The odd free surface flows of a colloidal chiral fluid, Nat. Phys. 15, 1188 (2019).
  * Vega Reyes  _et al._ [2022] F. Vega Reyes, M. A. López-Castaño, and Á. Rodríguez-Rivas, Diffusive regimes in a two-dimensional chiral fluid, Commun. Phys. 5, 1 (2022).
  * Bililign  _et al._ [2022] E. S. Bililign, F. Balboa Usabiaga, Y. A. Ganan, A. Poncet, V. Soni, S. Magkiriadou, M. J. Shelley, D. Bartolo, and W. T. M. Irvine, Motile dislocations knead odd crystals into whorls, Nat. Phys. 18, 212 (2022).
  * Czopnik and Garbaczewski [2001] R. Czopnik and P. Garbaczewski, Brownian motion in a magnetic field, Phys. Rev. E 63, 021105 (2001).
  * Chun  _et al._ [2018] H.-M. Chun, X. Durang, and J. D. Noh, Emergence of nonwhite noise in langevin dynamics with magnetic lorentz force, Phys. Rev. E 97, 032117 (2018).
  * Vuijk  _et al._ [2019] H. D. Vuijk, J. M. Brader, and A. Sharma, Anomalous fluxes in overdamped brownian dynamics with lorentz force, J. Stat. Mech. Theory Exp. 2019, 063203 (2019).
  * Abdoli  _et al._ [2020a] I. Abdoli, H. D. Vuijk, J.-U. Sommer, J. M. Brader, and A. Sharma, Nondiffusive fluxes in a brownian system with lorentz force, Phys. Rev. E 101, 012120 (2020a).
  * Abdoli  _et al._ [2020b] I. Abdoli, E. Kalz, H. D. Vuijk, R. Wittmann, J.-U. Sommer, J. M. Brader, and A. Sharma, Correlations in multithermostat brownian systems with lorentz force, New J. Phys. 22, 093057 (2020b).
  * Shinde  _et al._ [2022] R. Shinde, J.-U. Sommer, H. Löwen, and A. Sharma, Strongly enhanced dynamics of a charged rouse dimer by an external magnetic field, PNAS Nexus 1, pgac119 (2022).
  * Schütte  _et al._ [2014] C. Schütte, J. Iwasaki, A. Rosch, and N. Nagaosa, Inertia, diffusion, and dynamics of a driven skyrmion, Phys. Rev. B 90, 174434 (2014).
  * Troncoso and Núñez [2014] R. E. Troncoso and Á. S. Núñez, Brownian motion of massive skyrmions in magnetic thin films, Ann. Phys. 351, 850 (2014).
  * Wiesendanger [2016] R. Wiesendanger, Nanoscale magnetic skyrmions in metallic films and multilayers: a new twist for spintronics, Nat. Rev. Mater. 1, 1 (2016).
  * Fert  _et al._ [2017] A. Fert, N. Reyren, and V. Cros, Magnetic skyrmions: advances in physics and potential applications, Nat. Rev. Mat. 2, 1 (2017).
  * Büttner  _et al._ [2018] F. Büttner, I. Lemesh, and G. S. D. Beach, Theory of isolated magnetic skyrmions: From fundamentals to room temperature applications, Sci. Rep. 8, 4464 (2018).
  * Weißenhofer  _et al._ [2021] M. Weißenhofer, L. Rózsa, and U. Nowak, Skyrmion dynamics at finite temperatures: Beyond thiele’s equation, Phys. Rev. Lett. 127, 047203 (2021).
  * [42] See Supplemental Material at $***link***$ in which we provide detailed calculations. The two-particle Smoluchowski equation for interacting off-diffusive hard disks in two dimensions is solved, we establish the connection from the velocity to the force autocorrelation function and explicitly solve for the elements of the force autocorrelation tensor in the Laplace domain.
  * Kümmel  _et al._ [2013] F. Kümmel, B. Ten Hagen, R. Wittkowski, I. Buttinoni, R. Eichhorn, G. Volpe, H. Löwen, and C. Bechinger, Circular motion of asymmetric self-propelling particles, Phys. Rev. Lett. 110, 198302 (2013).
  * Van Teeffelen and Löwen [2008] S. Van Teeffelen and H. Löwen, Dynamics of a brownian circle swimmer, Phys. Rev. E 78, 020101 (2008).
  * Vuijk  _et al._ [2022] H. D. Vuijk, S. Klempahn, H. Merlitz, J.-U. Sommer, and A. Sharma, Active colloidal molecules in activity gradients, Phys. Rev. E 106, 014617 (2022).
  * Muzzeddu  _et al._ [2022] P. L. Muzzeddu, H. D. Vuijk, H. Löwen, J.-U. Sommer, and A. Sharma, Active chiral molecules in activity gradients, J. Chem. Phys. 157, 134902 (2022).
  * Reichhardt and Reichhardt [2022] C. J. O. Reichhardt and C. Reichhardt, Active rheology in odd-viscosity systems, Europhys. Lett. 137, 66004 (2022).
  * Kählert  _et al._ [2012] H. Kählert, J. Carstensen, M. Bonitz, H. Löwen, F. Greiner, and A. Piel, Magnetizing a complex plasma without a magnetic field, Phys. Rev. Lett. 109, 155003 (2012).
  * Abdoli  _et al._ [2020c] I. Abdoli, H. D. Vuijk, R. Wittmann, J.-U. Sommer, J. M. Brader, and A. Sharma, Stationary state in brownian systems with lorentz force, Phys. Rev. Res. 2, 023381 (2020c).
  * Dhont [1996] J. K. G. Dhont,  _An Introduction to Dynamics of Colloids_ (Elsevier, 1996).
  * Franosch  _et al._ [2010] T. Franosch, F. Höfling, T. Bauer, and E. Frey, Persistent memory for a brownian walker in a random array of obstacles, Chem. Phys. 375, 540 (2010).
  * Hull and Froese [1955] T. Hull and C. Froese, Asymptotic behaviour of the inverse of a laplace transform, Can. J. Math. 7, 116 (1955).
  * Abramowitz and Stegun [1968] M. Abramowitz and I. A. Stegun,  _Handbook of mathematical functions with formulas, graphs, and mathematical tables_ , Vol. 55 (US Government printing office, 1968).
  * Ernst and Weyland [1971] M. Ernst and A. Weyland, Long time behaviour of the velocity auto-correlation function in a lorentz gas, Phys. Lett. A 34, 39 (1971).
  * van Beijeren [1982] H. van Beijeren, Transport properties of stochastic lorentz models, Rev. Mod. Phys. 54, 195 (1982).
  * Alder and Alley [1978] B. Alder and W. Alley, Long-time correlation effects on displacement distributions, J. Stat. Phys. 19, 341 (1978).
  * Lewis and Tjon [1978] J. C. Lewis and J. Tjon, Evidence for slowly-decaying tails in the velocity autocorrelation function of a two-dimensional lorentz gas, Phys. Lett. A 66, 349 (1978).
  * Frenkel [1987] D. Frenkel, Velocity auto-correlation functions in a 2d lattice lorentz gas: Comparison of theory and computer simulation, Phys. Lett. A 121, 385 (1987).
  * Höfling and Franosch [2007] F. Höfling and T. Franosch, Crossover in the slow decay of dynamic correlations in the lorentz model, Phys. Rev. Lett. 98, 140601 (2007).
  * Jacobs and Harris [1977] G. Jacobs and S. Harris, Macromolecular self-diffusion and momentum autocorrelation functions in dilute solutions, J. Chem. Phys. 67, 5655 (1977).
  * Ackerson and Fleishman [1982] B. J. Ackerson and L. Fleishman, Correlations for dilute hard core suspensions, J. Chem. Phys. 76, 2675 (1982).
  * Felderhof and Jones [1983a] B. Felderhof and R. Jones, Cluster expansion of the diffusion kernel of a suspension of interacting brownian particles, Phys. A: Stat. Mech. 121, 329 (1983a).
  * Felderhof and Jones [1983b] B. U. Felderhof and R. B. Jones, Diffusion in hard sphere suspensions, Phys. A: Stat. Mech. 122, 89 (1983b).
  * Yasuda  _et al._ [2022] K. Yasuda, K. Ishimoto, A. Kobayashi, L.-S. Lin, I. Sou, Y. Hosaka, and S. Komura, Time-correlation functions for odd langevin systems, J. Chem. Phys. 157, 095101 (2022).
  * Kalz [2022] E. Kalz,  _Diffusion under the Effect of Lorentz Force_ (Springer Spektrum Wiesbaden, 2022).
  * Avron [1998] J. E. Avron, Odd viscosity, J. Stat. Phys. 92, 543 (1998).


# Fast and accurate method for computing non-smooth solutions to constrained
control problems

Lucian Nita1, Eduardo M. G. Vila1, Marta A. Zagorowska1, Eric C. Kerrigan1,
Yuanbo Nie2, Ian McInerney1, Paola Falugi1 *This work has received funding
from the Engineering and Physical Sciences Research Council under a Doctoral
Training Grant (reference number: EP/T51780X/1), the Active Building Centre
project (reference number: EP/V012053/1) and the Centre for Doctoral Training
in High Performance Embedded and Distributed Systems (HiPEDS, Grant Reference
EP/L016796/1).1 Department of Electrical & Electronic Engineering, Imperial
College London, SW7 2AZ London, UK, email: {lucian.nita16, emg216,
m.zagorowska, e.kerrigan, i.mcinerney17<EMAIL_ADDRESS>Department
of Automatic Control and Systems Engineering, University of Sheffield, S1 3JD
Sheffield, UK, email<EMAIL_ADDRESS>
###### Abstract

Introducing flexibility in the time-discretisation mesh can improve
convergence and computational time when solving differential equations
numerically, particularly when the solutions are discontinuous, as commonly
found in control problems with constraints. State-of-the-art methods use fixed
mesh schemes, which cannot achieve superlinear convergence in the presence of
non-smooth solutions. In this paper, we propose using a flexible mesh in an
integrated residual method. The locations of the mesh nodes are introduced as
decision variables, and constraints are added to set upper and lower bounds on
the size of the mesh intervals. We compare our approach to a uniform fixed
mesh on a real-world satellite reorientation example. This example
demonstrates that the flexible mesh enables the solver to automatically locate
the discontinuities in the solution, has superlinear convergence and faster
solve times, while achieving the same accuracy as a fixed mesh.

##  I Introduction

The differential equations that model a physical system fully describe the
dynamical behavior of the system and enable development of control strategies
to make the system achieve a desired behavior. Due to the complexity of
physical systems, the differential equations describing the dynamics are often
impossible to solve analytically, necessitating the use of numerical solvers.
The challenge in numerically solving complex differential equations lies in
accurately capturing the nature of the solution while using a discrete number
of time points, called the time mesh. In this paper, we propose a new method
for solving differential equations subject to inequality constraints to a
user-specified accuracy by using a flexible time mesh with an integrated
residual method.

Differential equations, together with initial values for the states, form an
Initial Value Problem (IVP), which can be efficiently solved using explicit
[1] or implicit methods [6]. When the differential equations contain boundary
conditions, they are called Boundary Value Problems (BVPs), and are
customarily solved using collocation methods that compute a polynomial
approximation to the solution [12]. Recently, the authors of [10] proposed
revisiting the use of integrated residual methods for solving dynamic
optimization problems, and they highlighted the advantages an integrated
residual method has when used to solve complex ordinary differential equations
in optimal control.

The quality of the numerical solution of differential equations depends on the
refinement strategy the solver uses to improve the accuracy of the numerical
solution. There are three main refinement methods used for improving the
solution accuracy [2]: (i) h-methods that focus on adjusting the size of the
time intervals (i.e. the time mesh), (ii) p-methods that focus on adjusting
the polynomial degree, and (iii) hp-methods that attempt simultaneous
refinement of the time mesh and the degree of the polynomial. State-of-the-art
refinement methods are mainly based on adaptive time mesh refinement [8, 11],
but most of the existing approaches are unable to achieve superlinear
convergence when discontinuities are present in the solution, and in some
cases the numerical solvers may fail to converge to a solution.

The main contribution of this paper is a new method for solving differential
equations and control problems. The method is based on an integrated residual
method with a flexible time mesh, which allows one to accurately capture
potential discontinuities in the solution. Our approach introduces the time-
mesh nodes as variables in an optimisation problem while numerically
integrating the residual function over each time interval. Using a least
squares approach to compute the residuals enables us to solve a wide range of
problems including ordinary differential equations (ODEs) and differential
algebraic equations (DAEs), BVPs, and consistently over-determined systems.
This new method can also be used to handle the dynamic constraints in optimal
control problems (OCPs), which often have discontinuous solutions.
Additionally, the proposed method can handle differential inclusions or
discontinuities that may appear in the dynamics. The proposed flexible time
mesh provides better performance than a fixed time mesh, with a numerical
example showing superlinear convergence and an improved solution accuracy,
while also decreasing the computational time.

In Section II, we present an overview of the integrated residual method for
solving differential equations and extend the method to the solution of
control problems. We introduce our flexible time mesh scheme in Section III.
In Section IV, we use our proposed flexible mesh scheme to solve a minimum-
time satellite reorientation problem. Section V concludes the paper and
presents an outline for future works.

##  II Problem Definition

###  II-A Differential Equations and Control Problems

Many control problems can be formulated as finding the solution to one or more
feasibility problems of the form

find $\displaystyle x(\cdot),u(\cdot)$ (1a) s.t. $\displaystyle
F(\dot{x}(t),x(t),u(t),t)=0,\quad\forall t\in\mathcal{T}$ (1b) $\displaystyle
G(\dot{x}(t),x(t),u(t),t)\leq 0,\quad\forall t\in\mathcal{T},$ (1c) where
$\mathcal{T}:=[t_{0},t_{f}]\subset\mathbb{R}$ is the time interval over which
the problem is defined, $x:\mathbb{R}\rightarrow\mathbb{R}^{N_{x}}$ are the
state differential variables and are forced to be continuous,
$\dot{x}:\mathbb{R}\rightarrow\mathbb{R}^{N_{x}}$ are the time derivatives of
the state, $u:\mathbb{R}\rightarrow\mathbb{R}^{N_{u}}$ are the algebraic
variables, which can model the control inputs, and are allowed to be
discontinuous. The function
$F:\mathbb{R}^{N_{x}}\times\mathbb{R}^{N_{x}}\times\mathbb{R}^{N_{u}}\times\mathbb{R}\rightarrow\mathbb{R}^{N_{F}}$,
which is typically called the dynamics function, defines a set of $N_{F}$
equality constraints that the controlled system have to satisfy and
$G:\mathbb{R}^{N_{x}}\times\mathbb{R}^{N_{x}}\times\mathbb{R}^{N_{u}}\times\mathbb{R}\rightarrow\mathbb{R}^{N_{g}}$
defines $N_{g}$ path inequality constraints. Moreover, the problem may include
boundary constraints $\displaystyle\Psi_{E}(x(t_{0}),x(t_{f}),t_{0},t_{f})$
$\displaystyle=0,$ (1d) $\displaystyle\Psi_{I}(x(t_{0}),x(t_{f}),t_{0},t_{f})$
$\displaystyle\leq 0,$ (1e) where
$\Psi_{E}:\mathbb{R}^{N_{x}}\times\mathbb{R}^{N_{x}}\times\mathbb{R}\times\mathbb{R}\rightarrow\mathbb{R}^{N_{E}}$
are the boundary equality constraints, and
$\Psi_{I}:\mathbb{R}^{N_{x}}\times\mathbb{R}^{N_{x}}\times\mathbb{R}\times\mathbb{R}\rightarrow\mathbb{R}^{N_{I}}$
are the boundary inequality constraints.

Differential equations are usually included in problem (1) as equality
constraints in (1b). Note that solutions to (1) are _non-unique_ and _non-
smooth_ , in general, even if $F$, $G$ $\Psi_{E}$ and $\Psi_{I}$ are all
smooth functions.

Of particular interest is to note that problems of the form (1) need to be
solved in certain classes of direct methods for solving optimal control
problems [10].

###  II-B Discretisation and Numerical Solution

The problem introduced in Section II-A is an infinite-dimensional problem
because of the dependence on time. To numerically solve an infinite-
dimensional feasibility problem using finite-dimensional optimisation methods,
the problem is discretized. The goal of the resulting problem is to find
approximating functions $\tilde{x}:\mathbb{R}\rightarrow\mathbb{R}^{N_{x}}$and
$\tilde{u}:\mathbb{R}\rightarrow\mathbb{R}^{N_{u}}$, such that the integrated
square of the residual

$\epsilon_{R}:=\frac{1}{(t_{f}-t_{0})\cdot
N_{F}}\int_{t_{0}}^{t_{f}}{\left\lVert
F(\dot{\tilde{x}}(t),\tilde{x}(t),\tilde{u}(t),t)\right\rVert^{2}_{2}}\mathop{}\\!\mathrm{d}t$
(2)

is minimised. The residual
$\|F(\dot{\tilde{x}}(t),\tilde{x}(t),\tilde{u}(t),t)\|^{2}_{2}$ captures how
close the numerical solution $(\tilde{x}(t),\tilde{u}(t))$ is to an exact
solution $(x(t),u(t))$ at each time instant $t$, since (1b) becomes zero at
the exact solution. A scaling factor is added in order to average the residual
over the domain $\mathcal{T}$ and over all components of the dynamics
equations. A necessary and sufficient condition for achieving convergence of
$(\tilde{x}(\cdot),\tilde{u}(\cdot))$ to an exact solution
$(x(\cdot),u(\cdot))$ is to make the integrated residual go to zero. Since the
exact solution $(x(\cdot),u(\cdot))$ is non-unique, in general, defining an
error metric based on the exact solution (such as
$\|\tilde{x}(t)-x(t)\|^{2}_{2}$ for example) would not be practical.

The approximating functions $\tilde{x}$ and $\tilde{u}$ are typically
represented using polynomial basis functions [4, Sect. 1.17.1]. Finding these
approximating functions consists of determining a finite number of
coefficients for the polynomials. In order to numerically enforce constraints
in (1b), we insert the approximating functions $\tilde{x},\tilde{u}$ into
function $F$. In the least squares approach we are proposing here, $F(\cdot)$
is replaced with a function of the coefficients of the polynomials. A
procedure to obtain the approximating function and approximate solutions to
the differential equations will be discussed in Section II-C.

The choice of the basis functions for $(\tilde{x}(\cdot),\tilde{u}(\cdot))$
determines the possible solution space, which implies that the exact solution
$(x(\cdot),u(\cdot))$ may not be representable in that solution space. This
means the integrated residual $\epsilon_{R}$ may be non-zero in general, and
the only case when $\epsilon_{R}$ can be zero is if the span of the chosen
basis functions contains an element of the solution set. When basis functions
span the solution set, it is then possible to numerically determine the exact
trajectories.

To reduce the approximation error in the numerical computation of (2) while
maintaining numerical stability, the interval $[t_{0},t_{f}]$ is split into
$N$ subintervals

$[t_{i},t_{i+1}]=:\mathcal{T}_{i}\subset\mathcal{T}\quad\forall
i\in\\{0,1,\dots,N-1\\}$ (3)

delimited by the time-mesh points $t_{i}$ with $t_{N}=t_{f}$. A diagram
showing the subintervals and the basis functions can be seen in Figure 1, with
each subinterval $\mathcal{T}_{i}$ having a component of $\tilde{x}$ given by
the basis function $\tilde{\chi}_{i}$ with 5 support points
$\tau_{i}^{j},j\in\\{0,1,2,3,4\\}$.

Figure 1: Approximation method and notation used in constructing a numerical
solution to problem (1).

In our implementation, we use a piecewise polynomial of degree $a$ chosen from
the Lagrange basis

$\displaystyle\tilde{x}(t)$ $\displaystyle=\tilde{\chi}_{i}(t)\quad\forall
t\in[t_{i},t_{i+1}]$ (4) $\displaystyle\tilde{\chi}_{i}(t)$
$\displaystyle\coloneqq\sum_{j=0}^{a}{s_{i}^{j}\cdot\prod_{\begin{subarray}{c}k=0\\\
k\neq
j\end{subarray}}^{a}{\frac{t-\tau_{i}^{k}}{\tau_{i}^{j}-\tau_{i}^{k}}}},$

where $\tau_{i}^{j}\in[t_{i},t_{i+1}],\forall j\in\\{0,\dots,a\\}$ are the
supports for the interpolating polynomial in the $i$-th interval, and
$s_{i}^{j}$ for $j\in\\{0,\dots,a\\}$ are the decision variables that need to
be obtained from the optimisation process for each interval
$i\in\\{0,\dots,N-1\\}$. Since the basis function is an interpolating
polynomial, the $s_{i}^{j}$ are the $\tilde{\chi}_{i}$ function values when
evaluated at a finite number of time points $\tau_{i}^{j}$, giving

$s_{i}^{j}\coloneqq\tilde{\chi}_{i}(\tau_{i}^{j})=\tilde{x}(\tau_{i}^{j}),\
\forall i\in\\{0,\dots,N-1\\},j\in\\{0,\dots,a\\}.$ (5)

To enforce state continuity at mesh nodes, we introduce the constraints

$\displaystyle\tilde{\chi}_{i}(t_{i+1})$
$\displaystyle=\tilde{\chi}_{i+1}(t_{i+1})\quad\forall i\in\\{0,\dots,N-2\\},$
(6a) $\displaystyle\tilde{\chi}_{0}(t_{0})$ $\displaystyle=x_{0},$ (6b)
$\displaystyle\tilde{\chi}_{N-1}(t_{f})$ $\displaystyle=x_{f}.$ (6c)

One way of ensuring the constraints in (6) are satisfied is by placing
internal supports at the interval boundaries $\tau_{i}^{0}=t_{i}$ and
$\tau_{i}^{a}=t_{i+1}\quad\forall i\in\\{0,\dots,N-1\\}$. With the internal
supports at the interval boundaries, the same decision variables
$s_{i}^{a}=s_{i+1}^{0}$ can be used to represent both
$\tilde{\chi}_{i}(t_{i+1})$ and $\tilde{\chi}_{i+1}(t_{i+1}),\forall
i\in\\{0,\dots,N-2\\}$.

A similar strategy is used for transforming a control input $u(\cdot)$ into a
piecewise polynomial $\tilde{u}(\cdot)$ based on a finite number of decision
variables $c_{i}^{j}$ where $i$ is the index of the subinterval and
$j\in\\{0,1,...,b\\}$ is the $j^{\text{th}}$ polynomial coefficient. Note that
the degree of polynomial approximation may be different for the states and
inputs ($a\neq b$), so the supports may not overlap. In this case, we will
evaluate all the functions at the supports used by the highest degree
polynomial approximation. For the functions which do not have a support
present at the point of evaluation, we will use the interpolated value.
Finally, $\dot{\tilde{x}}(\cdot)$ is obtained by differentiating the state
approximation function with respect to time.

A typical design choice is to fix the values $t_{i}$ in the time-mesh before
solving problem (1) to create a fixed mesh. The most common approach is to
place mesh nodes $t_{i}$ at equidistant locations to create a uniform mesh.
However, there are also other options for selecting the predefined position of
the nodes [5]. A uniform time-mesh is given by

$t_{i}=t_{0}+i\cdot\frac{t_{f}-t_{0}}{N},$ (7)

where $t_{i}$ are the time-mesh nodes and $i\in\\{0,\dots,N\\}$ is the time-
mesh node index.

The choice of the number of intervals for a fixed mesh depends on the required
solution accuracy with respect to the error metric $\epsilon_{R}$ from (2). An
in-depth discussion on the influence of $N$ on the error metric can be found
in Section III.

###  II-C Integrated Residual Methods

Integrated residual methods (IRMs) can be used to solve differential equations
of the form described in (1b). Additionally, they can be used to solve optimal
control problems with inequality constraints using the method in [10].
Classical collocation methods enforce constraint (1b) only at a finite number
of points, called collocation points. A drawback of the collocation approach
is that one does not have direct control over the numerical error encountered
in-between the collocation points. In contrast, IRMs aim to enforce constraint
(1b) by directly minimising the integrated residual over the entire desired
time-span $\mathcal{T}$.

To efficiently solve feasibility and control problems with discontinuous
solutions, which are otherwise difficult to tackle, we will use an integrated
residual method for solving differential equations that is similar to [7] and
[9]. The idea of this approach is to transform the feasibility problem (1)
into a constrained minimisation problem of the form

$\displaystyle\min_{\tilde{x}(\cdot),\tilde{u}(\cdot)}\quad$
$\displaystyle\epsilon_{R}$ (8a) s.t. $\displaystyle
G(\dot{\tilde{x}}(t),\tilde{x}(t),\tilde{u}(t),t)\leq 0,$
$\displaystyle\forall t\in\mathbf{\tau}$ (8b)
$\displaystyle\Psi_{E}(\tilde{x}(t_{0}),\tilde{x}(t_{f}),t_{0},t_{f})=0,$ (8c)
$\displaystyle\Psi_{I}(\tilde{x}(t_{0}),\tilde{x}(t_{f}),t_{0},t_{f})\leq 0,$
(8d)

with the integrated residual (2) as the objective to be minimised. To keep the
implementation simple, path inequality constraints are only enforced at the
support time points

$\displaystyle\mathbf{\tau}:=\\{\tau_{i}^{j}\mid i=0,1,\dots N-1;j=0,1,\dots
a\\}.$ (9)

In practice, the integrals in (2) and (8a) are solved numerically using
quadrature rules. In our implementation we used Gaussian quadrature with $Q$
quadrature points [13]. The integral in (2) is computed numerically and since
$F$ is a general nonlinear function, this will introduce an additional
numerical error, namely the _quadrature error_

$\epsilon_{Q}\coloneqq\Big{|}\epsilon_{R}-\sum_{i=0}^{N-1}{\sum_{k=1}^{Q}{w_{i}^{k}\cdot\left\lVert
F(\dot{\tilde{x}}(\rho_{i}^{k}),\tilde{x}(\rho_{i}^{k}),\tilde{u}(\rho_{i}^{k}),\rho_{i}^{k})\right\rVert^{2}_{2}\Big{|}}}$
(10)

where $w_{i}^{k}$ for $k\in\\{1,\dots,Q\\}$ are the $Q$ Gaussian quadrature
weights associated with the integration interval $[t_{i},t_{i+1}]$ and
$\rho_{i}^{k}$ are the quadrature nodes for the interval $[t_{i},t_{i+1}]$.
After introducing the quadrature rule, the discretisation method from Section
II-B can be used to turn problem (8) into the discrete formulation

$\displaystyle\min_{\mathbf{s},\mathbf{c}}$
$\displaystyle\sum_{i=0}^{N-1}{\sum_{k=1}^{Q}{w_{i}^{k}\cdot\left\lVert
F(\dot{\tilde{x}}(\rho_{i}^{k}),\tilde{x}(\rho_{i}^{k}),\tilde{u}(\rho_{i}^{k}),\rho_{i}^{k})\right\rVert^{2}_{2}}}$
(11) s.t.
$\displaystyle\eqref{eq:constra},\eqref{eq:varplsq},\eqref{eq:fixedplsq},\eqref{eq:infixedplsq},$

where $\mathbf{s}\in\mathbb{R}^{N(a+1)}$ and
$\mathbf{c}\in\mathbb{R}^{N(b+1)}$ are vectors consisting of all the
polynomial coefficients $s_{i}^{j},c_{i}^{j}$. Problem (11) can now be solved
directly using available NLP solvers.

To make sure the quadrature error is negligible, we need to check that the
number of quadrature points $Q$ is sufficiently high. This can be done after
the original solve by recomputing $\epsilon_{R}$ (to ensure the numerical
integration has converged) using a larger $Q$ than was used to solve (11)
(e.g. with $2Q$ quadrature points) and evaluating $\epsilon_{Q}$ according to
(10). If the difference between the converged $\epsilon_{R}$ and the objective
value obtained from optimization problem (11) is above a certain tolerance
$\varepsilon_{\text{quad,tol}}$ (i.e.
$\epsilon_{Q}>\varepsilon_{\text{quad,tol}}$), this implies that the
quadrature error is significant and the solution needs to be recomputed using
more quadrature points $Q$.

##  III Improving Accuracy Through Mesh Design Updates

###  III-A Flexible mesh

When selecting a suitable number of mesh intervals $N$, one should recall that
the solution $(x(\cdot),u(\cdot))$ is in a particular space, which may be
different from the one spanned by the selected basis functions of the
discretisation. Thus, the residual for the discretisation might be non-zero.

In order to improve the solution accuracy, the conventional mesh refinement
process relies on increasing the number of time nodes or changing the
polynomial degree according to a predefined algorithm based on error metrics
[4]. One popular strategy is to start with a coarse mesh (small $N$), evaluate
the solution accuracy with respect to the error metric $\epsilon_{R}$ while
ensuring the quadrature error $\epsilon_{Q}$ remains below a threshold
$\varepsilon_{\text{quad,tol}}$, and then recompute a new solution for a finer
mesh (large $N$) if the error is above a certain tolerance level.

For many classes of problems, mesh refinement can be ineffective and
inefficient. Consider for example a problem with a discontinuous solution with
discontinuities at unknown locations. Unless it happens by chance for a node
to be placed exactly at a discontinuity, we can consider without loss of
generality that the discontinuity lies inside the interval $(t_{i},t_{i+1})$.
This means that one is trying to approximate a discontinuous function
$x(\cdot)$ by a continuous function $\tilde{\chi}_{i}(\cdot)$. Increasing the
polynomial degree will not help, since a ringing phenomenon will start to
occur. Equivalently, increasing the number of mesh nodes will not improve the
solution for the interior of the interval where the discontinuity is located.

To automate this process of designing the time mesh and capturing
discontinuities and regions of high gradients more accurately, as well as to
improve the overall solution for a given number $N$ of mesh intervals, we
propose to include the mesh nodes $t_{i}$ as decision variables in our
optimization problem.

To ensure the quadrature error $\epsilon_{Q}$ remains insignificant as the
mesh nodes move, we introduce a rule restricting the allowable change in the
individual mesh interval length. Consider a flexibility parameter
$\phi\in[0,1)$, then the constraints

$\displaystyle t_{i+1}-t_{i}$
$\displaystyle\leq(1+\phi)\cdot\frac{t_{f}-t_{0}}{N},\ \forall
i\in\\{1,\dots,N-1\\}$ (12a) $\displaystyle t_{i+1}-t_{i}$
$\displaystyle\geq(1-\phi)\cdot\frac{t_{f}-t_{0}}{N},\ \forall
i\in\\{1,\dots,N-1\\}$ (12b)

will be included in the optimisation problem. It follows that if
$\phi\in(0,1)$ then $t_{0}<t_{1}<\dots<t_{N-1}<t_{N}$ and hence no overlaps
are possible.

Hence, the optimal node locations can be determined automatically by solving
the following optimisation problem

$\displaystyle\min_{\mathbf{s},\mathbf{c},\mathbf{t}}$
$\displaystyle\sum_{i=0}^{N-1}{\sum_{k=1}^{Q}{w_{i}^{k}\cdot\left\lVert
F(\dot{\tilde{x}}(\rho_{i}^{k}),\tilde{x}(\rho_{i}^{k}),\tilde{u}(\rho_{i}^{k}),\rho_{i}^{k})\right\rVert^{2}_{2}}}$
(13) s.t.
$\displaystyle\eqref{eq:constra},\eqref{eq:varplsq},\eqref{eq:fixedplsq},\eqref{eq:infixedplsq},\eqref{eq:quasi_mes}$

where $\mathbf{t}\in\mathbb{R}^{N-1}$ is the vector of mesh points
$t_{i}\quad\forall i\in\\{1,\dots,N-1\\}$ where the initial and final times
($t_{0}$ and $t_{f}$) are excluded, since they are fixed. Note that quadrature
points $\rho_{i}^{k}$ and internal supports $\tau_{i}^{j}$ become a function
of $t_{i}$ and $t_{i+1}$, hence they are shifted and scaled versions of the
original values.

###  III-B Error and Performance Considerations

When implementing the flexible mesh together with the integrated residual
approach, we ensure the accuracy for the trajectories $\tilde{x}(\cdot)$ and
$\tilde{u}(\cdot)$ in-between mesh nodes. Compared to classical time-marching
schemes (shooting methods) or point-wise residual minimisation (collocation),
integrated residual methods have the benefit of producing a solution with a
more uniform error over the whole time domain.

When comparing the problem formulations in (11) and (13), note that the only
difference is represented by the variables we are optimising over. Observe
also that $\phi=0$ leads to the fixed uniform mesh solution, i.e.
$t_{i+1}-t_{i}=(t_{f}-t_{0})/N$. As a result, problem (13) includes the fixed
uniform mesh problem (11). By enlarging the solution space, successfully
solving (13) with a flexible mesh should provide a solution with a residual
error that is not worse than the one obtained using a fixed mesh in (11),
provided the quadrature error is sufficiently small.

For a given value for the flexibility parameter $\phi$, an algorithm that
implements (13) simplifies the mesh refinement procedure, since the mesh node
locations are automatically determined and the main parameters that can be
changed are the number of nodes $N$ and the polynomial degree. This increases
the accuracy of the obtained solution $\tilde{x}(\cdot)$ and
$\tilde{u}(\cdot)$. Additionally, the proposed strategy speeds up the
computation of a solution, as demonstrated in Section IV.

##  IV Results

An integrated residual solver was developed in the Julia v1.6 programming
language, implementing fixed and flexible mesh schemes for the solution of
constrained DAEs. Dynamic variables were parameterised by Lagrange polynomials
in the barycentric form [3]. These were discretised across $N$ time intervals,
as illustrated in Figure 1. Residuals were integrated with Gauss-Legendre
quadrature of a sufficiently high order. Derivatives were evaluated using
forward and reverse automatic differentiation. The least squares method
described in Section II-C was applied using Ipopt [14] as the NLP solver
(relative convergence tolerance set to $10^{-8}$). All tests were performed on
an Intel® Core™ i7-1065G7 at 1.3 GHz with 16 GB of RAM.

###  IV-A Satellite Reorientation Example

The NASA X-ray Timing Explorer spacecraft is modelled as a 3D rigid body with
moments of inertia $I_{xx}>I_{yy}>I_{zz}$. The orientation is defined by
quaternions $q=[q_{1},q_{2},q_{3},q_{4}]$, where $\left\lVert
q\right\rVert_{2}=1$ must always hold. The spacecraft must be reoriented
($\Delta q_{1}=150^{\circ}$) in minimum time, subject to the saturation of
control torques $\left\lVert u\right\rVert_{\infty}\leq 50$ N m. To model the
dynamics, the least squares method allows for the DAE formulation to be used
to directly enforce the quaternion magnitude, as opposed to the ODE
formulation, which requires high-accuracy integration, hence we used the DAE
formulation [4, Sect. 6.8].

Despite the fact that the dynamic functions are continuous (albeit nonlinear),
this control problem yields a discontinuous solution with switches at non-
trivial times. Figure 2 shows a solution to (1) when the final time is set to
the optimum $t_{f}=28.630408\,\mathrm{s}$ [4, Sect. 6.8].

Figure 2: Control solution to satellite reorientation using fixed and flexible
meshes with 15 intervals. Dotted vertical lines indicate the location of the
fixed mesh points, blue crosses indicate the location of the flexible mesh
points.

To perform the maneuver in minimum time, control torques 2 and 3 ($u_{2}$ and
$u_{3}$, respectively) introduce coupled rotations, reducing the effective
moment of inertia.

As expected, the flexible mesh is able to capture the discontinuities by
shifting the mesh points to the switch times. On the other hand, the fixed
mesh has to compromise accuracy at discontinuities, with the (accidental)
exception of control torque 1 ($u_{1}$) where the uniform distribution (7) has
a time node coinciding with the discontinuity at $t=t_{f}/2$.

###  IV-B Convergence

To assess the order of convergence for the proposed scheme, the problem was
solved for an increasing number of intervals using both flexible and fixed
mesh points. The initial guess for the NLP time-mesh decision variables was
set to be the uniform mesh, while the state and control decision variables
$\mathbf{s},\mathbf{c}$ were cold started (matrices full of zeros were given
as the initial guess). Figure 3 shows the value of the integrated residuals
plotted against $N$.

Figure 3: Log-log plot of integrated residual for solutions to the satellite
reorientation problem for both mesh schemes.

The slopes of linear fits for $N\geq 5$ report orders of convergence of 0.70
and 5.55 for fixed and flexible meshes, respectively. Lagrange polynomials of
degree 4 were used as basis functions, with a sufficiently high quadrature
order, $Q=7$.

Once the number of intervals is larger than the number of discontinuities, the
flexible scheme is able to locate the position of the discontinuities. Hence,
for a larger number of intervals, the flexible mesh is able to achieve a
higher order of convergence than the fixed mesh.

###  IV-C Computation time

Adding flexible mesh nodes increases the size of the resulting NLP, and
introduces nonlinearities. On the other hand, this reduces the number of NLP
iterations to achieve a given integrated residual. Computation time was
assessed by comparing the solution time averaged with 5 repeated solves.
Figure 4 shows that for low accuracy solutions ($\epsilon_{R}>10^{-5}$), there
is little difference in computational time between the two schemes. However,
the flexible mesh is able to achieve a significantly better accuracy given the
same solve time as a fixed mesh. Conversely, faster solve times can be
achieved for a given accuracy.

Figure 4: Log-log plot of performance against integrated residuals for both
mesh schemes with different number of intervals $N$ between 2 and 24 and
$\phi=0.5$. The solve times are normalised by the highest value.

##  V Conclusions and future works

Finding discontinuous solutions to constrained control problems is a
challenge, especially if a fixed time-mesh is used. In this paper, we proposed
using a flexible mesh in an integrated residual method to capture possible
discontinuities in the solution. By including the mesh points as decision
variables in the integrated residual method, we are able to capture
discontinuities in the solution, while preserving accuracy of the solution.
The proposed scheme was tested on a challenging aerospace satellite
reorientation problem with nonlinear dynamics and a discontinuous solution. In
contrast to using a fixed mesh, the new method with a flexible mesh
successfully identified the location of the discontinuities, provided enough
mesh intervals were chosen. The flexible scheme showed super-linear
convergence once all discontinuities are captured, whereas the fixed mesh
scheme plateaued. Compared to a fixed mesh scheme, the underlying NLP solver
converged to a solution in fewer iterations, resulting in the faster overall
computation of a solution at a higher accuracy, despite the added overhead of
augmenting the underlying optimisation problem.

These preliminary numerical results on order of convergence and computational
performance are encouraging. Future work could include developing formal
proofs on the order of convergence. The use of a flexible mesh should also be
compared numerically against state-of-the-art adaptive mesh refinement
methods, which do not include mesh points as decision variables, but
iteratively add fixed mesh points to intervals with a large residual. Similar
sophisticated mesh refinement schemes could also be explored when using a
flexible mesh.

The flexible mesh method presented in this paper can be adapted to solve
optimal control problems using an iterative algorithm. This can be done, for
example, by solving a sequence of feasibility problems (13) with an increasing
number of mesh intervals $N$ until $\epsilon_{R}$ is below a given threshold.
The optimal control problem is then solved by defining and solving another
suitable optimisation problem as in [10].

## References

  * [1] R. R. Ahmad, N. Yaacob, and A.-H. Mohd Murid.  Explicit methods in solving stiff ordinary differential equations.  International Journal of Computer Mathematics, 81(11):1407–1415, November 2004.
  * [2] I. Babuška and B. Q. Guo.  The h, p and h-p version of the finite element method; basis theory and applications.  Advances in Engineering Software, 15(3):159–174, January 1992.
  * [3] Jean-Paul Berrut and Lloyd N. Trefethen.  Barycentric Lagrange Interpolation.  SIAM Review, 46(3):501–517, January 2004.
  * [4] John T. Betts.  Practical Methods for Optimal Control and Estimation Using Nonlinear Programming, Second Edition.  Advances in Design and Control. Society for Industrial and Applied Mathematics, January 2010.
  * [5] John T. Betts and William P. Huffman.  Mesh refinement in direct transcription methods for optimal control.  Optimal Control Applications and Methods, 19(1):1–21, 1998.
  * [6] Luigi Brugnano and Cecilia Magherini.  Blended implicit methods for solving ODE and DAE problems, and their extension for second-order problems.  Journal of Computational and Applied Mathematics, 205(2):777–790, August 2007.
  * [7] Ernest D. Eason.  A review of least-squares methods for solving partial differential equations.  International Journal for Numerical Methods in Engineering, 10(5):1021–1046, 1976.  _eprint: https://onlinelibrary.wiley.com/doi/pdf/10.1002/nme.1620100505.
  * [8] Fengjin Liu, William W. Hager, and Anil V. Rao.  Adaptive mesh refinement method for optimal control using nonsmoothness detection and mesh size reduction.  Journal of the Franklin Institute, 352(10):4081–4106, October 2015\.
  * [9] Daniele Mortari, Hunter Johnston, and Lidia Smith.  High accuracy least-squares solutions of nonlinear differential equations.  Journal of Computational and Applied Mathematics, 352:293–307, May 2019.
  * [10] Yuanbo Nie and Eric C Kerrigan.  Solving dynamic optimization problems to a specified accuracy: An alternating approach using integrated residuals.  IEEE Transactions on Automatic Control, 2022.
  * [11] Luís Tiago Paiva and Fernando A. C. C. Fontes.  Adaptive time–mesh refinement in optimal control problems with state constraints.  Discrete & Continuous Dynamical Systems, 35(9):4553, 2015.
  * [12] R. D. Russell and L. F. Shampine.  A collocation method for boundary value problems.  Numerische Mathematik, 19(1):1–28, February 1972.
  * [13] S. P. Venkateshan and Prasanna Swaminathan.  Chapter 9 - Numerical Integration.  In S. P. Venkateshan and Prasanna Swaminathan, editors, Computational Methods in Engineering, pages 317–373. Academic Press, Boston, January 2014.
  * [14] Andreas Wächter and Lorenz T. Biegler.  On the implementation of an interior-point filter line-search algorithm for large-scale nonlinear programming.  Mathematical Programming, 106(1):25–57, March 2006.