Shoja et al. Norwegian Sea target-oriented LSRTM

# Target-oriented least-squares reverse-time migration with Marchenko
redatuming and double-focusing: Field data application

Aydin Shoja11footnotemark: 1 Joost van der Neut22footnotemark: 2 and Kees
Wapenaar11footnotemark: 1

###### Abstract

Recently, the focus of reflection seismologists has shifted to applications
where a high-resolution image of the subsurface is required. Least-Squares
Reverse-Time Migration (LSRTM) is a common tool used to compute such images.
Still, its high computational costs have led seismologists to use target-
oriented LSRTM for imaging only a small target of interest within a larger
subsurface block. Redatuming the data to the upper boundary of the target of
interest is one approach to target-oriented LSRTM. Still, many redatuming
methods cannot account for multiple scatterings within the overburden. This
paper presents a target-oriented least-squares reverse time migration
algorithm which integrates Marchenko redatuming and double-focusing. This
special redatuming method accounts for all orders of multiple scattering in
the overburden for target-oriented LSRTM. Additionally, the paper demonstrates
that a double-focusing algorithm can further reduce the size of the data by
reducing both spatial and temporal dimensions. This algorithm is applied to
field data acquired in the Norwegian Sea.

##  1 Introduction

Seismic imaging and inversion are a set of techniques used by geophysicists to
estimate parameters related to wave propagation, such as reflectivity,
velocity, and density, within the Earth’s subsurface. A network of sources and
receivers is positioned on the Earth’s surface to determine these parameters
to produce and record seismic waves. Geophysicists typically assume a model of
the subsurface that separates into a background model ($m_{0}$) for longer
wavelengths and a short-wavelength reflectivity model ($\delta m$) based on a
weak-scattering assumption (Schuster,, 2017; Claerbout,, 1985). The primary
objective of seismic imaging is to generate a structural image of the short-
wavelength reflectivity model ($\delta m$).

Reverse-Time Migration (RTM) is more popular among different imaging
techniques since it can produce higher-resolution images and better handles
complex geological structures (Baysal et al.,, 1983; McMECHAN,, 1983; Zhou et
al.,, 2018; Zhang et al.,, 2019). RTM creates images by cross-correlating the
forward-propagated wavefield and its back-propagated counterpart based on Born
approximation. However, improving the resolution and quality of RTM images is
still possible by inverting them with a least-squares algorithm (Dutta et
al.,, 2017; Tang,, 2009; Liu et al.,, 2016). This process is known as Least-
Squares Reverse-time migration (LSRTM).

However, LSRTM is a computationally expensive algorithm (Dai et al.,, 2012;
Tang,, 2009; Herrmann and Li,, 2012; Farshad and Chauris,, 2021). To reduce
the computational cost of LSRTM, one can reduce the model’s dimensions by
focusing on a small area inside the big block of the subsurface model. To
compute the image of this smaller region, the wavefield on the boundary of
this region is needed. The process of computing the wavefield on the boundary
of this target from surface recorded data is called redatuming (Valenciano et
al.,, 2006; Haffinger et al.,, 2013; Willemsen et al.,, 2016; Yuan et al.,,
2017; Zhao and Sen,, 2018; Guo and Alkhalifah,, 2020; Ravasi et al.,, 2016).
One prominent redatuming technique is Marchenko redatuming (Wapenaar et al.,,
2014, 2021; Diekmann and Vasconcelos,, 2021).

Marchenko redatuming can create virtual receivers on the boundary of the
target of interest while accounting for all orders of internal multiple
scatterings and reflections. Since Marchekno redatuming and Green’s functions
retrieval are powerful tools, researchers use them to address issues in
seismic imaging and inversion (Cui et al.,, 2020; Zhang et al.,, 2019;
Diekmann et al.,, 2023). Moreover, with Marchenko double-focusing, it is
possible to create virtual sources and redatum both sources and receivers to
the boundary of the target (Staring et al.,, 2018; Shoja et al.,, 2023).
Marchenko double-focused wavefields account for all orders of internal
multiples generated inside the overburden, enabling us to create images with
less impact from internal multiples. Moreover, Marchenko double-focusing
compacts the data’s time axis, reducing the data’s size even more.

This paper combines the Marchenko double-focusing and target-oriented LSRTM
algorithm to create high-resolution artifact-free images of a marine data set
from the Vøring region in the Norwegian Sea. First, we review the theory of
target-oriented LSRTM with Marchenko double-focusing. Second, we apply this
algorithm to a marine dataset, and finally, we discuss the results and
conclude the paper.

##  2 Theory

###  2.1 Least-squares reverse-time migration

Dai et al., (2012) show that classical RTM can be derived from the Born
approximation of seismic reflection data. In the Born approximation, the
incident wavefield ($P^{inc}$) can be represented using the background Green’s
function. The reflectivity model is expressed as $\delta
m=(\frac{1}{c^{2}}-\frac{1}{c_{0}^{2}})$ where $c$ represents the modified
velocity. This equation links $\delta m$ to the scattered data ($P^{scat}$)
through a linear equation (Born et al.,, 1999; Schuster,, 2017; van den Berg,,
2021):

$P^{scat}_{pred}(\textbf{x}_{r},\textbf{x}_{s},\delta
m,\omega)=\frac{\omega^{2}}{\rho_{0}}\int_{V}G_{0}(\textbf{x}_{r},\textbf{x},\omega)\delta
m(\textbf{x})G_{0}(\textbf{x},\textbf{x}_{s},\omega)W(\omega)\,d\textbf{x}.$
(1)

The integral in Equation 1 is computed throughout the model’s volume ($V$),
with $\omega$ being the angular frequency, and the subscripts ”$r$” and ”$s$”
indicating the receiver and source, respectively. This equation can be
expressed in an operator format as follows:

$P^{scat}_{pred}(\textbf{x}_{r},\textbf{x}_{s},\delta
m,\omega)=\mathcal{L}\delta m.$ (2)

Here $\mathcal{L}$ is the forward Born operator.

The standard method of reverse-time migration involves obtaining an
approximate reflectivity model by taking the adjoint of $\mathcal{L}$ and
applying it to the observed scattered data:

$\delta m^{mig}(\textbf{x})=\mathcal{L}^{\dagger}P^{scat}_{obs}.$ (3)

Due to the fact that the adjoint of this kernel is merely an approximation of
its inverse, the resolution of the reflectivity model obtained through this
process is limited.

To tackle the problem of limited resolution, scholars have adopted a least-
squares strategy in which the adjoint operator ($\mathcal{L}^{\dagger}$) is
substituted with a damped least-squares solution, as suggested in references
(Marquardt,, 1963; Dai et al.,, 2012; Dutta et al.,, 2017):

$\delta
m^{mig}=[\mathcal{L}^{\dagger}\mathcal{L}+\epsilon]^{-1}\mathcal{L}^{\dagger}P^{scat}_{obs}.$
(4)

Unfortunately, calculating the Hessian matrix
($\mathcal{L}^{\dagger}\mathcal{L}$) and its inverse is computationally
infeasible. As an alternative, an iterative algorithm that minimizes the
L2-norm of the discrepancy between the observed and anticipated data is often
used to update the reflectivity model:

$C(\delta m)=\frac{1}{2}\big{\|}P^{scat}_{pred}(\delta
m)-P^{scat}_{obs}\big{\|}_{2}^{2}.$ (5)

One potential way to tackle this optimization problem is by utilizing a
conjugate gradient algorithm (Nocedal and Wright,, 2006). In least-squares
reverse-time migration, the background velocity model ($c_{0}(\textbf{x}$)) is
not changed, and only the reflectivity model ($\delta m$) is updated,
resulting in the Green’s functions of Equation 1 being calculated only once.
To learn more about least-squares reverse-time migration, please see Schuster,
(2017).

###  2.2 Marchenko redatuming and double-focusing

Marchenko redatuming is an innovative data-driven technique that can recover
the Green’s function above the target area’s surface, including all orders of
multiple-scattered events. This method only requires the reflection response
at the surface and a smooth background velocity model of the overburden
capable of predicting the direct arrival from the surface to the redatuming
level.

To retrieve the Green’s functions at the redatuming level, the coupled
Marchenko-type representations are solved iteratively, as shown in the
following equations (Wapenaar et al.,, 2014):

$G^{-}_{Mar}(\textbf{x}_{v},\textbf{x}_{r},\omega)=\int_{\mathcal{D}_{acq}}R(\textbf{x}_{r},\textbf{x}_{s},\omega)f_{1}^{+}(\textbf{x}_{s},\textbf{x}_{v},\omega)\,d\textbf{x}_{s}-f_{1}^{-}(\textbf{x}_{r},\textbf{x}_{v},\omega),$
(6)

and

$G^{+}_{Mar}(\textbf{x}_{v},\textbf{x}_{r},\omega)=-\int_{\mathcal{D}_{acq}}R(\textbf{x}_{r},\textbf{x}_{s},\omega)f_{1}^{-}(\textbf{x}_{s},\textbf{x}_{v},\omega)^{*}\,d\textbf{x}_{s}+f_{1}^{+}(\textbf{x}_{r},\textbf{x}_{v},\omega)^{*}.$
(7)

In these equations, $\mathcal{D}_{acq}$ represents the acquisition surface
where $\textbf{x}_{s}$ and $\textbf{x}_{r}$ are situated. $G^{-}_{Mar}$ and
$G^{+}_{Mar}$ denote the up-going and down-going components of the Marchenko
redatumed Green’s function, respectively (see Fig.
LABEL:fig:images/Gpluspic.eps and LABEL:fig:images/Gminpic.eps). Additionally,
$f_{1}^{-}(\textbf{x}_{s},\textbf{x}_{v},\omega)$ and
$f_{1}^{+}(\textbf{x}_{s},\textbf{x}_{v},\omega)$ denote the up-going and
down-going parts of the focusing function, respectively, with the subscript
”$v$” denoting a virtual point situated on the redatuming level denoted by
$\mathcal{D}_{tar}$. Furthermore, $R(\textbf{x}_{r},\textbf{x}_{s},\omega)$
refers to the dipole response of the medium at the acquisition surface, and it
is related to the up-going Green’s function ($G^{-}$) via the following
relationship:

$R(\textbf{x}_{r},\textbf{x}_{s},\omega)=\frac{\partial_{3,s}G^{-}(\textbf{x}_{r},\textbf{x}_{s},\omega)}{\frac{1}{2}i\omega\rho(\textbf{x}_{s})},$
(8)

The partial derivative in the downward direction taken at $\textbf{x}_{s}$ is
denoted by $\partial_{3,s}$. It is important to remove horizontally
propagating waves and surface-related multiples before inserting
$R(\textbf{x}_{r},\textbf{x}_{s},\omega)$ into Equations 6 and 7. The detailed
derivation of these integrals and their solution for computing the focusing
functions and Green’s functions can be found in Wapenaar et al., (2014) and
Thorbecke et al., (2017).

The above-mentioned equations correspond to single-sided redatuming. To
perform a double-sided redatuming, a convolution operation on the up-going and
down-going parts of the Marchenko redatumed Green’s function is proposed by
Staring et al., (2018). This operation involves filtering the down-going
focusing function in a multi-dimensional manner:

$G^{-,+}_{df}(\textbf{x}_{v},\textbf{x}^{\prime}_{v},\omega)=\int_{\mathcal{D}_{acq}}G^{-}_{Mar}(\textbf{x}_{v},\textbf{x}_{r},\omega)\mathcal{F}_{1}^{+}(\textbf{x}_{r},\textbf{x}^{\prime}_{v},\omega)\,d\textbf{x}_{r},$
(9)

and

$G^{+,+}_{df}(\textbf{x}_{v},\textbf{x}^{\prime}_{v},\omega)=\int_{\mathcal{D}_{acq}}G^{+}_{Mar}(\textbf{x}_{v},\textbf{x}_{r},\omega)\mathcal{F}_{1}^{+}(\textbf{x}_{r},\textbf{x}^{\prime}_{v},\omega)\,d\textbf{x}_{r},$
(10)

where

$\mathcal{F}_{1}^{+}(\textbf{x}_{r},\textbf{x}^{\prime}_{v},\omega)=\frac{\partial_{3,r}f_{1}^{+}(\textbf{x}_{r},\textbf{x}^{\prime}_{v},\omega)}{\frac{1}{2}i\omega\rho(\textbf{x}_{r})}.$
(11)

Here the vertical derivative is taken with respect to $\textbf{x}_{r}$.
Equations 9 and 10 use superscripts to indicate the direction of propagation
at the receiver and source locations, respectively. The term ”$df$” stands for
”double-focused.” This process is referred to as ”Marchenko double-focusing.”

The Marchenko double-focusing technique yields two Green’s functions, namely a
down-going ($G^{+,+}_{df}$) and an up-going ($G^{-,+}_{df}$) Green’s function
(fig. LABEL:fig:images/Gplusdfpic.eps and LABEL:fig:images/Gmindfpic.eps). The
down-going Green’s function has a band-limited delta function and interactions
between the target and the overburden. $G^{-,+}_{df}$ can be interpreted as
the continuation of propagation of $G^{+,+}_{df}$ through the target and
recording the up-going part of it at the redatuming level. This up-going
wavefield includes interactions between the target and the overburden on the
source side. In contrast, the conventional double-focusing approach involves
using the inverse of the direct arrival of the transmission response of the
overburden instead of the down-going Marchenko focusing function. However,
this approach cannot predict and remove the multiples generated by the
overburden. In the subsequent sections, the term ”double-focusing” is a
general term that refers to both methods, and it is explicitly mentioned where
a distinction between the methods is necessary.

*4images/Gpluspic.eps,images/Gminpic.eps,images/Gplusdfpic.eps,images/Gmindfpic.epswidth=0.4 The Green’s functions resulting from Marchenko redatuming and double-focusing. a) down-going part of Marchenko Green’s function, b) up-going part of Marchenko Green’s function, c) down-going Marchenko double-focused Green’s function, and d) up-going Marchenko double-focused Green’s function.

###  2.3 Target-oriented LSRTM by Marchenko double-focusing

The following integral is the base for target-oriented LSRTM by Marchenko
double-focusing:

$\hat{P}^{scat}_{pred}(\textbf{x}^{\prime}_{vr},\textbf{x}^{\prime}_{vs},\delta
m,\omega)=\frac{\omega^{2}}{\rho_{0}}\int_{\nu}\hat{G}_{0}(\textbf{x}^{\prime}_{vr},\textbf{x},\omega)\delta
m(\textbf{x})P^{inc}_{df}(\textbf{x},\textbf{x}^{\prime}_{vs},\omega)d\textbf{x}.$
(12)

Here, $\nu$ is the target volume, x is a point inside the target, and
$\textbf{x}^{\prime}_{vs}$ and $\textbf{x}^{\prime}_{vr}$ are the virtual
source and virtual receiver locations on the upper boundary of the target,
respectively. Moreover,

$P^{inc}_{df}(\textbf{x},\textbf{x}^{\prime}_{vs},\omega)=\int_{\mathcal{D}_{tar}}\frac{\partial_{3,vs}G_{0}(\textbf{x},\textbf{x}_{vs},\omega)}{\frac{1}{2}i\omega\rho(\textbf{x}_{vs})}G^{+,+}_{df}(\textbf{x}_{vs},\textbf{x}^{\prime}_{vs},\omega)W(\omega)\,d\textbf{x}_{vs},$
(13)

and

$\hat{G}_{0}(\textbf{x}^{\prime}_{vr},\textbf{x},\omega)=\int_{\mathcal{D}_{tar}}\Gamma(\textbf{x}^{\prime}_{vr},\textbf{x}_{vr},\omega)G_{0}(\textbf{x}_{vr},\textbf{x},\omega)\,d\textbf{x}_{vr},$
(14)

where $\Gamma(\textbf{x}^{\prime}_{vr},\textbf{x}_{vr},\omega)$ is a point-
spread function, which acts as a band limitation filter on the predicted data.
For a complete derivation of the above equations and the definition of the
point-spread function
($\Gamma(\textbf{x}^{\prime}_{vr},\textbf{x}_{vr},\omega)$), we refer to Shoja
et al., (2023). Thus, the new cost function is:

$C(\delta m)=\frac{1}{2}\big{\|}\hat{P}^{scat}_{pred}(\delta
m)-\hat{P}^{scat}_{obs}\big{\|}_{2}^{2},$ (15)

which we solve with a conjugate gradient algorithm.

##  3 Field data example

###  3.1 Field data explanation

This part of the paper shows the results of the Marchenko-based target-
oriented LSRTM on a field dataset provided by Equinor, which was acquired in
the Norwegian Sea in 1994. The depth of the water bottom is $1.5km$, which is
deep enough to separate the free-surface multiple reflections from the primary
and internal multiple reflections. The field dataset contains 399 shot gathers
with 180 traces per gather, and the spatial sampling of sources and receivers
is $25m$. The field dataset was processed according to Davydenko and
Verschuur, 2018a  methodology, which involved muting the direct wave,
estimating near-offset traces through the parabolic Radon transform (Kabir and
Verschuur,, 1995), compensating for 3D effects by multiplying with $\sqrt{t}$,
and deconvolving the source wavelet. Source-receiver reciprocity is also
applied to create offsets in the positive direction to prepare the dataset for
the Estimation of Primaries through the Sparse Inversion (EPSI) method to
remove free-surface multiples (van Groenestijn and Verschuur,, 2009). After
source-receiver reciprocity, each gather contains 371 receivers. Table 1 shows
the acquisition parameters of this dataset.

Parameter | Value
---|---
Number of source positions | 399
Source spacing | 25 m
First source position | 5,000 m
Final source position | 14,950
Number of receiver positions per source | 180
Receiver spacing | 25m
Minimum source-receiver offset | 150 m
Maximum source-receiver offset | 4,625 m
Number of time samples | 2001
Sampling rate | 0.004 s
High-cut frequency | 90 Hz
Table 1: Acquisition Parameters for the dataset

Figure 1 shows the surface reflection response after prepossessing, with a
source located at $\textbf{x}_{s}=(5000m,0m)$. A Ricker wavelet with a
dominant frequency of $30Hz$ is convolved with the data. We choose two
different targets inside the medium.

Figure 1: Reflection response with a source located at
$\textbf{x}_{s}=(5000m,0m)$

###  3.2 LSRTM with double-focusing

####  3.2.1 Target of interest 1

Figure 2 shows the smooth velocity model provided by Equinor for migration.
The red rectangle inside the velocity model indicates the target area, and the
virtual sources and receivers’ positions are at the upper boundary of this
target area.

Figure 2: The smooth velocity model provided by Equinor for migration. The red
rectangle inside the velocity model indicates the first target area. The
virtual sources and receivers’ positions are at the upper boundary of this
target area.

We apply the double-focusing algorithm to the field data for this target. For
this, we define 241 virtual sources and 241 virtual receivers with a spacing
of $12.5m$ at $2500m$ depth extending from $9000m$ to $12000m$ over the upper
boundary of the first target area. The up-going wavefield resulting from
double-focusing is used as input for LSRTM and is called ’observed data’ in
the following. Figure 3 shows the observed, and predicted data, and the
residuals of Marchenko double-focusing target-oriented LSRTM. Moreover, Figure
4 shows the same but for a conventional double-focusing approach. Conventional
means using the inverse of the direct arrival between the target and the
surface as the redatuming operator instead of the Marchenko focusing
functions. The non-physical noises inside the data are caused by imperfect
surface multiple elimination in this part of the data. The computational
advantage of target-oriented LSRTM with double-focused data is twofold. First,
this algorithm reduces the spatial dimension of the problem, and second, it
reduces the temporal dimension of the problem as well. The original recording
time of the data at the surface is 8 seconds, whereas the temporal length of
the double-focused data is 0.5 seconds.

Figure 3: Marchenko double-focused data with a virtual source located at
$\textbf{x}_{vs}=(10500m,2500m)$ and virtual receivers at the same depth as
virtual sources. a) observed data, b) predicted data after 35 iterations of
LSRTM, and c) residuals after 35 iterations of LSRTM. Figure 4: Conventional
double-focused data with a virtual source located at
$\textbf{x}_{vs}=(10500m,2500m)$. a) observed data, b) predicted data after 35
iterations of LSRTM, and c) residuals after 35 iterations of LSRTM.

Figure 5 compares the LSRTM images of using Marchenko and conventional double-
focused data as input. Figure 5 shows some improvements from using Marchenko
double-focused wavefields compared to conventional double-focused ones.

Figure 5: Comparison of images obtained with Marchenko target-oriented LSRTM
(a) and Conventional one (b). Red lines in panel (a) delineate some trends
that are not visible in panel (b), and the black arrows and rectangles in
panel (b) show some events that may be internal multiple reflection artifacts
that are suppressed in panel (a).

Moreover, Figure 6 compares the RTM and LSRTM images of Marchenko double-
focused data as input. The LSRTM algorithm improved the quality of the image.

Figure 6: Comparison of images obtained with Marchenko target-oriented RTM (a)
and LSRTM (b) of the first target.

####  3.2.2 Target of interest 2

Here we choose another target. This target is located between depths of
$2100m$ and $2600m$ and lateral extension from $7000m$ to $10000m$ as shown in
Figure 7. Virtual sources and receivers are located at the upper boundary of
this target area.

Figure 7: The smooth velocity model provided by Equinor for migration. The red
rectangle inside the velocity model indicates the second target area, and the
virtual sources and receivers’ positions are at the upper boundary.

Figure 8 shows the observed, and predicted data, and the residuals of the
Marchenko double-focused approach. Figure 9 shows the same for the
conventional double-focused approach.

Figure 8: Marchenko double-focused data with a virtual source located at
$\textbf{x}_{vs}=(8500m,2100m)$ and virtual receivers at the same depth as
virtual sources. a) observed data, b) predicted data after 35 iterations of
LSRTM, and c) residuals after 35 iterations of LSRTM. Figure 9: Convetional
double-focused data with a virtual source located at
$\textbf{x}_{vs}=(8500m,2100m)$. a) observed data, b) predicted data after 35
iterations of LSRTM, and c) residuals after 35 iterations of LSRTM.

Moreover, Figure 10 shows the LSRTM images of the target-oriented algorithm
with Marchenko and conventional double-focusing. The black arrows and ellipse
indicate the internal multiple reflections that are suppressed by our method.
Figure 11 shows the RTM and LSRTM images of the target-oriented algorithm with
Marchenko double-focusing. The quality and resolution of the image are
increased by the LSRTM algorithm.

Figure 10: Comparison of images obtained with Marchenko target-oriented LSRTM
(a) and Conventional one (b). The black arrows and the ellipse in panel (b)
indicate some of the internal multiple reflection artifacts that are
suppressed in panel (a) Figure 11: Comparison of images obtained with
Marchenko target-oriented RTM (a) and LSRTM (b) of the second target.

##  4 Discussion

In section 2 of this paper, we derive a target-oriented LSRTM algorithm based
on double-focusing that can significantly reduce the dimensions of the
problem, which also reduces the computational costs of the LSRTM algorithm. We
also integrate the Marchenko double-focusing algorithm with our target-
oriented LSRTM algorithm to reduce the artifacts caused by internal multiple
reflections.

To demonstrate the advantages of our proposed algorithm, we applied it to a
dataset acquired by Equinor in the Norwegian Sea in 1994. We chose two
different target zones. Figures 2 and 7 show our targets of interest embedded
in the entire domain of the region. This spatial dimension reduction is to
validate the first advantage we mentioned above. Figures 4a, 5a, 9a, and 10a
show the double-focused data with a recording duration of $0.5s$, whereas the
recording time of the original data is $8s$.

To move forward with our investigation, we showed the imaging results with
double-focusing for both targets. Figure 5 compares the imaging results of the
conventional and Marchenko double-focusing target-oriented LSRTM. The first
panel (fig. 6a) shows the LSRTM result of our proposed algorithm with
Marchenko double-focused data, and the second panel (fig. 6b) shows the LSRTM
results with conventional double-focused data. Comparing these two panels
reveals that using Marchenko wavefields leads to better visualization of true
events and fewer artifacts due to internal multiples delineated by the lines
and arrows in those panels. Moreover, Figure 6 shows the resolution and
quality improvement of target-oriented LSRTM compared to target-oriented RTM
with Marchenko double-focusing.

The same discussion stands for the second target. Figure 10 shows a comparison
between conventional and Marchenko double-focusing target-oriented LSRTM
images where the internal multiple suppression is visible and indicated by
arrows and an ellipse, and Figure 11 shows the RTM and LSRTM images of
Marchenko double-focusing target-oriented LSRTM. The quality and resolution of
the image are increased noticeably. The internal multiple suppression results
of this paper are in agreement with the findings of Davydenko and Verschuur,
2018b  and Ypma and Verschuur, (2013).

##  5 Conclusion

This paper discusses a target-oriented LSRTM algorithm based on double-
focusing. The advantages of this algorithm are: 1) reduction of the spatial
dimensions of the problem by choosing a smaller target of interest, and 2)
reduction of the temporal dimension of the problem by creating both virtual
sources and receivers at the boundary of the target, which leads to lower
computational costs. One can also opt for more sophisticated redatuming
algorithms such as Marchenko redatuming and double-focusing to create virtual
sources and receivers. The advantage of using Marchenko double-focusing
compared to a more conventional redatuming algorithm is the ability to predict
the internal multiple reflections inside the overburden and a reduction of
artifacts due to these multiple reflections.

Present-day seismic imaging and inversion applications need more accurate and
higher-resolution images. Computing higher-resolution images demands
significant amounts of computational power and time. Thus, devising algorithms
that can reduce this computational burden is essential. Our proposed target-
oriented algorithm is not only able to greatly reduce the spatial and temporal
dimensions of the problem but also can reduce the artifacts due to internal
multiple reflections by integrating Marchenko double-focusing with LSRTM
algorithm. Consequently, our algorithm enables us to create higher-resolution
images with fewer artifacts at a lower computational cost.

##  6 ACKNOWLEDGMENTS

The authors would like to thank Eric Verschuur for fruitful discussions and
feedback. This work was supported by the European Union’s Horizon 2020
Research and Innovation Program: European Research Council under Grant 742703.

## References

  * Baysal et al., (1983) Baysal, E., D. D. Kosloff, and J. W. C. Sherwood, 1983, Reverse time migration: GEOPHYSICS, 48, 1514–1524.
  * Born et al., (1999) Born, M., E. Wolf, A. B. Bhatia, P. C. Clemmow, D. Gabor, A. R. Stokes, A. M. Taylor, P. A. Wayman, and W. L. Wilcock, 1999, Principles of optics: Electromagnetic theory of propagation, interference and diffraction of light, 7 ed.: Cambridge University Press.
  * Claerbout, (1985) Claerbout, J. F., 1985, Imaging the earth’s interior: BlackWell Scientific Publications.
  * Cui et al., (2020) Cui, T., J. Rickett, I. Vasconcelos, and B. Veitch, 2020, Target-oriented full-waveform inversion using Marchenko redatumed wavefields: Geophysical Journal International, 223, 792–810.
  * Dai et al., (2012) Dai, W., P. Fowler, and G. T. Schuster, 2012, Multi-source least-squares reverse time migration: Geophysical Prospecting, 60, 681–695.
  * (6) Davydenko, M., and D. J. Verschuur, 2018a, Including and using internal multiples in closed-loop imaging – field data examples: GEOPHYSICS, 83, R297–R305.
  * (7) ——–, 2018b, Including and using internal multiples in closed-loop imaging â€” field data examples: GEOPHYSICS, 83, R297–R305.
  * Diekmann and Vasconcelos, (2021) Diekmann, L., and I. Vasconcelos, 2021, Focusing and Green’s function retrieval in three-dimensional inverse scattering revisited: A single-sided marchenko integral for the full wave field: Phys. Rev. Research, 3, no. 1, 013206.
  * Diekmann et al., (2023) Diekmann, L., I. Vasconcelos, and T. van Leeuwen, 2023, A note on Marchenko-linearised full waveform inversion for imaging: Geophysical Journal International, 234, 228–242.
  * Dutta et al., (2017) Dutta, G., M. Giboli, C. Agut, P. Williamson, and G. T. Schuster, 2017, Least-squares reverse time migration with local Radon-based preconditioning: GEOPHYSICS, 82, S75–S84.
  * Farshad and Chauris, (2021) Farshad, M., and H. Chauris, 2021, Sparsity-promoting multiparameter pseudoinverse Born inversion in acoustic media: GEOPHYSICS, 86, S205–S220.
  * Guo and Alkhalifah, (2020) Guo, Q., and T. Alkhalifah, 2020, Target-oriented waveform redatuming and high-resolution inversion: Role of the overburden: GEOPHYSICS, 85, R525–R536.
  * Haffinger et al., (2013) Haffinger, P., A. Gisolf, and P. M. v. d. Berg, 2013, Towards high resolution quantitative subsurface models by full waveform inversion: Geophysical Journal International, 193, 788–797.
  * Herrmann and Li, (2012) Herrmann, F. J., and X. Li, 2012, Efficient least-squares imaging with sparsity promotion and compressive sensing: Geophysical Prospecting, 60, 696–712.
  * Kabir and Verschuur, (1995) Kabir, M. N., and D. Verschuur, 1995, Restoration of missing offsets by parabolic radon transform1: Geophysical Prospecting, 43, 347–368.
  * Liu et al., (2016) Liu, Y., X. Liu, A. Osen, Y. Shao, H. Hu, and Y. Zheng, 2016, Least-squares reverse time migration using controlled-order multiple reflections: GEOPHYSICS, 81, S347–S357.
  * Marquardt, (1963) Marquardt, D. W., 1963, An algorithm for least-squares estimation of nonlinear parameters: Journal of the Society for Industrial and Applied Mathematics, 11, 431–441.
  * McMECHAN, (1983) McMECHAN, G. A., 1983, Migration by extrapolation of time-dependent boundary values*: Geophysical Prospecting, 31, 413–420.
  * Nocedal and Wright, (2006) Nocedal, J., and S. J. Wright, 2006, Numerical optimization: Springer.
  * Ravasi et al., (2016) Ravasi, M., I. Vasconcelos, A. Kritski, A. Curtis, C. A. d. C. Filho, and G. A. Meles, 2016, Target-oriented Marchenko imaging of a North Sea field: Geophysical Journal International, 205, 99–104.
  * Schuster, (2017) Schuster, G. T., 2017, Seismic inversion: Society of Exploration Geophysicists.
  * Shoja et al., (2023) Shoja, A., J. van der Neut, and K. Wapenaar, 2023, Target-oriented least-squares reverse-time migration using Marchenko double-focusing: reducing the artefacts caused by overburden multiples: Geophysical Journal International, 233, 13–32.
  * Staring et al., (2018) Staring, M., R. Pereira, H. Douma, J. van der Neut, and K. Wapenaar, 2018, Source-receiver Marchenko redatuming on field data using an adaptive double-focusing method: GEOPHYSICS, 83, S579–S590.
  * Tang, (2009) Tang, Y., 2009, Target-oriented wave-equation least-squares migration/inversion with phase-encoded Hessian: GEOPHYSICS, 74, WCA95–WCA107.
  * Thorbecke et al., (2017) Thorbecke, J. W., E. Slob, J. Brackenhoff, J. van der Neut, and K. Wapenaar, 2017, Implementation of the Marchenko method: GEOPHYSICS, 82, WB29–WB45.
  * Valenciano et al., (2006) Valenciano, A. A., B. Biondi, and A. Guitton, 2006, Target-oriented wave-equation inversion: GEOPHYSICS, 71, A35–A38.
  * van den Berg, (2021) van den Berg, P. M., 2021, Acoustic waves, in Forward and Inverse Scattering Algorithms based on Contrast Source Integral Equations: John Wiley & Sons, Ltd, 2, 79–179.
  * van Groenestijn and Verschuur, (2009) van Groenestijn, G. J., and D. J. Verschuur, 2009, Estimating primaries by sparse inversion and application to near-offset data reconstruction: GEOPHYSICS, 74, A23–A28.
  * Wapenaar et al., (2021) Wapenaar, K., J. Brackenhoff, M. Dukalski, G. Meles, C. Reinicke, E. Slob, M. Staring, J. Thorbecke, J. van der Neut, and L. Zhang, 2021, Marchenko redatuming, imaging, and multiple elimination and their mutual relations: GEOPHYSICS, 86, WC117–WC140.
  * Wapenaar et al., (2014) Wapenaar, K., J. Thorbecke, J. van der Neut, F. Broggini, E. Slob, and R. Snieder, 2014, Marchenko imaging: GEOPHYSICS, 79, WA39–WA57.
  * Willemsen et al., (2016) Willemsen, B., A. Malcolm, and W. Lewis, 2016, A numerically exact local solver applied to salt boundary inversion in seismic full-waveform inversion: Geophysical Journal International, 204, 1703–1720.
  * Ypma and Verschuur, (2013) Ypma, F., and D. Verschuur, 2013, Estimating primaries by sparse inversion, a generalized approach: Geophysical Prospecting, 61, 94–108.
  * Yuan et al., (2017) Yuan, S., N. Fuji, S. Singh, and D. Borisov, 2017, Localized time-lapse elastic waveform inversion using wavefield injection and extrapolation: 2-D parametric studies: Geophysical Journal International, 209, 1699–1717.
  * Zhang et al., (2019) Zhang, L., J. Thorbecke, K. Wapenaar, and E. Slob, 2019, Transmission compensated primary reflection retrieval in the data domain and consequences for imaging: GEOPHYSICS, 84, Q27–Q36.
  * Zhao and Sen, (2018) Zhao, Z., and M. K. Sen, 2018, Fast image-domain target-oriented least-squares reverse time migration: GEOPHYSICS, 83, A81–A86.
  * Zhou et al., (2018) Zhou, H.-W., H. Hu, Z. Zou, Y. Wo, and O. Youn, 2018, Reverse time migration: A prospect of seismic imaging methodology: Earth-Science Reviews, 179, 207–227.


# The Lorenz system as a gradient-like system

Jeremy P. Parker<EMAIL_ADDRESS>Division of Mathematics, University
of Dundee, Dundee DD1 4HN, United Kingdom

###### Abstract

We formulate, for continuous-time dynamical systems, a sufficient condition to
be a gradient-like system, i.e. that all bounded trajectories approach
stationary points and therefore that periodic orbits, chaotic attractors, etc.
do not exist. This condition is based upon the existence of an auxiliary
function defined over the state space of the system, in a way analogous to a
Lyapunov function for the stability of an equilibrium. For polynomial systems,
Lyapunov functions can be found computationally by using sum-of-squares
optimisation. We demonstrate this method by finding such an auxiliary function
for the Lorenz system. We are able to show that the system is gradient-like
for $0\leq\rho\leq 12$ when $\sigma=10$ and $\beta=8/3$, significantly
extending previous results. The results are rigorously validated by a novel
procedure: First, an approximate numerical solution is found using finite-
precision floating-point sum-of-squares optimisation. We then prove that there
exists an exact solution close to this using interval arithmetic.

##  I Introduction

Nonlinear dynamical systems exhibit incredibly rich behaviour. A simple set of
precisely deterministic differential equations can give rise to completely
unpredictable chaos, while in other cases similar equations give no
interesting dynamics at all, with all trajectories converging on the same
equilibrium. One approach to understanding such systems is to find
topologically simple trajectories, such as stationary points and periodic
orbits, with the idea that more complex trajectories must approximately follow
these when close. Even these ‘simple’ trajectories, however, can be very
difficult to detect, and the fact that they have not been found should not be
taken as evidence that they do not exist.

Finding stationary points in systems of ordinary differential equations is
straightforward, and for the polynomial systems studied in this work, this
reduces to the well-studied problem of finding roots of polynomials. Often,
one can be certain that they have found all the stationary points of a system.
Finding periodic orbits, however, is a much harder problem. Indeed, the famous
16th problem of David Hilbert [1], unsolved for over a century, is concerned
with understanding periodic orbits in two-dimensional polynomial systems. In
three-dimensions, the problem is even more acute, since periodic orbits can
form together in a chaotic attractor. The existence or otherwise of chaos
within a system is a difficult problem with practical applications to
everything from engineering design to planetary sciences.

Some periodic orbits are born from Hopf bifurcations of stationary points, and
these are relatively easy to find. Others, however, appear in so-called global
bifurcations. If dynamically unstable, these are can be incredible difficult
to detect without a very close guess. One example of this is in the celebrated
Lorenz [2] equations:

$\displaystyle\begin{split}\frac{\mathrm{d}x_{1}}{\mathrm{d}t}&=\sigma(x_{2}-x_{1}),\\\
\frac{\mathrm{d}x_{2}}{\mathrm{d}t}&=x_{1}(\rho-x_{3})-x_{2},\\\
\frac{\mathrm{d}x_{3}}{\mathrm{d}t}&=x_{1}x_{2}-\beta x_{3}.\end{split}$ (1)

These were originally derived as a model for Rayleigh-Bénard convection, and
the parameter $\rho$ represents a rescaled Rayleigh number. Like in the fluid
dynamical problem, complexity increases as $\rho$ increases (for moderate
values), initially having a single globally stable stationary point and
eventually exhibiting chaos.

The Lorenz system is most often studied at parameter values $\beta=8/3$,
$\sigma=10$, and $\rho>0$. For $0<\rho<1$, there is a single globally
attracting equilibrium at the origin. At $\rho=1$ this undergoes a
supercritical pitchfork bifurcation to produce two initially stable
equilibria. These remain stable until a subcritical Hopf bifurcation at
$\rho=\sigma(\sigma+\beta+3)/(\sigma-\beta-1)\approx 24.7$ [3]. When
$\rho=28$, this system exhibits the famous ‘butterfly wing’ chaotic attractor,
which persists as a strange invariant set, along with its infinity of unstable
periodic orbits, as low as $\rho=13.927$ [3, 4]. At this value, the
heteroclinic connections between the equilibria undergo a global bifurcation,
and an ‘explosion’ of periodic orbits occurs.

No periodic orbits have been found for $0<\rho<13.926$, but results
guaranteeing their nonexistence are limited. Most of the research has been
concerned with guaranteeing that chaos does not exist, which is a stronger
condition, but has more practical relevance, and can be investigated through
Lyapunov exponents. Martini  _et al._ [5][6] use a method of Angeli  _et al._
[7] to show that positive Lyapunov exponents on the attractor, i.e. chaos,
cannot exist for $\rho\leq
2\sqrt{\beta-1}(\beta+\sigma)(\beta+1)/\sigma\beta\approx 4.5$. Leonov  _et
al._ [8] proved that if $\rho<(\beta+1)(\beta+\sigma)/\sigma\approx 4.6$, all
bounded trajectories converge to one of the equilibria. Above this, they
proved an upper bound on the dimension of the global attractor which is
strictly greater than 2, and thus does not prohibit chaos. In this paper, we
will show that periodic orbits and chaos cannot exist for $0\leq\rho\leq 12$,
a significantly higher maximum value than those mentioned above. In fact, we
will prove global stability. This term, coined in this context by Kuznetsov
_et al._ [9], means that all trajectories in the system converge to one of the
system’s equilibria.

The paper is laid out as follows: In section II, we present a general Lyapunov
function method, and specialise this to the Lorenz system. Section III
presents the sum-of-squares relaxation of this method which is computationally
tractable. In section IV, we derive and apply a new method for the rigorous
validation of these results using interval arithmetic. Section V summarises
the work and discusses potential extensions to the methods described. An
appendix presents an analogous method for maps.

##  II Mathematical results

Let $\frac{\mathrm{d}x}{\mathrm{d}t}=f(x)$ be a dynamical system for a
continuously differentiable function $f:\mathbb{R}^{n}\to\mathbb{R}^{n}$, such
that the corresponding flow
$\varphi:\mathbb{R}^{n}\times\mathbb{R}\to\mathbb{R}^{n}$ with
$\frac{\mathrm{d}}{\mathrm{d}t}\varphi(x,t)=f(x)$ exists and is continuously
differentiable for all $x\in\mathbb{R}_{>0}^{n}$ and $t\in\mathbb{R}$. Recall
the following definition [10]:

######  Definition 1.

The $\omega$-limit set of a point $x$ with respect to a flow $\varphi$ is the
subset of state space which is visited at arbitrary times,

$\omega(x,\varphi)=\\{y\in\mathbb{R}^{n}\mid\exists(t_{n})\text{ such that
}t_{n}\to\infty,\,\varphi(x,t_{n})\to y\text{ as }n\to\infty\\}.$ (2)

A limit point of the system is $y\in\mathbb{R}^{n}$ such that
$y\in\omega(x,\varphi)$ for some $x$. The limit set of the system is the
collection of all such points.

######  Lemma 1.

Let $g:\mathbb{R}^{n}\to\mathbb{R}_{\geq 0}$ be a continuous, nonnegative
function defined over state space. Suppose there exists a continuously
differentiable function $V:\mathbb{R}^{n}\to\mathbb{R}$ which satisfies

$f(x)\cdot\nabla V(x)\geq g(x)$ (3)

for all $x\in\mathbb{R}^{n}$. Then $g(y)=0$ for all limit points $y$.

###### Proof.

Suppose that $V:\mathrm{R}^{n}\to\mathbb{R}$ satisfies the condition of the
theorem.

Choose any $x_{0}\in\mathbb{R}^{n}$ and $y\in\omega(x_{0},\varphi)$. In
particular, let $(t_{n})$ be a sequence of times such that
$\lim_{n\to\infty}\varphi(x_{0},t_{n})=y$. Then

$\displaystyle V(\varphi(x_{0},t_{n}))-V(\varphi(x_{0},t_{n-1}))$
$\displaystyle=\int_{t_{n-1}}^{t_{n}}\frac{\mathrm{d}}{\mathrm{d}t}V(\varphi(x_{0},t))\mathrm{d}t$
$\displaystyle=\int_{t_{n-1}}^{t_{n}}f(\varphi(x_{0},t))\cdot\nabla
V(\varphi(x_{0},t))\mathrm{d}t$
$\displaystyle\geq\int_{t_{n-1}}^{t_{n}}g(\varphi(x_{0},t))\mathrm{d}t$
$\displaystyle\geq 0.$

But $V(\varphi(x_{0},t_{n}))-V(\varphi(x_{0},t_{n-1}))\to 0$ as $n\to\infty$
by the continuity of $V$, so

$\int_{t_{n-1}}^{t_{n}}g(\varphi(x_{0},t))\mathrm{d}t\to 0\text{ as
}n\to\infty.$

Since $g(x)\geq 0$ for all $x$, we must therefore have $g(\varphi(x_{0},t))\to
0$ as $t\to\infty$, and in particular,

$g(y)=\lim_{n\to\infty}g(\varphi(x_{0},t_{n}))=0.$

∎

An immediate application of this result is the following:

######  Theorem 1.

Suppose there exists a smooth function $V:\mathbb{R}^{n}\to\mathbb{R}$ which
satisfies

$f(x)\cdot\nabla V(x)\geq\|f(x)\|^{2}$ (4)

for all $x\in\mathbb{R}^{n}$. Then the limit set of the system $f$ is the set
of stationary points.

###### Proof.

For every limit point $x$, we know from lemma 1 that $\|f(x)\|^{2}=0$ and thus
$f(x)=0$, i.e. the point is a stationary point. Conversely, any stationary
point is its own limit point. ∎

Notice that this implies $f\cdot\nabla V$ is strictly positive away from
stationary points, so this condition is a special case of a strict Lyapunov
function as discussed by Bárta  _et al._ [11], and they prove that this means
the system is a gradient system when restricted to the subset of non-
stationary points. Such systems are known by various names [9]. We prefer the
term gradient-like systems [12], since an explicit gradient is not easily
found. Nevertheless, the behaviour of such a system is particularly simple.
Certainly no periodic orbits can exist, since any periodic point is contained
in its own $\omega$-limit set, and these points must necessarily be non-
stationary. Similarly, quasiperiodic motion is forbidden in a gradient-like
system, though this becomes more relevant in higher dimensions than 3. There
is no universally agreed definition of chaos, but the most widely cited is
that of Devaney [13], who takes as an axiom that periodic points are dense
within a chaotic set. Therefore, gradient-like systems are non-chaotic. Note
that, however, a gradient-like system does not imply that all limit points are
isolated stationary points.

###  II.1 Application to the Lorenz system

If there are no unbounded solutions to a gradient-like system, for example
because of the existence of globally absorbing compact set (sometimes called
dissipativity in the sense of Levinson), we can further deduce that all
trajectories are attracted to the stationary set. This property has been
called global stability [9], and it is what we shall attempt to prove for the
Lorenz system. First recall that [2]:

######  Theorem 2.

The system 1 is dissipative in the sense of Levinson.

In fact, various different simple, globally attracting sets can be found, the
compact intersection of which can be used to estimate the shape of the
attractor [14].

For the remainder of this work – though there is no reason to suppose these
methods cannot be easily applied to other parameter values – we shall
concentrate on the classic case of $\sigma=10$ and $\beta=8/3$. This poses a
slight problem for the numerical validation, since $8/3$ cannot be exactly
represented using IEEE 754 floating point numbers. Therefore, we rescale 1 so
that $x_{1}^{\prime}=x_{1}/M$, $x_{2}^{\prime}=x_{2}/M$,
$x_{3}^{\prime}=x_{3}/M$ and $t^{\prime}=t/3$, and also let $r=(\rho-1)/M$, so
that (immediately dropping the primes),

$\displaystyle\begin{split}\frac{\mathrm{d}x_{1}}{\mathrm{d}t}&=30x_{2}-30x_{1},\\\
\frac{\mathrm{d}x_{2}}{\mathrm{d}t}&=3Mx_{1}r-3Mx_{1}x_{3}-3x_{2}+3x_{1},\\\
\frac{\mathrm{d}x_{3}}{\mathrm{d}t}&=3Mx_{1}x_{2}-8x_{3},\end{split}$ (5)

where $M$ is an arbitrary integer factor chosen the rescale the numerics [15,
16], which will be helpful later. The following result is an alternative to
theorem 1 which will lead to a less computationally complex sum-of-squares
program.

######  Proposition 1.

Let $f:\mathbb{R}^{3}\to\mathbb{R}^{3}$ be defined by

$f(x_{1},x_{2},x_{3})=\left(30x_{2}-30x_{1},\,3Mx_{1}r-3Mx_{1}x_{3}-3x_{2}+3x_{1},\,3Mx_{1}x_{2}-8x_{3}\right).$
(6)

If there exists $V:\mathbb{R}^{3}\to\mathbb{R}^{3}$ such that

$f(x)\cdot\nabla V(x)\geq(x_{2}-x_{1})^{2}$ (7)

then the system 5 is globally stable.

###### Proof.

The system is equivalent to 1 and therefore dissipative by theorem 2.

Suppose such a $V$ exists. By lemma 1, we then have that $x_{2}-x_{1}=0$ for
every limit point $x$. Then at a limit point

$0=\frac{\mathrm{d}}{\mathrm{d}t}(x_{2}-x_{1})=f_{2}(x)-f_{1}(x)=f_{2}(x)-30(x_{2}-x_{1})=f_{2}(x)=3Mx_{1}r-3Mx_{1}x_{3}-3x_{2}+3x_{1}.$

Therefore $x_{1}(r-x_{3})=0$.

If $x_{1}=x_{2}=0$, we have a one-dimensional system
$\frac{\mathrm{d}x_{3}}{\mathrm{d}t}=-8x_{3}$. Taking the antiderivative as a
gradient, all one-dimensional systems are gradient systems and thus all
bounded trajectories are attracted to stationary points.

Otherwise, $x_{3}=r$, and
$0=\frac{\mathrm{d}x_{3}}{\mathrm{d}t}=x_{1}^{2}-8r$, so
$x_{1}=x_{2}=\pm\sqrt{8r}$, which are the two non-trivial equilibria of the
system when $r>0$. So all limit points are stationary points.

Furthermore, the system is dissipative in the sense of Levinson so it is
globally stable. ∎

As stated, proposition 1 could only prove global stability for a particular
value of the parameter $r$. We would like to numerically prove global
stability for a range of $r$, which suggests treating it as a fourth variable
in the system. This leads to the following result:

######  Proposition 2.

Let $f:\mathbb{R}^{3}\times\mathbb{R}\to\mathbb{R}^{3}$ be defined by

$f(x_{1},x_{2},x_{3};r)=\left(30x_{2}-30x_{1},\,3Mx_{1}r-3Mx_{1}x_{3}-3x_{2}+3x_{1},\,3Mx_{1}x_{2}-8x_{3}\right).$
(8)

If there exists $V:\mathbb{R}^{3}\times\mathbb{R}\to\mathbb{R}$ and
$s:\mathbb{R}^{3}\times\mathbb{R}\to\mathbb{R}_{\geq 0}$ such that

$f(x;r)\cdot\nabla
V(x;r)\geq(x_{2}-x_{1})^{2}+(Mr+1-\rho_{min})(\rho_{max}-Mr-1)s(x;r)$ (9)

then the system 1 with $\beta=8/3$ and $\sigma=10$ is globally stable for
$\rho_{min}\leq\rho\leq\rho_{max}$.

###### Proof.

System 1 is equivalent to system 5 with $r=(\rho-1)/M$. If
$\rho_{min}\leq\rho\leq\rho_{max}$, then

$(Mr+1-\rho_{min})(\rho_{max}-Mr-1)\geq 0,$

and so

$f(x;r)\cdot\nabla V(x;r)\geq(x_{2}-x_{1})^{2}$

since $s\geq 0$. Therefore, by proposition 1, system 5 is globally stable and
hence also system 1. ∎

Introducing the second auxiliary function $s$ is an example of the generalised
S-procedure, which is frequently used in sum-of-squares optimization [17, 18].
Notice that with $r$ as a variable, the system 8 is a polynomial system with
integer coefficients which can be exactly represented by a computer (assuming
$M$ is not unreasonably large).

##  III Sum-of-squares relaxation

Finding Lyapunov functions for stability problems can be difficult, but one
systematic approach is through sum-of-squares optimization [19]. If the system
and the auxiliary functions are both restricted to be polynomials of
relatively low degree, the problem often becomes computationally tractable.
Verifying that a polynomial is non-negative is NP-hard [20], but verifying
that it is a sum of squares of polynomials, a stronger but sufficient
condition, is possible through semi-definite programming [21].

Let $P_{n}$ be the space of polynomials in $n$ real variables (of any degree).
For a polynomial $f\in P_{n}$, we say that $f$ is a sum-of-squares, written
$f\in\Sigma_{n}$, if there exist polynomials $p_{1},...,p_{m}\in P_{n}$ such
that

$f=\sum_{k=1}^{m}p_{k}^{2}.$

Clearly $f\in\Sigma_{n}$ is a non-negative polynomial, though in general the
converse is not true. Let $\mathbb{S}_{m}$ be the space of $m\times m$
symmetric matrices which are positive semi-definite. An important result [22]
is that $f\in\Sigma_{n}$ if and only if there exists a (non-unique)
$Q\in\mathbb{S}_{m}$ and a choice of ‘basis’ functions $b_{1},\dots,b_{m}\in
P_{n}$ such that $f=\mathbf{b}^{T}Q\mathbf{b}$. The matrix $Q$ is called the
Gram matrix.

Figure 1: Dual objective value reported by the SDP solver for a naïve
implementation of proposition 1 (with $M=16$). A value close to zero indicates
that the SDP is likely to be feasible, i.e. that the Lorenz system is indeed
gradient-like at that value of $\rho$.

A straightforward relaxation of proposition 1 gives a sum-of-squares program,
which itself is equivalent to a semi-definite program (SDP) which can be
solved numerically with a number of existing software packages. We employ the
Julia package SumOfSquares.jl [23, 24], on top of the MOSEK SDP optimizer
[25]. Figure 1 shows the dual objective value, effectively a residual,
returned by MOSEK for this SDP, with different choices of $\rho$ and
polynomial degree. A dual objective value close to zero indicates that the
problem is likely to be feasible. For $d=4$ and $d=6$, a significant jump in
the value is observed around $\rho=4$ and $\rho=5.5$ respectively, indicating
where a Lyapunov function of that degree is no longer sufficient. For $d=8$, a
clear jump is observed around $\rho=12$, and increasing the degree further
does not improve this significantly. These results suggest that the Lorenz
system is gradient-like up to at least $\rho=12$. However, these numerical
results are only indicative and are not proofs, and furthermore they are valid
only at a a given fixed value of $\rho$ and are not useful for the general
results we want to prove. Therefore we relax proposition 2 to the following:

######  Proposition 3.

Let $f_{1},f_{2},f_{3}\in P_{4}$ be defined by

$f(x_{1},x_{2},x_{3};r)=\left(30x_{2}-30x_{1},\,3Mx_{1}r-3Mx_{1}x_{3}-3x_{2}+3x_{1},\,3Mx_{1}x_{2}-8x_{3}\right).$
(10)

If there exist polynomials $V,b_{1},\dots,b_{m_{2}}\in P_{4}$ and
$s,a_{1},\dots,a_{m_{1}}\in P_{3}$ and matrices $Q\in\mathbb{S}_{m_{1}}$ and
$R\in\mathbb{S}_{m_{2}}$ such that

$\displaystyle s$ $\displaystyle=\mathbf{a}^{T}Q\mathbf{a},$ (11)
$\displaystyle
f_{1}\partial_{1}V+f_{2}\partial_{2}V+f_{3}\partial_{3}V-(x_{2}-x_{1})^{2}-M^{2}(r-r_{min})(r_{max}-r)s$
$\displaystyle=\mathbf{b}^{T}R\mathbf{b},$ (12)

then the system 1 with $\beta=8/3$ and $\sigma=10$ is globally stable for
$\rho_{min}\leq\rho\leq\rho_{max}$, where $\rho_{min}=Mr_{min}+1$ and
$\rho_{max}=Mr_{max}+1$.

For fixed bases $\mathbf{a}$ and $\mathbf{b}$, this is a semi-definite
program, which can be solved numerically for the coefficients of $Q$ and $R$
and the coefficients of the polynomials $V$ and $s$ with respect to the
monomial basis. Nevertheless, such optimisations are slow, imprecise and
difficult to validate if applied naively, and can be improved by a careful
choice of bases.

The first observation is that 5 is invariant under the symmetry transformation
$(x_{1},x_{2},x_{3},r)\mapsto(-x_{1},-x_{2},x_{3},r)$, and so if this is true
for $V$ and $s$, the left-hand sides of 11 is also invariant. In order for the
right-hand sides to be invariant under this symmetry, the Gram matrices are
necessarily block diagonal, and we can subdivide the bases and matrices such
that

$\displaystyle\mathbf{a}^{T}Q\mathbf{a}=\mathbf{a_{e}}^{T}Q_{e}\mathbf{a_{e}}+\mathbf{a_{o}}^{T}Q_{o}\mathbf{a_{o}},$
$\displaystyle\mathbf{b}^{T}R\mathbf{b}=\mathbf{b_{e}}^{T}R_{e}\mathbf{b_{e}}+\mathbf{b_{o}}^{T}R_{o}\mathbf{b_{o}},$

where $\mathbf{a_{e}}$ and $\mathbf{b_{e}}$ are invariant under the symmetry
transformation and $\mathbf{a_{o}}$ and $\mathbf{b_{o}}$ change sign. Since
$\mathbf{a_{e}}$ and $\mathbf{a_{o}}$ are roughly half the size of
$\mathbf{a}$, we have reduced the number of unknown matrix entries by a factor
of four, which results in a very significant speedup of both SDP and
validation.

Secondly, notice that at a stationary point of the system,
$f_{1}\partial_{1}V+f_{2}\partial_{2}V+f_{3}\partial_{3}V-(x_{2}-x_{1})^{2}=0$.
Therefore, for the equalities and definiteness to hold, we must have that
$s=0$ and also $\mathbf{a}^{T}Q\mathbf{a}=\mathbf{b}^{T}R\mathbf{b}=0$ at
stationary points. This suggests choosing
$\mathbf{a_{e}},\mathbf{b_{e}},\mathbf{a_{o}},\mathbf{b_{o}}$ which themselves
vanish at the stationary points, as was done (at much lower degree) in the
analytic results of Goluskin [15].

One further optimisation would be to choose a specific form of the polynomials
$V$ and $s$, following Goluskin [15]. This was found to make little difference
to performance, so for clarity of code we simply used a set monomials
invariant under the symmetry transformations. The Lyapunov function $V$ was
taken to be a polynomial of degree at most $d_{x}$ in $x_{1}$, $x_{2}$ and
$x_{3}$ combined, and of maximum degree $d_{r}$ in $r$ (so that terms like
$x_{1}^{d_{x}}r$ are allowed, but not $x_{1}^{d_{x}}x_{2}$). The function $s$
was taken to be a polynomial of maximum degree $d_{x}$ in $x$ and maximum
degree $d_{r}-1$ in $r$.

For the validation procedure described in section IV to be successful, it will
be necessary that the Gram matrices are strictly positive definite, rather
than simply semi-definite. In other words, we need non-singular Gram matrices.
We can attempt to enforce this by a parsimonious choice of basis, so that no
unnecessary zero eigenvectors are included. Successful bases were found
through trial-and-error, using the principles above, and these are listed in
table 1 for $d_{x}=8$. Notice that all of the polynomials listed vanish at
each of the system’s three equilibria. Partial results, with lower
$\rho_{max}$ and lower $d_{x}$ were also successful; the bases for these are
included in the source code.

$\mathbf{a}_{e}$ | $(x_{2}-x_{1})(x_{1},x_{2},x_{1}x_{3},x_{2}x_{3},x_{1}x_{3}^{2},x_{2}x_{3}^{2}),(3Mx_{1}x_{2}-8x_{3})(1,x_{3},x_{1}^{2},x_{1}x_{2},x_{2}^{2},x_{3}^{2})$
---|---
$\mathbf{a}_{o}$ | $(x_{2}-x_{1})(1,x_{3},x_{1}^{2},x_{1}x_{2},x_{2}^{2},x_{3}^{2},x_{1}^{2}x_{3},x_{1}x_{2}x_{3},x_{2}^{2}x_{3}),(3Mx_{1}x_{2}-8x_{3})(x_{1},x_{2},x_{1}x_{3},x_{2}x_{3})$
|
$(x_{2}-x_{1})(x_{1},x_{2},x_{1}x_{3},x_{2}x_{3},x_{1}x_{3}^{2},x_{2}x_{3}^{2},rx_{1},rx_{2},rx_{1}x_{3},x_{2}x_{3},rx_{1}x_{3}^{2},rx_{2}x_{3}^{2}),$
$\mathbf{b}_{e}$ |  $(r-x_{3})(x_{3},x_{1}^{2},x_{1}x_{2},x_{2}^{2},x_{3}^{2},x_{1}^{2}x_{3},x_{1}x_{2}x_{3},x_{2}^{2}x_{3})$,
|
$(3Mrx_{2}^{2}-8x_{3}^{2})(1,x_{3},x_{1}^{2},x_{1}x_{2},x_{2}^{2},x_{3}^{2})$
|
$(x_{2}-x_{1})(1,x_{3},x_{1}^{2},x_{1}x_{2},x_{2}^{2},x_{3}^{2},x_{1}^{2}x_{3},x_{1}x_{2}x_{3},x_{2}^{2}x_{3},r,rx_{3},rx_{1}^{2},rx_{1}x_{2},rx_{2}^{2},rx_{3}^{2},rx_{1}^{2}x_{3},rx_{1}x_{2}x_{3})$,
$\mathbf{b}_{o}$ |  $(r-x_{3})(x_{1},x_{2},x_{1}x_{3},x_{2}x_{3},x_{1}x_{3}^{2},x_{2}x_{3}^{2})$,
| $(3Mrx_{2}^{2}-8x_{3}^{2})(x_{1},x_{2},x_{1}x_{3},x_{2}x_{3})$
Table 1: Bases of polynomials used for Gram matrices, for Lyapunov function of
degree $(d_{x},d_{r})=(8,1)$. Notice that $\mathbf{a}_{e}$ and
$\mathbf{a}_{o}$ have maximum degree 4 in $x$ and do not depend on $r$, and
$\mathbf{b}_{e}$ and $\mathbf{b}_{o}$ have maximum degree 4 in $x$ and 1 in
$r$.

##  IV Validation with interval arithmetic

The rigorous validation of numerical results of dynamical systems theory
through the use of interval arithmetic have been some of the major successes
of the field in the last few decades [26, 27, 28]. However, once it has been
demonstrated to be possible, it is not usually considered necessary to
validate every single result. For example, if one has converged a periodic
orbit to machine precision for a given system, and the location and period of
this periodic orbit appears to be robust to variations in the timestep or
discretisation of the numerical solution, then there is every reason to expect
that this is a genuine solution to the underlying differential equations
without the laborious process of validation through interval arithmetic.
Unfortunately this is not the case for sum-of-squares optimisation. The
underlying semidefinite programs are known to be imprecise and very sensitive
to numerical artifacts. Furthermore, sum-of-squares programs are observed to
generate particularly poorly-conditioned SDPs, likely due to the fact that
very high degree exponents are applied to finite precision floating number.
Solvers such as MOSEK can report success and return solutions with small
residuals, even when the problem is in fact infeasible, or report failure but
return plausible-seeming results. Therefore it is imperative that we can
validate our numerics.

Several different procedures for validating sum-of-squares relaxations of this
type are found in the literature. The simplest uses the numerical result to
write down an analytic Lyapunov function whose validity is confirmed
explicitly [15]. This is practical only in relatively simple cases, although
specially designed computer-algebra packages make this easier [29]. A
different method is to project the problem onto the space of sums-of-squares
of polynomials with rational coefficients [30], which can be computed exactly
and for which a package is also available [31]. A third option is to use
interval arithmetic to validate an approximate numerical result. For
optimisation problems, one can find rigorous, but not sharp, upper- and lower-
bounds on a quantity of interest by perturbing the result of an SDP solver
[32]. This has been applied to find rigorous bounds for dynamical systems
[15]. However, the method discussed in the present work is equivalent to
finding upper- and lower-bounds of exactly zero, so this cannot be employed
directly. We therefore derive a different, simpler method to validate an
approximate result for a sum-of-squares feasibility problem with interval
arithmetic, which formulates the problem as a linear system.

###  IV.1 A semi-definite program as a linear system

Consider the problem of finding $s\in P_{n}$ and $Q\in\mathbb{S}_{m}$ such
that

$s=\mathbf{a}^{T}Q\mathbf{a},$ (13)

for some given set of polynomials $a_{1},...,a_{m}\in P_{n}$. (With no other
conditions this example is of course trivial.) Let
$\mathbf{q}\in\mathbb{R}^{m(m+1)/2}$ be the vector of unique elements of the
symmetric Gram matrix Q, and let $\mathbf{s}\in\mathbb{R}^{N}$ be the
coefficients of $s$ with respect to the monomial basis. Equation 13 is linear
in these coefficients, and so expanding it in monomials and comparing
coefficients, we can formulate this as a system

$A\left(\begin{matrix}\mathbf{s}\\\ \mathbf{q}\end{matrix}\right)=\mathbf{c}$
(14)

for some matrix of coefficients $A\in\mathbb{R}^{N\times(N+m(m+1)/2)}$ and
vector of constant terms $\mathbf{c}\in\mathbb{R}^{N}$. (In this example
$\mathbf{c}=0$.) Suppose that the rank of $A$ is $M$, and that we have some
approximate solution $\mathbf{s}_{0}$, $\mathbf{q}_{0}$ such that
$A\left(\begin{matrix}\mathbf{s_{0}}\cr\mathbf{q_{0}}\end{matrix}\right)-\mathbf{c}=\mathbf{\epsilon}$
is small. Then we can project onto the solution space of the system using

$\left(\begin{matrix}\mathbf{s}\cr\mathbf{q}\end{matrix}\right)=A^{+}c+(I-A^{+}A)\left(\begin{matrix}\mathbf{s_{0}}\cr\mathbf{q_{0}}\end{matrix}\right)=\left(\begin{matrix}\mathbf{s_{0}}\cr\mathbf{q_{0}}\end{matrix}\right)-A^{+}\mathbf{\epsilon}$
(15)

where $A^{+}$ is the Moore-Penrose right inverse, i.e. $AA^{+}=I$.
Reassembling this into a polynomial $s$ with coefficients $\mathbf{s}$ and a
Gram matrix $Q$ with coefficients $\mathbf{q}$ , we have a solution to eq. 13.
There is no guarantee that the resulting $Q$ is positive semi-definite, but
since the matrix $Q_{0}$ with coefficients $\mathbf{q_{0}}$ is positive semi-
definite, there is reason to hope $Q$ will be too if the initial error
$\mathbf{\epsilon}$ is sufficiently small. Our strategy will be to find a
vector of intervals which is guaranteed to contain a solution to the linear
system, and prove that all $\mathbf{q}$ in these intervals lead to a positive
semi-definite Gram matrix, as summarised in the following:

######  Proposition 4.

Let $\mathbf{v}$ be the coefficients of the unknown polynomial functions
$V(x,r)$ with respect to the basis of monomials up to degree $d_{x}$ in $x$
and $d_{r}$ in $r$. Let $\mathbf{q}_{e}$, $\mathbf{r}_{e}$, $\mathbf{q}_{o}$
and $\mathbf{r}_{o}$ be the unique elements of the unknown real symmetric
matrices $Q_{e}$, $R_{e}$, $Q_{o}$ and $R_{o}$ respectively, and let $A$ and
$\mathbf{c}$ be the (known) matrix and vector to represent the polynomial
equation

$f_{1}\partial_{1}V+f_{2}\partial_{2}V+f_{3}\partial_{3}V-(x_{2}-x_{1})^{2}=M^{2}(r-r_{min})(r_{max}-r)\left(\mathbf{a_{e}}^{T}Q_{e}\mathbf{a_{e}}+\mathbf{a_{o}}^{T}Q_{o}\mathbf{a_{o}}\right)+\left(\mathbf{b_{e}}^{T}R_{e}\mathbf{b_{e}}+\mathbf{b_{o}}^{T}R_{o}\mathbf{b_{o}}\right)$
(16)

as the linear system

$A\left(\begin{matrix}\mathbf{v}\cr\mathbf{q}_{e}\cr\mathbf{r}_{e}\cr\mathbf{q}_{o}\cr\mathbf{r}_{o}\end{matrix}\right)=\mathbf{c}.$
(17)

If $A$ has full row-rank, and, for some vector $\mathbf{y}_{0}$, the solution

$\left(\begin{matrix}\mathbf{v}\cr\mathbf{q}_{e}\cr\mathbf{r}_{e}\cr\mathbf{q}_{o}\cr\mathbf{r}_{o}\end{matrix}\right)=A^{+}\mathbf{c}+(I-A^{+}A)\mathbf{y_{0}}$
(18)

leads to symmetric matrices $Q_{e}$, $R_{e}$, $Q_{o}$ and $R_{o}$ which are
positive semi-definite, then the system 5 is globally stable for $r_{min}\leq
r\leq r_{max}$.

Note that we have entirely eliminated the additional auxiliary function $s$
from this problem, it is used only to find an initial guess $\mathbf{y_{0}}$
through approximate sum-of-squares optimisation.

One key point here is that the matrix $A$ must be of full row-rank. It is
actually the case that by naïvely assembling $A$ by comparing the coefficients
of each of the possible monomials, one arrives at a matrix which has
degenerate rows. At a stationary point, where $f(x_{1},x_{2},x_{3},r)=0$,
every term in the polynomial expression 16 vanishes. This is by design, by
construction of $\mathbf{a}_{e}$, $\mathbf{b}_{e}$, $\mathbf{a}_{o}$ and
$\mathbf{b}_{o}$ as discussed in section III. At a stationary point, either
$x_{1}=x_{2}=x_{3}=0$ (which we can ignore by not considering monomials for
which none of these appear), or $x_{2}=x_{1}$ and $r=x_{3}=3Mx_{1}^{2}/8$.
Substituting these expressions into eq. 16, all terms vanish. However,
subsituting these expressions in to the set of monomials which are used to
generate eq. 17, we arrive at a set of monomials in $x_{1}$ whose coefficients
must vanish. For example, the monomials $x_{1}^{2}$, $x_{1}x_{2}$,
$x_{2}^{2}$, $x_{3}$ and $r$ all reduce to a monomial $x_{1}^{2}$ at a
stationary point, and therefore the sum (up to a constant factor of $3M/8$) of
the rows in $A$ which correspond to these monomials must be zero. This is
straightforward to remedy: simply eliminate from eq. 17 one row for each power
of $x_{1}$ in this reduced set of monomials – since each row must be a linear
combination of the others for that power of $x_{1}$.

###  IV.2 Interval arithmetic

The basic idea of numerical interval arithmetic is, rather than using
imprecise floating-point numbers directly for computations, to use intervals
whose end-points are floating-point numbers but which certainly contain the
exact result. For example, adding two intervals results in an interval which
is typically slightly larger than the interval given by the approximate sum of
the endpoints of the two intervals. A comprehensive background on interval
arithmetic is outside the scope of this work. The textbooks by Jaulin  _et
al._ [33] and Tucker [34] give a good introduction to the field.

We have carefully constructed the problem to be validated in proposition 4 so
that the matrix $A$ and right-hand side $\mathbf{c}$ can be exactly expressed
either as integers or as IEEE double-precision floating point numbers.
However, constructing a right-inverse for $A$ with floating-point arithmetic
will almost certainly lead to an imprecise result. For that reason, this step
is performed with interval arithmetic, leading to an interval matrix which is
guaranteed to contain $A^{+}$.

Let $\mathbb{IR}$ be the space of closed but not necessarily bounded intervals
on the real line, specified by their endpoints. Addition and subtraction are
defined in the obvious ways, by addition and subtraction of the endpoints, and
multiplication is defined to properly account for signs when intervals are
negative or contain zero [34]. These operations are implemented numerically
using the Julia package IntervalArithmetic.jl [35], so that the result of any
calculation is guaranteed to contain the exact result. We denote by square
brackets $[\cdot]$ an interval vector in $\mathbb{IR}^{n}$ (or interval matrix
in $\mathbb{IR}^{n\times m}$) which contains the true value of a quantity in
$\mathbb{R}^{n}$ (resp. $\mathbb{R}^{n\times m}$).

The right-inverse for a real matrix is usually expressed as

$A^{+}=A^{T}(AA^{T})^{-1},$

where $A^{T}$ is the usual transpose. In fact, it is easier to calculate as

$A^{+}=\left(\left(AA^{T}\right)^{-1}A\right)^{T}.$

First we calculate an interval matrix $[AA^{T}]$ strictly enclosing the
product $AA^{T}$ (both of which are known exactly), and we then use interval
Gaussian elimination [33] with $[AA^{T}]$ and $A$ to find an interval matrix
strictly containing $\left(A^{+}\right)^{T}$, which can be exactly transposed
to give $[A^{+}]$. One important detail, as mentioned in the previous section,
is that for this to be possible, the rectangular matrix $A$ must have full
row-rank, otherwise the Gaussian elimination will give an error or result in
infinite intervals and the validation will ultimately fail.

The final step of our validation will be to check that the Gram interval
matrices are positive definite. There are a number of different equivalent
conditions for a real symmetric matrix to be positive definite. The easiest to
implement with interval arithmetic is to check that the determinants of all
the principal minors, i.e. the upper-left $1\times 1$ submatrix, the upper-
left $2\times 2$ submatrix etc., are positive. Since determinants are
calculated with addition, subtraction and multiplication, no new code is
required to compute these for interval matrices. Notice that for a semi-
definite matrix which is not strictly definite, the determinant is zero. When
the determinant is calculated using interval arithmetic, the result will be an
interval containing zero which is not strictly non-negative, and thus our
validation will be unsuccessful. This is the reason that we must restrict to
strictly positive definite matrices.

The sum-of-squares optimisation was performed using 64-bit ‘double’ precision
floating point numbers, to provide an initial guess. However, for validating
these results, using 64-bit endpoints to represent the intervals was
insufficient and resulted in determinants which were slightly negative.
128-bit floating point numbers were sufficient for almost all our results, but
for higher values of $\rho$ it was necessary to use 256 bits. This is
straightforward, if rather slow, with Julia.

###  IV.3 The full validation for the Lorenz system

Combining all the parts of the previous sections, the algorithm proceeds as
follows:

  1. 1.

The sum-of-squares program described in proposition 3 is formulated and solved
approximately, with careful choices of basis for the right-hand side as listed
in table 1, leading to a polynomial Lyapunov function $V$ with coefficients
$\mathbf{v}$ and symmetric Gram matrices $Q_{e}$, $R_{e}$, $Q_{o}$ and $R_{o}$
with coefficients $\mathbf{q}_{e}$, $\mathbf{r}_{e}$, $\mathbf{q}_{o}$ and
$\mathbf{r}_{o}$ respectively.

  2. 2.

A large rectangular linear system $A\mathbf{y}=\mathbf{c}$ is constructed as
in proposition 4, where $A$ and $\mathbf{c}$ are known exactly, and with
approximate solution
$\mathbf{y}_{0}=\left(\mathbf{v}^{T}\,\mathbf{q}_{e}^{T}\,\mathbf{r}_{e}^{T}\,\mathbf{q}_{o}^{T}\,\mathbf{r}_{o}^{T}\right)^{T}$
given by the result of the SOS program.

  3. 3.

Degenerate rows are removed from this linear system to produce an equivalent
system with full row rank $\hat{A}\mathbf{y}=\hat{\mathbf{c}}$.

  4. 4.

An interval matrix $[\hat{A}^{+}]$ enclosing the Moore-Penrose right-inverse
$\hat{A}^{+}$ is calculated. This step also verifies that $\hat{A}$ has full
row-rank.

  5. 5.

An interval vector $[\mathbf{y}]$ enclosing the true solution of
$\hat{A}\mathbf{y}=\hat{\mathbf{c}}$ is computed as
$[\mathbf{y}]=[\hat{A}^{+}]\hat{\mathbf{c}}+(I-[\hat{A}^{+}]A)\mathbf{y}_{0}$.

  6. 6.

The elements of $[\mathbf{y}]$ corresponding to the Gram matrices of the sum-
of-squares program are used to assemble symmetric interval matrices $[Q_{e}]$,
$[R_{e}]$, $[Q_{o}]$ and $[R_{o}]$.

  7. 7.

The determinant intervals of the principal minors of these Gram interval
matrices are calculated. If they are all strictly positive, the matrices are
positive-definite and the validation is successful.

The bottleneck during our computations were the calculations using large
interval matrices with very high precision - steps 4 and 5. The fully
validated results are listed in table 2. At least with the 8th degree
polynomials used, it was necessary to subdivide the ranges of $\rho$ into
progressively smaller intervals as the values increased. Combining these
results, we have proved:

######  Theorem 3.

The Lorenz system 1 with $\sigma=10$ and $\beta=8/3$ is globally stable for
all $\rho\in[0,12]$.

$\rho_{min}$ | $\rho_{max}$ | $d_{x}$ | $d_{r}$
---|---|---|---
0 | 2 | 4 | 1
0 | 4 | 6 | 1
0 | 6 | 8 | 1
6 | 10 | 8 | 1
10 | 11 | 8 | 1
11 | 11.5 | 8 | 1
11.5 | 11.75 | 8 | 1
11.75 | 11.875 | 8 | 1
11.875 | 12 | 8 | 1
Table 2: Successfully validated globally stable regions for
$\rho\in[\rho_{min},\rho_{max}]$, with the corresponding degrees used for the
Lyapunov function $V$.

##  V Conclusion

This work presents three new results. Firstly, a new Lyapunov function method
has been described which can be used to prove a system is gradient-like, and
in particular the absence of chaos, in continuous time dynamical systems.
Secondly, this method has been applied to the Lorenz system up to $\rho=12$,
using the inexact computational method of sum-of-squares optimisation, greatly
extending the known parameter values for which the system is non-chaotic.
Thirdly, the result has been rigorously validated using a new procedure based
on interval arithmetic.

It is interesting to note that both the fixed-$\rho$ numerical results and the
fully validated results failed to demonstrate gradient-like behaviour above
$\rho=12$, even with significantly higher degree polynomials. There are
several possible reasons for this. One is that this method would be possible
in that range, with more carefully controlled numerics and higher-precision
floating point numbers. The second is that the behaviour there is gradient-
like, but that it is not possible to use our method to prove it, being only a
sufficient condition. The third possibility is that there are some as-yet
undiscovered periodic orbits or other non-trivial behaviour which exists above
$\rho=12$ but below the known periodic orbits arising in the homoclinic
explosion at $\rho\approx 13.93$.

The validation procedure described in section IV is essentially proving a
sharp bound of 0 for the method of Fantuzzi  _et al._ [18]. This could
straightforwardly be extended to prove sharp bounds which are rational
numbers, if care is taken to formulate the problem so that the bound is
exactly specified as a floating-point number. However, it must be stressed
that this method took a great deal of trial-and-error to implement, in order
to choose the correct basis of polynomials, and so further development in this
area is to be encouraged.

One potential extension of this work is to partial differential equations.
Previous authors [36, 37] have applied sum-of-squares optimisation methods to
dissipative PDEs using a suitably chosen function space to construct auxiliary
functions. Initial experimentation suggests that the present method is equally
applicable in such high-dimensional systems, though transforming the numerical
results into rigorous ones becomes significantly more complicated.

## Acknowledgements

The author is indebted to David Goluskin who provided encouragement and advice
on this work. Thanks are also due to Spyridon Pougkakiotis for useful
comments. Olivier Hénot gave helpful suggestions for the Julia implementation.
The preliminary phase of this project was supported by the European Research
Council (ERC) under the European Union’s Horizon 2020 research and innovation
programme (grant no. 865677), while the author was a member of the Emergent
Complexity in Physical Systems laboratory at the École Polytechnique Fédérale
de Lausanne.

The full source code to produce the figures in this paper and to perform the
validation required for theorem 3 is available at
https://github.com/jeremypparker/lorenzstability.

##  Appendix A An auxiliary function method for maps

Consider a discrete-time dynamical system, that is a continuous map
$f:\mathbb{R}^{n}\to\mathbb{R}^{n}$, and define a trajectory $(x_{n})$ by

$x_{n+1}=f(x_{n}).$ (19)

######  Definition 2.

The $\omega$-limit set of a point $x$ with respect to a map $f$ is the subset
of state space which is visited at arbitrarily high iterations,

$\omega(x,f)=\\{y\in\mathbb{R}^{n}\mid\exists(n_{k})\in\mathbb{N}\text{ with
}n_{k+1}>n_{k}\text{ such that }f^{n_{k}}(x)\to y\text{ as }k\to\infty\\}.$
(20)

A limit point of the system is $y\in\mathbb{R}^{n}$ such that
$y\in\omega(x,f)$ for some $x$.

######  Lemma 2.

Let $g:\mathbb{R}^{n}\to\mathbb{R}_{\geq 0}$ be some non-negative function
defined continuously over state space. Suppose there exists a continuous
function $V:\mathbb{R}^{n}\to\mathbb{R}$ such that, for all $x$,

$V(f(x))-V(x)\geq g(x).$ (21)

Then $g(y)=0$ for all limit points $y$.

###### Proof.

Let $y$ be a limit point and $(x_{n})$ be a trajectory with a sequence of
iterations $(n_{k})$ such that $x_{n_{k}}\to y$. Since $g(x)\geq 0$, we have
that $V(f(x))\geq V(x)$. Then

$\displaystyle V(x_{n_{k+1}})-V(x_{n_{k}})$
$\displaystyle=V\left(f^{n_{k+1}-n_{k}}(x_{n_{k}})\right)-V(x_{n_{k}})$
$\displaystyle\geq V\left(f(x_{n_{k}})\right)-V(x_{n_{k}})$ $\displaystyle\geq
g(x_{n_{k}})\geq 0.$

As $k\to\infty$, $V(x_{n_{k+1}})-V(x_{n_{k}})\to 0$, and $g(x_{n_{k}})\to
g(y)$. So $0\geq g(y)\geq 0$, and thus $g(y)=0$. ∎

Now taking $g(x)=\left\|f(x)-x\right\|^{2}$, if such a $V$ exists we could use
this result to show that all limit points are fixed points, and no other
periodic points can exist. Furthermore, we could, for example, show that
periodic points of period 1 (i.e. equilibria) or period $2$ can exist, but no
higher, with $g(x)=\left\|f^{2}(x)-x\right\|^{2}$. Any such result would show
that the system is non-chaotic. In general let
$g_{k}(x)=\left\|f^{k}(x)-x\right\|^{2}$.

A sum-of-squares relaxation is again possible for polynomial $f$. We would
need to find a $V$ such that

$V(f(x))-V(x)-g(x)\in\Sigma_{n}.$ (22)

Due to the composition of $V$ with $f$, the polynomial degrees rapidly become
very large as the degree of $V$ is increased, which means that approximate
solution and rigorous validation is both computationally expensive and
requires high precision. Nevertheless, we can demonstrate the method in the
simple case of the Hénon map, defined by

$f(x_{1},x_{2})=\left(1-ax_{1}^{2}+x_{2},bx_{1}\right).$ (23)

We fix $b=0.3$, the standard value, and vary $a$. As $a$ increases from $0$ to
the classical value of $1.4$, the initial fixed point undergoes a cascade of
period-doubling bifurcations until the onset of chaos around $a\approx 1.1$.

Figure 2 shows the results of an implementation of this method with different
possible degrees $d$ of the polynomial auxiliary function $V$ and different
numbers $k$ of iterations of the map. Even for $k=3$ the problem becomes very
difficult and slow, but for $k=1$ we see a clear increase in the residual
reported by the semi-definite solver (MOSEK) at the location of the first
bifurcation, and for $k=2$ we see a clear increase at the second bifurcation
when $d$ is sufficiently large. The SDP solver reports nonzero residuals, so
an exact solution has not been found and these results can only be taken as
numerical evidence that it would be possible to find a suitable auxiliary
function. Nevertheless, we have every reason to believe that a similar
approach to that in section IV would be possible here, to rigorously validate
these results.

Figure 2: Bifurcation diagram (black, left axis) and SOS results (colours) for
the Hénon map at $b=0.3$. For each choice of polynomial degree $d$ and map
iterate $k$, the right vertical axis shows the dual objective value reported
by MOSEK. A value close to zero indicates that the problem is likely to be
feasible.

## References

  * Hilbert [1902] D. Hilbert, Bulletin of the American Mathematical Society 8, 437 (1902).
  * Lorenz [1963] E. N. Lorenz, Deterministic nonperiodic flow, Journal of Atmospheric Sciences 20, 130 (1963).
  * Sparrow [2012] C. Sparrow,  _The Lorenz equations: bifurcations, chaos, and strange attractors_ , Vol. 41 (Springer Science & Business Media, 2012).
  * Chen [1996] X. Chen, Lorenz equations part i: existence and nonexistence of homoclinic orbits, SIAM Journal on Mathematical Analysis 27, 1057 (1996).
  * Martini  _et al._ [2022] D. Martini, D. Angeli, G. Innocenti, and A. Tesi, Ruling out positive Lyapunov exponents by using the Jacobian’s second additive compound matrix, IEEE Control Systems Letters 6, 2924 (2022).
  * Martini  _et al._ [2023] D. Martini, D. Angeli, G. Innocenti, and A. Tesi, Bounding Lyapunov exponents through second additive compound matrices: Case studies and application to systems with first integral, International Journal of Bifurcation and Chaos 33, 2350114 (2023).
  * Angeli  _et al._ [2021] D. Angeli, M. A. Al-Radhawi, and E. D. Sontag, A robust Lyapunov criterion for nonoscillatory behaviors in biological interaction networks, IEEE Transactions on Automatic Control 67, 3305 (2021).
  * Leonov  _et al._ [2016] G. A. Leonov, N. V. Kuznetsov, N. Korzhemanova, and D. Kusakin, Lyapunov dimension formula for the global attractor of the Lorenz system, Communications in Nonlinear Science and Numerical Simulation 41, 84 (2016).
  * Kuznetsov  _et al._ [2020] N. V. Kuznetsov, M. Y. Lobachev, M. V. Yuldashev, R. V. Yuldashev, E. V. Kudryashova, O. A. Kuznetsova, E. Rosenwasser, and S. Abramovich, The birth of the global stability theory and the theory of hidden oscillations, in  _2020 European Control Conference (ECC)_ (IEEE, 2020) pp. 769–774.
  * Glendinning [1994] P. Glendinning,  _Stability, instability and chaos: an introduction to the theory of nonlinear differential equations_ (Cambridge University Press, 1994).
  * Bárta  _et al._ [2012] T. Bárta, R. Chill, and E. Fašangová, Every ordinary differential equation with a strict Lyapunov function is a gradient system, Monatshefte für Mathematik 166, 57 (2012).
  * Hale [2004] J. K. Hale, Stability and gradient dynamical systems., Revista matemática complutense 17, 7 (2004).
  * Devaney [2018] R. Devaney,  _An introduction to chaotic dynamical systems_ (CRC press, 2018).
  * Doering and Gibbon [1995] C. R. Doering and J. Gibbon, On the shape and dimension of the Lorenz attractor, Dynamics and Stability of Systems 10, 255 (1995).
  * Goluskin [2018] D. Goluskin, Bounding averages rigorously using semidefinite programming: mean moments of the Lorenz system, Journal of Nonlinear Science 28, 621 (2018).
  * Parker  _et al._ [2021] J. P. Parker, D. Goluskin, and G. M. Vasil, A study of the double pendulum using polynomial optimization, Chaos: An Interdisciplinary Journal of Nonlinear Science 31 (2021).
  * Tan [2006] W. Tan,  _Nonlinear Control Analysis and Synthesis using Sum-of-Squares Programming_ , Phd thesis, University of California, Berkley (2006).
  * Fantuzzi  _et al._ [2016] G. Fantuzzi, D. Goluskin, D. Huang, and S. I. Chernyshenko, Bounds for deterministic and stochastic dynamical systems using sum-of-squares optimization, SIAM Journal on Applied Dynamical Systems 15, 1962 (2016).
  * Papachristodoulou and Prajna [2002] A. Papachristodoulou and S. Prajna, On the construction of Lyapunov functions using the sum of squares decomposition, in  _Proceedings of the 41st IEEE Conference on Decision and Control, 2002._ , Vol. 3 (IEEE, 2002) pp. 3482–3487.
  * Murty and Kabadi [1985] K. G. Murty and S. N. Kabadi,  _Some NP-complete problems in quadratic and nonlinear programming_ , Tech. Rep. (1985).
  * Parrilo [2003] P. A. Parrilo, Semidefinite programming relaxations for semialgebraic problems, Mathematical programming 96, 293 (2003).
  * Choi  _et al._ [1995] M.-D. Choi, T. Y. Lam, and B. Reznick, Sums of squares of real polynomials, in  _Proceedings of Symposia in Pure mathematics_ , Vol. 58 (American Mathematical Society, 1995) pp. 103–126.
  * Legat  _et al._ [2017] B. Legat, C. Coey, R. Deits, J. Huchette, and A. Perry, Sum-of-squares optimization in Julia, in  _The First Annual JuMP-dev Workshop_ (2017).
  * Weisser  _et al._ [2019] T. Weisser, B. Legat, C. Coey, L. Kapelevich, and J. P. Vielma, Polynomial and moment optimization in Julia and JuMP, in  _JuliaCon_ (2019).
  * MOSEK ApS [2023] MOSEK ApS,  _The MOSEK optimizer for C manual. Version 10.0._ (2023).
  * Tucker [1999] W. Tucker, The Lorenz attractor exists, Comptes Rendus de l’Académie des Sciences-Series I-Mathematics 328, 1197 (1999).
  * Figueras and de la Llave [2017] J.-L. Figueras and R. de la Llave, Numerical computations and computer assisted proofs of periodic orbits of the Kuramoto–Sivashinsky equation, SIAM Journal on Applied Dynamical Systems 16, 834 (2017).
  * Wilczak and Zgliczyński [2020] D. Wilczak and P. Zgliczyński, A geometric method for infinite-dimensional chaos: Symbolic dynamics for the Kuramoto-Sivashinsky PDE on the line, Journal of Differential Equations 269, 8509 (2020).
  * Henrion  _et al._ [2019] D. Henrion, S. Naldi, and M. Safey El Din, SPECTRA–a Maple library for solving linear matrix inequalities in exact arithmetic, Optimization Methods and Software 34, 62 (2019).
  * Peyrl and Parrilo [2008] H. Peyrl and P. A. Parrilo, Computing sum of squares decompositions with rational coefficients, Theoretical Computer Science 409, 269 (2008).
  * Cifuentes  _et al._ [2020] D. Cifuentes, T. Kahle, and P. Parrilo, Sums of squares in macaulay2, Journal of Software for Algebra and Geometry 10, 17 (2020).
  * Jansson  _et al._ [2008] C. Jansson, D. Chaykin, and C. Keil, Rigorous error bounds for the optimal value in semidefinite programming, SIAM Journal on Numerical Analysis 46, 180 (2008).
  * Jaulin  _et al._ [2001] L. Jaulin, M. Kieffer, O. Didrit, and Éric Walter,  _Applied Interval Analysis: With Examples in Parameter and State Estimation, Robust Control and Robotics_ (Springer London, 2001).
  * Tucker [2011] W. Tucker,  _Validated numerics: a short introduction to rigorous computations_ (Princeton University Press, 2011).
  * Sanders and Benet [2014] D. P. Sanders and L. Benet, IntervalArithmetic.jl (2014).
  * Goluskin and Fantuzzi [2019] D. Goluskin and G. Fantuzzi, Bounds on mean energy in the Kuramoto–Sivashinsky equation computed using semidefinite programming, Nonlinearity 32, 1705 (2019).
  * Fuentes  _et al._ [2022] F. Fuentes, D. Goluskin, and S. Chernyshenko, Global stability of fluid flows despite transient growth of energy, Physical Review Letters 128, 204502 (2022).