11institutetext: Department of Aeronautics, Imperial College London, London
SW7 2AZ, UK 22institutetext: The Alan Turing Institute, London NW1 2DB, UK
33institutetext: Isaac Newton Institute for Mathematical Sciences, Cambridge,
CB3 0EH, UK (visiting)
33email<EMAIL_ADDRESS><EMAIL_ADDRESS>
# Data-driven stability analysis of a chaotic time-delayed system ††thanks:
This research has received financial support from the ERC Starting Grant No.
PhyCo 949388.

Georgios Margazoglou 11 Luca Magri 11 2 2 3 3

###### Abstract

Systems with time-delayed chaotic dynamics are common in nature, from control
theory to aeronautical propulsion. The overarching objective of this paper is
to compute the stability properties of a chaotic dynamical system, which is
time-delayed. The stability analysis is based only on data. We employ the echo
state network (ESN), a type of recurrent neural network, and train it on
timeseries of a prototypical time-delayed nonlinear thermoacoustic system. By
running the trained ESN autonomously, we show that it can reproduce (i) the
long-term statistics of the thermoacoustic system’s variables, (ii) the
physical portion of the Lyapunov spectrum, and (iii) the statistics of the
finite-time Lyapunov exponents. This work opens up the possibility to infer
stability properties of time-delayed systems from experimental observations.

###### Keywords:

Echo State Networks Chaos Time-delayed systems.

##  1 Introduction

Chaotic systems with time-delayed dynamics appear in a range of scientific
fields [7]. Because of their dependence on both the present and past states,
these systems have rich and intricate dynamics. Their chaotic behaviour can be
assessed with stability analysis, which is a mathematical tool that quantifies
the system’s response to infinitesimal perturbations. Stability analysis
relies on the linearization of the time-delayed dynamical equations, which
spawns the Jacobian of the system. From the Jacobian, we compute the Lyapunov
Exponents (LEs), which are the key quantities to quantifying chaos [2].

A data-driven method, which considers the sequential nature of the dataset
(e.g. timeseries) to infer chaotic dynamics, is the recurrent neural network
(RNN). Such networks have been successfully applied to learn chaotic dynamics
(with no time delay) for different applications [12, 16, 11]. The majority of
RNNs require backpropagation through time for training, which can lead to
vanishing or exploding gradients, as well as long training times [16].
Instead, the echo state network (ESN), which we employ here, is trained via
ridge regression, which eliminates backpropagation and provides a faster
training [13, 6, 8]. The objective of this paper is to train an ESN with data
from a prototypical chaotic thermoacoustic system, which is a nonlinear time-
delayed wave equation. We further assess the capabilities of the ESN to
accurately learn the ergodic and stability properties of the thermoacoustic
system, by calculating fundamental quantities, such as the Lyapunov exponents.
We briefly review stability analysis for time-delayed systems in Sec. 2 and
present the considered thermoacoustic system in Sec. 2.1. In Sec. 3, we
discuss the ESN architecture and properties. We present the results in Sec. 4
and conclude in Sec. 5.

##  2 Stability Analysis for time-delayed systems

We consider a physical state ${\bm{x}}(t)\in\mathbb{R}^{D}$, which is the
solution of a nonlinear time-delayed dynamical system

$\frac{d{\bm{x}}}{dt}=f({\bm{x}}(t),\,{\bm{x}}(t-\tau)),\quad{\bm{x}}(t)={\bm{x}}_{0},\,\,\,\forall\,\,\,t\leq
0,$ (1)

where $\tau$ is a constant time-delay. We analyse the system’s stability by
perturbing the state with infinitesimal perturbations
${\bm{u}}\sim\mathcal{O}(\epsilon),\;\epsilon\to 0$, as ${\bm{x}}+{\bm{u}}$,
with ${\bm{x}}\sim\mathcal{O}(1)$. Hence, we obtain the tangent linear
equation

$\frac{d{\bm{U}}}{dt}={\bm{J}}\left({\bm{x}}(t),\,{\bm{x}}(t-\tau)\right){\bm{U}},$
(2)

which involves the time-marching of $K\leq D$ tangent vectors,
${\bm{u}}_{i}\in\mathbb{R}^{D}$, as columns of the matrix
${\bm{U}}\in\mathbb{R}^{D\times K}$,
${\bm{U}}=[{\bm{u}}_{1},{\bm{u}}_{2},\dots,{\bm{u}}_{K}]$. This is a linear
basis of the tangent space. The linear operator
${\bm{J}}\in\mathbb{R}^{D\times D}$ is the Jacobian of the system, which is
time-dependent in chaotic attractors. As shown in [11], we can extract the
Jacobian of a reservoir computer by linearizing Eqs. (7).

We periodically orthonormalize the tangent space basis during time evolution
by using a QR-decomposition of ${\bm{U}}$, as
${\bm{U}}(t)={\bm{Q}}(t){\bm{R}}(t,\Delta t)$ and by updating the columns of
${\bm{U}}$ with the columns of ${\bm{Q}}$, i.e. ${\bm{U}}\leftarrow{\bm{Q}}$
[2]. The matrix ${\bm{R}}(t,\Delta t)\in\mathbb{R}^{K\times K}$ is upper-
triangular and its diagonal elements $[{\bm{R}}]_{i,i}$ are the local growth
rates over a time span $\Delta t$ of ${\bm{U}}$. The LEs are the time averages
of the logarithms of the diagonal of $[{\bm{R}}]_{i,i}$, i.e.,

$\lambda_{i}=\lim\limits_{T\to\infty}\frac{1}{T}\int_{t_{0}}^{T}\ln[{\bm{R}}(t,\Delta
t)]_{i,i}dt.$ (3)

The FTLEs are defined as $\Lambda_{i}=\frac{1}{\Delta t}\ln[{\bm{R}}]_{i,i}$,
which quantify the expansion and contraction rates of the tangent space on
finite-time intervals, $\Delta t=t_{2}-t_{1}$.

###  2.1 Time-delayed thermoacoustic system

As a practical application of time-delayed systems, we consider a
thermoacoustic system, which is composed of three interacting subsystems, the
acoustics, the flame and the hydrodynamics (see, e.g., [9]). The interaction
of these sub systems can result in a positive feedback loop, which manifests
itself as a thermoacoustic instability. If uncontrolled, this instability can
lead to structural failure. We consider a prototypical time-delayed
thermoacoustic system with a longitudinal acoustic cavity and a heat source
modelled with a time-delayed model, following the same setup as in [10, 4, 5].
The system is governed by the conservation of momentum, mass, and energy. Upon
re-arrangement [9], thermoacoustic dynamics are governed by the nondimensional
partial differential equations

$\frac{\partial u}{\partial t}+\frac{\partial p}{\partial
x}=0,\qquad\frac{\partial p}{\partial t}+\frac{\partial u}{\partial x}+\zeta
p-\dot{q}\delta(x-x_{f})=0,$ (4)

where $u$, $p$, $\dot{q}$, $x\in[0,1]$ and $t$ are the non-dimensional
velocity, pressure, heat-release rate, axial coordinate and time,
respectively; and $\zeta$ is the damping coefficient, which takes into account
all the acoustic dissipation. The heat source is assumed to be small compared
to the acoustic wavelength, and it is modelled as a point in the grid, via the
Dirac delta distribution $\delta(x-x_{f})$, located at $x_{f}=0.2$. The heat-
release rate is provided by a modified King’s law,
$\dot{q}(t)=\beta\left(\sqrt{|1+u_{f}(t-\tau)|}-1\right)$, which is a
nonlinear time-delayed model. For the numerical studies of this paper, we set
$\beta=7.0$ for the heat parameter, and $\tau=0.2$ for the time delay. Those
values ensure chaotic evolution, and encapsulate all information about the
heat source, base velocity and ambient conditions.

As in [4], we transform the time-delayed problem into an initial value
problem. This is mathematically achieved by modelling the advection of a
perturbation $v$ with velocity $\tau^{-1}$ as

$\frac{\partial v}{\partial t}+\frac{1}{\tau}\frac{\partial v}{\partial
X}=0,\quad 0\leq X\leq 1,\quad v(X=0,t)=u_{f}(t).$ (5)

We discretise Eqs. (4) by a Galerkin method. First, we separate the acoustic
variables in time and space as $u(x,t)=\sum_{j=1}^{N_{g}}\eta_{j}(t)\cos(j\pi
x)$, and $p(x,t)=-\sum_{j=1}^{N_{g}}\mu_{j}(t)\sin(j\pi x)$, in which the
spatial functions are the acoustic eigenfunctions of the configuration under
investigation. Then, we project Eqs. (4) onto the Galerkin spatial basis
$\\{\cos(\pi x),$ $\cos(2\pi x),$ $\dots,\cos(N_{g}\pi x)\\}$ to obtain

$\dot{\eta}_{j}-j\pi\mu_{j}=0,\qquad\dot{\mu}_{j}+j\pi\eta_{j}+\zeta_{j}\mu_{j}+2\dot{q}\sin(j\pi
x_{f})=0.$ (6)

The system has $2N_{g}$ degrees of freedom. The time-delayed velocity becomes
$u_{f}(t-\tau)=\sum_{k=1}^{N_{g}}\eta_{k}(t-\tau)\cos(k\pi x_{f})$, and the
damping, $\zeta_{j}$, is modelled by $\zeta_{j}=c_{1}j^{2}+c_{2}j^{1/2}$,
where $c_{1}=0.1$ and $c_{2}=0.06$. The equation for linear advection, Eq.
(5), is discretised using $N_{c}+1$ points with a Chebyshev spectral method.
This discretisation adds $N_{c}$ degrees of freedom, thus a total of
$D=2N_{g}+N_{c}=30$ in our case, as $N_{g}=N_{c}=10$. We integrate Eqs. (6)
with a fourth order Runge-Kutta scheme and timestep $dt=0.01$.

##  3 Echo State Network

By applying the method of [11] to time-delayed problems, we linearize the Echo
State Network (ESN) [6] to calculate the stability properties of chaotic
systems. ESNs are proven effective for accurate learning of chaotic dynamics
(see e.g. [12, 16, 3, 1, 13, 11]). The ESN is a reservoir computer [6]. It has
a sparsely-connected single-layer hidden state, which is termed “reservoir”.
The reservoir weights, $\mathbf{W}$, as well as the input-to-reservoir
weights, $\mathbf{W}_{\mathrm{in}}$, are randomly assigned and remain fixed
through training and testing. The reservoir-to-output weights,
$\mathbf{W}_{\mathrm{out}}$, are trained via ridge regression. The evolution
equations of the reservoir and output are, respectively

${\bm{r}}({\rm{t}}_{i+1})=\tanh\left([\hat{{\mathbf{y}}}_{\mathrm{in}}({\rm{t}}_{i});b_{\mathrm{in}}]^{T}\mathbf{W}_{\mathrm{in}}+{\bm{r}}({\rm{t}}_{i})^{T}\mathbf{W}\right),\qquad{\mathbf{y}}_{\mathrm{p}}({\rm{t}}_{i+1})=[{\bm{r}}({\rm{t}}_{i+1});1]^{T}\mathbf{W}_{\mathrm{out}},$
(7)

where at any discrete time ${\rm{t}}_{i}$ the input vector,
${\mathbf{y}}_{\mathrm{in}}({\rm{t}}_{i})\in\mathbb{R}^{N_{y}}$, is mapped
into the reservoir state ${\bm{r}}\in\mathbb{R}^{N_{r}}$, by the input matrix,
$\mathbf{W}_{\mathrm{in}}$, where $N_{r}\gg N_{y}$[6, 13]. Here,
$\hat{(\;\;)}$ indicates normalization by the component-wise maximum-minus-
minimum range of the target in training set, T indicates matrix transposition,
and the semicolon indicates array concatenation. The dimensions of the weight
matrices are $\mathbf{W}_{\mathrm{in}}\in\mathbb{R}^{(N_{y}+1)\times N_{r}}$,
$\mathbf{W}\in\mathbb{R}^{N_{r}\times N_{r}}$ and
$\mathbf{W}_{\mathrm{out}}\in\mathbb{R}^{(N_{r}+1)\times N_{y}}$. The
hyperparameter input bias, $b_{\mathrm{in}}=1$, is selected to have the same
order of magnitude as the normalized inputs, $\mathbf{\hat{y}}_{\mathrm{in}}$.
The dimensions of the input and output vectors are equal to the dimension of
the dynamical system; here described by Eqs. (6), i.e. $N_{y}\equiv D$.
Furthermore, $\mathbf{W}_{\mathrm{out}}$ is trained via the minimization of
the mean square error
$\textrm{MSE}=\frac{1}{N_{\mathrm{tr}}N_{y}}\sum_{i=0}^{N_{\mathrm{tr}}}||{\mathbf{y}}_{\mathrm{p}}({\rm{t}}_{i})-{\mathbf{y}}_{\mathrm{in}}({\rm{t}}_{i})||^{2}$
between the outputs and the data over the training set, where $||\cdot||$ is
the $L_{2}$ norm, $N_{\mathrm{tr}}+1$ is the total number of data in the
training set, and ${\mathbf{y}}_{\mathrm{in}}$ the input data on which the ESN
is trained.

Training the ESN is performed by solving with respect to
$\mathbf{W}_{\mathrm{out}}$ via ridge regression of the equation
$(\mathbf{R}\mathbf{R}^{T}+\beta\mathbb{I})\mathbf{W}_{\mathrm{out}}=\mathbf{R}\mathbf{Y}_{\mathrm{d}}^{T}$.
In the previous expression $\mathbf{R}\in\mathbb{R}^{(N_{r}+1)\times
N_{\mathrm{tr}}}$ and $\mathbf{Y}_{\mathrm{d}}\in\mathbb{R}^{N_{y}\times
N_{\mathrm{tr}}}$ are the horizontal concatenation of the reservoir states
with bias, $[{\bm{r}}({\rm{t}}_{i});1]$,
${\rm{t}}_{i}\in[0,T_{\mathrm{train}}]$, and of the output data, respectively;
$\mathbb{I}$ is the identity matrix and $\beta$ is the Tikhonov regularization
parameter [15]. Therefore the ESN does not require backpropagation. The ESN
can run in two configurations, either open-loop or closed-loop. In open-loop,
which is necessary for the training stage, the input data is given at each
step, allowing for the calculation of the reservoir timeseries
${\bm{r}}({\rm{t}}_{i})$, ${\rm{t}}_{i}\in[0,T_{\mathrm{train}}]$. In closed-
loop the output ${\mathbf{y}}_{\mathrm{p}}$ at time step ${\rm{t}}_{i}$, is
recurrently used as an input at time step ${\rm{t}}_{i+1}$, allowing for the
autonomous temporal evolution of the network. The closed-loop configuration is
used for validation (i.e. hyperparameter tuning) and testing, but not for
training.

Regarding validation, we use the chaotic recycle validation (RVC), as
introduced in [13]. It has proven to be a robust strategy, providing enhanced
performance of the ESN, compared to standard strategies, as recently
successfully applied in [14, 11]. Briefly, in RVC the network is trained only
once on the entire training dataset (in open-loop), and validation is
performed on multiple intervals already used for training (but now in closed-
loop). The validation interval simply shifts as a small multiple of the first
Lyapunov exponent, $N_{\mathrm{val}}=3\lambda_{1}$ here. The key
hyperparameters that we tune are the input scaling $\sigma_{\mathrm{in}}$ of
the input matrix $\mathbf{W}_{\mathrm{in}}$, the spectral radius $\rho$ of the
matrix $\mathbf{W}$, and the Tikhonov parameter $\beta$. Furthermore,
$\sigma_{\mathrm{in}}$ and $\rho$ are tuned via Bayesian Optimization in the
hyperparameter space $[\sigma_{\mathrm{in}},\rho]=[0.1,5]\times[0.1,1]$ in
logarithmic scale, while for $\beta$ we perform a grid search
$\\{10^{-6},10^{-8},10^{-10},10^{-12}\\}$ within the optimal
$[\sigma_{\mathrm{in}},\rho]$. The reservoir size is $N_{r}=400$. The
connectivity of matrix $\mathbf{W}$ is set to $d=3$. We further add to the
training and validation data a Gaussian noise with zero mean and standard
deviation, $\sigma_{n}=0.0006\sigma_{y}$, where $\sigma_{y}$ is the standard
deviation of the data component-wise (noise regularizes the problem, see [16,
13, 11] for more details). The ESN is trained on a training set (open-loop) of
size $200\tau_{\lambda}$, and is tested on a test set (closed-loop) of size
$4000\tau_{\lambda}$, where $\tau_{\lambda}=1/\lambda_{1}$ is the Lyapunov
time, which is the inverse of the maximal Lyapunov exponent
$\lambda_{1}\approx 0.13$ in our case.

##  4 Results

We analyse the statistics produced by the autonomous temporal evolution of the
ESN and the target time-delayed system. The selected observables are the
statistics of the system’s chaotic variables, the Lyapunov exponents (LEs),
and the statistics of the finite-time Lyapunov exponents (FTLEs).

First, we test the capabilities of the ESN to learn the long-term statistical
properties of the thermoacoustic system by measuring the probability density
function (PDF) of the learned variables, ${\mathbf{y}}$. In Fig. 1, we show
the PDF of the first three components of the Galerkin modes (see Eq. (6))
$\eta_{i}$, $\mu_{i}$ for $i={1,2,3}$, in which the black line corresponds to
the target and the dashed red to the ESN. The ESN predictions are in agreement
with the target, including the variables that are not shown here.

Figure 1: Probability density functions (via histograms) of the three first
Galerkin modes, $\eta_{i}$, $\mu_{i}$, $i={1,2,3}$. Black is for target and
dashed red for ESN. The statistics for ESN are collected in autonomous
evolution on unseen data, after training and validation.

Testing the accuracy of the calculated LEs in autonomous evolution is a harder
consistency check for the ESN. Indeed, the ESN has been trained only on
timeseries of the variables $\eta_{i}$, $\mu_{i}$, and $v$. Therefore, a good
agreement of the LEs means that the ESN is capable to accurately reproduce
intrinsic chaotic properties of the system’s attractor. The LEs of the ESN are
calculated following [11], and in Fig. 2 we compare the first $K=14$ LEs. We
also add an inner plot showing, in additional detail, the first 6 LEs. Each LE
is the average of the measured LEs from five selected independent ESNs used
for the analysis. We train the ESNs independently on different chaotic target
sets. The shaded region corresponds to the standard deviation per
$\lambda_{i}$ from those five ESNs. There is close agreement for the first 8
exponents. In particular, we measure the leading, and only positive, exponent
$\lambda^{\rm targ}_{1}\approx 0.130$ and $\lambda^{\rm ESN}_{1}\approx 0.124$
for ESN, which gives a 4.7% absolute error. The ESN also provides an accurate
estimate of the neutral exponent ($\lambda_{2}=0$) with $\lambda^{\rm
ESN}_{2}\approx 0.008$. The rest of the exponents, $\lambda_{i},\,i\geq 3$,
are negative and the ESN achieves a small $8.2\%$ mean absolute percentage
error for all. Note that a gradual disagreement of the negative exponents,
which are sensitive due to the numerical method, between ESN and target has
also been reported in [12, 16] for the one-dimensional Kuramoto-Sivashinsky
equation.

Figure 2: The first 14 (outer plot) and the first 6 (inner plot) Lyapunov
exponents. The shaded red region indicates the error based on the ensemble of
5 ESNs. The statistics for ESN are collected in closed-loop mode.

Figure 3 shows the PDF of the first six FTLEs of the target (black line) and
ESN (dashed red line). We collect the statistics from 5 independent ESNs, thus
creating 5 histograms of FTLEs. We then average those histograms bin-wise, and
the standard deviation of each averaged bin is given by the shaded regions,
which are more pronounced at the tails. In Fig. 3 the plots in each column of
bottom row are identical to the ones in upper row, with the difference that
the y-axis is in logarithmic scale to emphasize the agreement of the
statistics also at the tails, as the agreement for the most probable
statistics close to the peak is good (upper row). Note that the mean of each
PDF should correspond to each Lyapunov exponent $\lambda_{i}$ (i.e. of Fig.
2), which is indeed the case.

Figure 3: Probability density functions of the first six finite-time Lyapunov
exponents. Black is for target and dashed red for ESN. The y-axis in bottom
row is in logarithmic scale.

We also report the values for the Kaplan-Yorke dimension for both ESN and
target. This dimension is an upper bound of the attractor’s fractal
dimension[2]. It is given by
$D_{\text{KY}}=k+\frac{\sum_{i=1}^{k}\lambda_{i}}{|\lambda_{i+1}|}$, where $k$
is such that the sum of the first $k$ LEs is positive and the sum of the first
$k+1$ LEs is negative. We obtain $D^{\rm targ}_{\text{KY}}\approx 3.37$ for
target and $D^{\rm ESN}_{\text{KY}}\approx 3.46$ for ESN, which results in a
2.7% absolute error. This observation further confirms the ability of the ESN
to accurately learn the properties of the chaotic attractor.

##  5 Conclusion

We propose a method to compute the stability properties of chaotic solutions
in time-delayed systems using only data. We use the echo state network (ESN)
as a surrogate model for learning the chaotic dynamics from the time series
observations of the system and determining its long-term statistical and
stability properties. By considering the ESN as a discrete dynamical system,
we linearize the map (7) to derive the tangent evolution of the attractor
through the Jacobian. When running the ESN in a long autonomous mode (closed-
loop), we show that (i) the long-term statistics of the variables are
correctly learned, (ii) the physical portion of the Lyapunov spectrum is
correctly predicted, and (iii) the finite-time Lyapunov exponents and their
statistics are correctly inferred. This work opens up the possibility to infer
the stability of nonlinear and time-delayed dynamical systems from data.

## References

  * [1] Doan, N.A.K., Polifke, W., Magri, L.: Short- and long-term predictions of chaotic flows and extreme events: a physics-constrained reservoir computing approach. Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences 477(2253), 20210135 (2021). https://doi.org/10.1098/rspa.2021.0135, https://royalsocietypublishing.org/doi/abs/10.1098/rspa.2021.0135
  * [2] Eckmann, J.P., Ruelle, D.: Ergodic theory of chaos and strange attractors. Rev. Mod. Phys. 57, 617–656 (Jul 1985). https://doi.org/10.1103/RevModPhys.57.617, https://link.aps.org/doi/10.1103/RevModPhys.57.617
  * [3] Huhn, F., Magri, L.: Learning ergodic averages in chaotic systems. In: Krzhizhanovskaya, V.V., Závodszky, G., Lees, M.H., Dongarra, J.J., Sloot, P.M.A., Brissos, S., Teixeira, J. (eds.) Computational Science – ICCS 2020. pp. 124–132. Springer International Publishing, Cham (2020)
  * [4] Huhn, F., Magri, L.: Stability, sensitivity and optimisation of chaotic acoustic oscillations. Journal of Fluid Mechanics 882, A24 (2020). https://doi.org/10.1017/jfm.2019.828
  * [5] Huhn, F., Magri, L.: Gradient-free optimization of chaotic acoustics with reservoir computing. Phys. Rev. Fluids 7, 014402 (Jan 2022). https://doi.org/10.1103/PhysRevFluids.7.014402, https://link.aps.org/doi/10.1103/PhysRevFluids.7.014402
  * [6] Jaeger, H., Haas, H.: Harnessing nonlinearity: Predicting chaotic systems and saving energy in wireless communication. Science 304(5667), 78–80 (2004). https://doi.org/10.1126/science.1091277, https://www.science.org/doi/abs/10.1126/science.1091277
  * [7] Lakshmanan, M., Senthilkumar, D.V.: Dynamics of nonlinear time-delay systems. Springer Science & Business Media (2011). https://doi.org/10.1007/978-3-642-14938-2, https://doi.org/10.1007/978-3-642-14938-2
  * [8] Lukoševičius, M.: A Practical Guide to Applying Echo State Networks, pp. 659–686. Springer Berlin Heidelberg, Berlin, Heidelberg (2012). https://doi.org/10.1007/978-3-642-35289-8_36, https://doi.org/10.1007/978-3-642-35289-8_36
  * [9] Magri, L.: Adjoint Methods as Design Tools in Thermoacoustics. Applied Mechanics Reviews 71(2) (03 2019). https://doi.org/10.1115/1.4042821, https://doi.org/10.1115/1.4042821, 020801
  * [10] Magri, L., Juniper, M.P.: Sensitivity analysis of a time-delayed thermo-acoustic system via an adjoint-based approach. Journal of Fluid Mechanics 719, 183–202 (2013). https://doi.org/10.1017/jfm.2012.639
  * [11] Margazoglou, G., Magri, L.: Stability analysis of chaotic systems from data. Nonlinear Dynamics pp. 1–21 (2023). https://doi.org/10.1007/s11071-023-08285-1, https://doi.org/10.1007/s11071-023-08285-1
  * [12] Pathak, J., Lu, Z., Hunt, B.R., Girvan, M., Ott, E.: Using machine learning to replicate chaotic attractors and calculate lyapunov exponents from data. Chaos: An Interdisciplinary Journal of Nonlinear Science 27(12), 121102 (2017). https://doi.org/10.1063/1.5010300, https://doi.org/10.1063/1.5010300
  * [13] Racca, A., Magri, L.: Robust optimization and validation of echo state networks for learning chaotic dynamics. Neural Networks 142, 252–268 (2021). https://doi.org/https://doi.org/10.1016/j.neunet.2021.05.004, https://www.sciencedirect.com/science/article/pii/S0893608021001969
  * [14] Racca, A., Magri, L.: Data-driven prediction and control of extreme events in a chaotic flow. Phys. Rev. Fluids 7, 104402 (Oct 2022). https://doi.org/10.1103/PhysRevFluids.7.104402, https://link.aps.org/doi/10.1103/PhysRevFluids.7.104402
  * [15] Tikhonov, A.N., Goncharsky, A., Stepanov, V., Yagola, A.G.: Numerical methods for the solution of ill-posed problems, vol. 328. Springer Science & Business Media (1995). https://doi.org/doi.org/10.1007/978-94-015-8480-7, https://doi.org/10.1007/978-94-015-8480-7
  * [16] Vlachas, P., Pathak, J., Hunt, B., Sapsis, T., Girvan, M., Ott, E., Koumoutsakos, P.: Backpropagation algorithms and Reservoir Computing in Recurrent Neural Networks for the forecasting of complex spatiotemporal dynamics. Neural Networks 126, 191–217 (2020). https://doi.org/https://doi.org/10.1016/j.neunet.2020.02.016, https://www.sciencedirect.com/science/article/pii/S0893608020300708


# Symplectic mean curvature flows in Kähler surfaces with positive holomorphic
sectional curvatures

Jiayu Li, Liuqing Yang  Jiayu Li, School of Mathematical Sciences, University
of Science and Technology of China Hefei 230026
AMSS CAS Beijing 100190, P. R. China  lijia@@amss.ac.cn Liuqing Yang, Academy
of Mathematics and Systems Sciences
Chinese Academy of Sciences
Beijing 100190, P. R. of China.  yangliuqing@@amss.ac.cn

###### Key words and phrases:

Symplectic mean curvature flow, Kähler angle, Holomorphic sectional curvature,
Holomorphic curve.

Abstract: In this paper, we mainly study the mean curvature flow in Kähler
surfaces with positive holomorphic sectional curvatures. We prove that if the
ratio of the maximum and the minimum of the holomorphic sectional curvatures
is less than $2$, then there exists a positive constant $\delta$ depending on
the ratio such that $\cos\alpha\geq\delta$ is preserved along the flow.

Mathematics Subject Classification (2000): 53C44 (primary), 53C21 (secondary).

##  1\. Introduction

Mean curvature flows were studied by many authors, for example Huisken ([16],
[17]), Ecker-Huisken ([7]), Huisken-Sinestrari [18], Carlo Ilmanen [19], Neves
[20], Smoczyk ([21]), Wang ([23]), White [24], etc.

In this paper we mainly concentrated on the symplectic mean curvature flows,
which were studied by Chen -Tian [5], Chen-Li [3], Chen-Li-Tian [4], Wang
[23], Han-Li [9], [10], [11], Han-Sun [18], Han-Li-Sun [13], and Han-Li-Yang
[14]. The basic fact is that the symplectic property is preserved by the mean
curvature flow if the ambient space $M$ is Kähler-Einstein, or if the ambient
Kähler surface evolves along the Kähler-Ricci flow [11].

Let $(M,J,\overline{\omega},\bar{g})$ be a Kähler surface. For a compact
oriented real surface $\Sigma$ which is smoothly immersed in $M$, the Kähler
angle [6] $\alpha$ of $\Sigma$ in $M$ was defined by

$\omega|_{\Sigma}=\cos\alpha d\mu_{\Sigma}$

where $d\mu_{\Sigma}$ is the area element of $\Sigma$ in the induced metric
from $\overline{g}$. We say that $\Sigma$ is a symplectic surface if
$\cos\alpha>0$; $\Sigma$ is a holomorphic curve if $\cos\alpha\equiv 1$.

Given an immersed $F_{0}:\Sigma\rightarrow M$, we consider a one-parameter
family of smooth maps $F_{t}=F(\cdot,t):\Sigma\rightarrow M$ with
corresponding images $\Sigma_{t}=F_{t}(\Sigma)$ immersed in $M$ and $F$
satisfies the mean curvature flow equation:

(1.1) $\begin{cases}\frac{\partial}{\partial t}F(x,t)=H(x,t)\\\
F(x,0)=F_{0}(x),\end{cases}$

where $H(x,t)$ is the mean curvature vector of $\Sigma_{t}$ at $F(x,t)$ in
$M$.

Choose an orthonormal basis $\\{e_{1},e_{2},e_{3},e_{4}\\}$ on ($M$,$\bar{g}$)
along $\Sigma_{t}$ such that $\\{e_{1},e_{2}\\}$ is the basis of $\Sigma_{t}$
and the symplectic form $\omega_{t}$ takes the form

(1.2) $\displaystyle\omega_{t}=\cos\alpha u_{1}\wedge u_{2}+\cos\alpha
u_{3}\wedge u_{4}+\sin\alpha u_{1}\wedge u_{3}-\sin\alpha u_{2}\wedge u_{4},$

where $\\{u_{1},u_{2},u_{3},u_{4}\\}$ is the dual basis of
$\\{e_{1},e_{2},e_{3},e_{4}\\}$. Then along the surface $\Sigma_{t}$ the
complex structure on $M$ takes the form ([3])

(1.3) $\displaystyle J=\begin{pmatrix}0&\cos\alpha&\sin\alpha&0\\\
-\cos\alpha&0&0&-\sin\alpha\\\ -\sin\alpha&0&0&\cos\alpha\\\
0&\sin\alpha&-\cos\alpha&0\end{pmatrix}.$

Recall the evolution equation of the Kähler angle along the mean curvature
flow deduced in [11],

######  Theorem 1.1.

The evolution equation for $\cos\alpha$ along $\Sigma_{t}$ is

(1.4) $\displaystyle(\frac{\partial}{\partial t}-\Delta)\cos\alpha$
$\displaystyle=$
$\displaystyle|\overline{\nabla}J_{\Sigma_{t}}|^{2}\cos\alpha+\sin^{2}\alpha
Ric(Je_{1},e_{2}).$

Here

(1.5)
$\displaystyle|\overline{\nabla}J_{\Sigma_{t}}|^{2}=|h^{4}_{1k}+h^{3}_{2k}|^{2}+|h^{4}_{2k}-h^{3}_{1k}|^{2}\geq\frac{1}{2}|H|^{2}.$

We want to see whether the symplectic property is preserved along the mean
curvature flow. In the case that $M$ is a Kähler-Einstein surface, we have
$Ric(Je_{1},e_{2})=\bar{\rho}\cos\alpha$, where $\bar{\rho}$ is the scalar
curvature of $M$, so the symplectic property is preserved. If the ambient
Kähler surface evolves along the Kähler-Ricci flow, Han-Li [11] derived the
evolution equation for $\cos\alpha$ and consequently they showed that the
symplectic property is also preserved. In this paper, we find another
condition to assure that along the flow, at each time the surface is
symplectic. Note that we don’t require $M$ to be Einstein. Denote the minimum
and maximum of holomorphic sectional curvatures of $M$ by $k_{1}$ and $k_{2}$.
We state our main theorem as follows:

Main Theorem Suppose $M$ is a Kähler surface with positive holomorphic
sectional curvatures. Set $\lambda=\frac{k_{2}}{k_{1}}$. If the flow satisfies
either

I. $1\leq\lambda<\frac{11}{7}$ and
$\cos\alpha(\cdot,0)\geq\delta>\frac{53(\lambda-1)}{\sqrt{(53\lambda-53)^{2}+(48-24\lambda)^{2}}}$,

or

II. $\frac{11}{7}\leq\lambda<2$ and
$\cos\alpha(\cdot,0)\geq\delta>\frac{8\lambda-5}{\sqrt{(8\lambda-5)^{2}+(12-6\lambda)^{2}}}$,

then along the flow

(1.6) $\displaystyle(\frac{\partial}{\partial
t}-\Delta)\cos\alpha\geq|\overline{\nabla}J_{\Sigma_{t}}|^{2}\cos\alpha+C\sin^{2}\alpha,$

where $C$ is a positive constant depending only on $k_{1}$, $k_{2}$ and
$\delta$. As a corollary, $\min_{\Sigma_{t}}\cos\alpha$ is increasing with
respect to $t$. In particular, at each time $t$, $\Sigma_{t}$ is symplectic.
Therefore, we call this flow the symplectic mean curvature flow.

Since we obtain (1.6), many theorems in “symplectic mean curvature flows in
Kähler-Einstein surfaces” still hold in our case. For example,

Arguing as in [6] by strong maximum principle, we have

######  Corollary 1.2.

I. Suppose $M$ is a Kähler surface with positive holomorphic sectional
curvatures and $1\leq\lambda<\frac{11}{7}$, then every symplectic minimal
surface satisfying

$\cos\alpha>\frac{53(\lambda-1)}{\sqrt{(53\lambda-53)^{2}+(48-24\lambda)^{2}}}$

in $M$ is a holomorphic curve.

II. Suppose $M$ is a Kähler surface with positive holomorphic sectional
curvatures and $\frac{11}{7}\leq\lambda<2$, then every symplectic minimal
surface satisfying

$\cos\alpha>\frac{8\lambda-5}{\sqrt{(8\lambda-5)^{2}+(12-6\lambda)^{2}}}$

in $M$ is a holomorphic curve.

Arguing exactly in the same way as in [3] or [23], we have

######  Theorem 1.3.

Under the same condition of the Main Theorem, the symplectic mean curvature
flow has no type I singularity at any $T>0$.

Acknowledgment: The research was supported by NSFC 11071236, 11131007 and
10421101.

##  2\. Curvature Tensor, Sectional Curvature and Holomorphic Sectional
Curvature

Denote the curvature tensor of $M$ by $K$. Set $K(X)=K(X,JX,X,JX)$ and
$K(X,Y)=K(X,Y,X,Y)$, where $X,Y$ are arbitrary vector fields on $M$. It is
known that (c.f. [1], [22]) we can express the sectional curvatures by
holomorphic sectional curvatures.

######  Theorem 2.1.

The sectional curvatures of $M$ can be determined by the holomorphic sectional
curvatures by

(2.1) $\displaystyle K(X,Y)$ $\displaystyle=$
$\displaystyle\frac{1}{32}[3K(X+JY)+3K(X-JY)-K(X+Y)-K(X-Y)$
$\displaystyle-4K(X)-4K(Y)].$

Using (2.1), it is easy to check that,

######  Theorem 2.2.

For any vector fields $X$, $Y$ and $Z$ on $M$,

(2.2) $\displaystyle K(X,Y,X,Z)$ $\displaystyle=$
$\displaystyle\frac{1}{2}[K(Y+Z,X)-K(X,Y)-K(X,Z)]$ $\displaystyle=$
$\displaystyle\frac{1}{64}[3K(Y+Z+JX)+3K(Y+Z-JX)-K(Y+Z+X)-K(Y+Z-X)$
$\displaystyle-3K(Y+JX)-3K(Y-JX)-3K(Z+JX)-3K(Z-JX)$
$\displaystyle-4K(Y+Z)+K(Y+X)+K(Y-X)+K(Z+X)+K(Z-X)$
$\displaystyle+4K(X)+4K(Y)+4K(Z)].$

Denote the minimum and the maximum of sectional curvatures by $K_{min}$ and
$K_{max}$ respectively, we have the following estimates.

######  Theorem 2.3.

$K_{min}$ and $K_{max}$ satisfy

(2.3) $\displaystyle K_{max}\leq\frac{3}{2}k_{2}-\frac{1}{2}k_{1}$

and

(2.4) $\displaystyle K_{min}\geq\frac{3}{4}k_{1}-\frac{1}{2}k_{2}$

Proof. Given any point $p\in M$ and any two unit orthogonal vectors $X$ and
$Y$ at $p$, we can find two vectors $Z$ and $W$ such that $\\{X,Y,Z,W\\}$ form
an orthonormal basis of $T_{p}M$. Suppose $JX=yY+zZ+wW$, then

(2.5) $\displaystyle{}\langle X+JY,X+JY\rangle=2-2y,$

and

(2.6) $\displaystyle{}\langle X-JY,X-JY\rangle=2+2y.$

Assume the Kähler form is anti-self-dual, it was shown in [14] that,
$y^{2}+z^{2}+w^{2}=1$ and $J$ has the form

(2.11) $\displaystyle J=\left(\begin{array}[]{clcr}0&y&z&w\\\ -y&0&w&-z\\\
-z&-w&0&y\\\ -w&z&-y&0\end{array}\right).$

Combining (2.1) with (2.5) and (2.6), we get

$\displaystyle K(X,Y)$ $\displaystyle\leq$
$\displaystyle\frac{1}{32}[3(2-2y)^{2}k_{2}+3(2+2y)^{2}k_{2}-2^{2}k_{1}-2^{2}k_{1}-4k_{1}-4k_{1}]$
$\displaystyle=$ $\displaystyle\frac{1}{4}[(3+3y^{2})k_{2}-2k_{1}]$
$\displaystyle\leq$ $\displaystyle\frac{1}{4}(6k_{2}-2k_{1})$ $\displaystyle=$
$\displaystyle\frac{3}{2}k_{2}-\frac{1}{2}k_{1},$

and similarly

$\displaystyle K(X,Y)$ $\displaystyle\geq$
$\displaystyle\frac{1}{4}[(3+3y^{2})k_{1}-2k_{2}]$ $\displaystyle\geq$
$\displaystyle\frac{1}{4}(3k_{1}-2k_{2})$ $\displaystyle=$
$\displaystyle\frac{3}{4}k_{1}-\frac{1}{2}k_{2}.$

This proves the theorem. Q.E.D.

##  3\. Proof of the Main Theorem

In this section, we will prove the Main Theorem of this paper.

Proof of the Main Theorem. In order to prove this theorem, we need to estimate
$Ric(Je_{1},e_{2})$. Using two different methods, we get two available
estimates. We now deduce the first one.

(3.1) $\displaystyle Ric(Je_{1},e_{2})$ $\displaystyle=$ $\displaystyle
K(Je_{1},e_{1},e_{2},e_{1})+K(Je_{1},e_{3},e_{2},e_{3})+K(Je_{1},e_{4},e_{2},e_{4})$
$\displaystyle=$ $\displaystyle K(\cos\alpha e_{2}+\sin\alpha
e_{3},e_{1},e_{2},e_{1})+K(\cos\alpha e_{2}+\sin\alpha
e_{3},e_{3},e_{2},e_{3})$ $\displaystyle+K(\cos\alpha e_{2}+\sin\alpha
e_{3},e_{4},e_{2},e_{4})$ $\displaystyle=$ $\displaystyle cos\alpha
R_{22}+\sin\alpha(K_{3121}+K_{3424}),$

where

(3.2) $\displaystyle R_{22}=K_{2121}+K_{2323}+K_{2424}.$

By (2.1), we have

$\displaystyle K_{2121}$ $\displaystyle=$
$\displaystyle\frac{1}{32}[3K(e_{1}+Je_{2})+3K(e_{1}-Je_{2})-K(e_{1}+e_{2})-K(e_{1}-e_{2})$
$\displaystyle-4K(e_{1})-4K(e_{2})].$

By our choice of the complex structure (1.3), we get

$\displaystyle\langle e_{1}+Je_{2},e_{1}+Je_{2}\rangle=2-2\cos\alpha,$

and

$\displaystyle\langle e_{1}-Je_{2},e_{1}-Je_{2}\rangle=2+2\cos\alpha.$

Hence $K_{2121}$ can be estimated by $k_{1}$ and $k_{2}$,

(3.3) $\displaystyle K_{2121}$ $\displaystyle\geq$
$\displaystyle\frac{1}{32}[3(2-2\cos\alpha)^{2}k_{1}+3(2+2\cos\alpha)^{2}k_{1}-2^{2}k_{2}-2^{2}k_{2}-4k_{2}-4k_{2}]$
$\displaystyle=$ $\displaystyle\frac{1}{4}[(3+3\cos^{2}\alpha)k_{1}-2k_{2}].$

Similarly, we get

(3.4) $\displaystyle K_{2323}\geq\frac{1}{4}(3k_{1}-2k_{2}),$

and

(3.5) $\displaystyle
K_{2424}\geq\frac{1}{4}[(3+3\sin^{2}\alpha)k_{1}-2k_{2}].$

Putting (3.3), (3.4) and (3.5) into (3.2), we obtain that

(3.6) $\displaystyle R_{22}\geq 3k_{1}-\frac{3}{2}k_{2}.$

Using (2.2) and (1.3), we can also estimate $K_{3121}$ and $K_{3424}$. We have

(3.7) $\displaystyle
K_{3121}\geq\frac{1}{32}[(53+48\sin\alpha\cos\alpha)k_{1}-53k_{2}],$

and

(3.8) $\displaystyle
K_{3424}\geq\frac{1}{32}[(53-48\sin\alpha\cos\alpha)k_{1}-53k_{2}].$

Adding (3.7) and (3.8) yields

(3.9) $\displaystyle K_{3121}+K_{3424}\geq-\frac{53}{16}(k_{2}-k_{1}).$

By a similar computation in the opposite direction, we get

(3.10) $\displaystyle|K_{3121}+K_{3424}|\leq\frac{53}{16}(k_{2}-k_{1}).$

Therefore by (3.1), (3.6), (3.10) and short time existence of the mean
curvature flow, we have

$\displaystyle Ric(Je_{1},e_{2})$ $\displaystyle\geq$
$\displaystyle\cos\alpha(3k_{1}-\frac{3}{2}k_{2})-\sqrt{1-\cos^{2}\alpha}\frac{53}{16}(k_{2}-k_{1})$
$\displaystyle=$
$\displaystyle(3\cos\alpha+\frac{53}{16}\sqrt{1-\cos^{2}\alpha})k_{1}-(\frac{3}{2}\cos\alpha+\frac{53}{16}\sqrt{1-\cos^{2}\alpha})k_{2}.$

If $1\leq\lambda<2$ and
$\cos\alpha>\frac{53(\lambda-1)}{\sqrt{(53\lambda-53)^{2}+(48-24\lambda)^{2}}}$,
then the RHS of (3) is positive.

Another estimate follows directly from Theorem 2.3 and Berger s inequality
(c.f. [8]) that

(3.12) $\displaystyle|K_{3121}+K_{3424}|\leq|K_{3121}|+|K_{3424}|\leq
K_{max}-K_{min}\leq 2k_{2}-\frac{5}{4}k_{1}.$

Putting the above estimate into $(\ref{Ric})$ yields

$\displaystyle Ric(Je_{1},e_{2})$ $\displaystyle\geq$
$\displaystyle\cos\alpha(3k_{1}-\frac{3}{2}k_{2})-\sqrt{1-\cos^{2}\alpha}(2k_{2}-\frac{5}{4}k_{1})$
$\displaystyle=$
$\displaystyle(3\cos\alpha+\frac{5}{4}\sqrt{1-\cos^{2}\alpha})k_{1}-(\frac{3}{2}\cos\alpha+\frac{5}{4}\sqrt{1-\cos^{2}\alpha})k_{2}.$

It follows that if $1\leq\lambda<2$ and
$\cos\alpha>\frac{8\lambda-5}{\sqrt{(8\lambda-5)^{2}+(12-6\lambda)^{2}}}$,
then the RHS of (3) is positive. Note that

$\displaystyle\frac{53(\lambda-1)}{\sqrt{(53\lambda-53)^{2}+(48-24\lambda)^{2}}}\leq\frac{8\lambda-5}{\sqrt{(8\lambda-5)^{2}+(12-6\lambda)^{2}}}$

for $1\leq\lambda<\frac{11}{7}$, and

$\displaystyle\frac{53(\lambda-1)}{\sqrt{(53\lambda-53)^{2}+(48-24\lambda)^{2}}}\geq\frac{8\lambda-5}{\sqrt{(8\lambda-5)^{2}+(12-6\lambda)^{2}}}$

for $\frac{11}{7}\leq\lambda<2$, we get the conclusion. Q.E.D.

## References

  * [1] M. Brozos-V zquez, E. Garc a-R o and P. Gilkey, Relating the curvature tensor and the complex Jacobi operator of an almost Hermitian manifold, Adv. Geom., 8 (2008), no. 3, 353-365.
  * [2] J. Cheeger and D. Ebin, Comparison Theorems in Riemannian Geometry, American Mathematical Society (1975).
  * [3] J. Chen and J. Li, Mean curvature flow of surfaces in 4-manifolds, Adv. Math., 163 (2001), 287-309.
  * [4] J. Chen, J. Li and G. Tian, Two-Dimensional graphs moving by mean curvature flow, Acta. Math. Sinica, English Series, 18 (2002), 209-224.
  * [5] J. Chen and G. Tian, Moving symplectic curves in Kähler-Einstein surfaces, Acta Math. Sinica, English Series, 16(2000), 541-548.
  * [6] S. S. Chern and J. Wolfson, Minimal surfaces by moving frams, Amer. J. Math., 105 (1983), 59-83.
  * [7] K. Ecker and G. Huisken, Mean curvature evolution of entire graphs, Ann. of Math., 130 (1989), 453-471.
  * [8] S. Goldberg, Curvature and Homology, Academic Press, London, 1962.
  * [9] X. Han and J. Li, The mean curvature flow approach to the symplectic isotopy problem, IMRN, 26 (2005), 1611-1620.
  * [10] X. Han, J. Li, Translating solitons to symplectic and Lagrangian mean curvature flows, Internat. J. Math. 20(4) (2009), 443 C458.
  * [11] X. Han and J. Li, The mean curvature flow along the Kähler-Ricci flow, arXiv: math.DG/.1105.1200v1.
  * [12] X. Han and J. Li, Simplectic critical surfaces in Kähler surfaces, J.Eur.Math.Soc.,12(2) (2010), 505-527.
  * [13] X. Han, J. Li and J. Sun, The second type singularities of symplectic and Lagrangian mean curvature flows, Chin. Ann. Math. Ser. B 32 (2011), no. 2, 223-240.
  * [14] X. Han,J. Li and L. Yang, symplectic mean curvature flow in $\textbf{C}P^{2}$, to appear in Cal. Var. PDE.
  * [15] X. Han, J. Sun, Translating solitons to symplectic mean curvature flows, Ann. Glob. Anal. Geom., 38 (2010), 161-169.
  * [16] G. Huisken, Flow by mean curvature of convex surfaces into spheres, J. Diff. Geom., 20(1984), 237-266.
  * [17] G. Huisken, Asymptotic behavior for singularities of the mean curvature flow, J. Diff. Geom. 31 (1990),
  * [18] G. Huisken and C. Sinestrari, Mean curvature flow with surgeries of two-convex hypersurfaces, Invent. Math. 175 (2009), 137 C221.
  * [19] T. Ilmanen, Elliptic Regularization and Partial Regularity for Motion by Mean Curvature , Mem. Am. Math. Soc., 520 (1994).
  * [20] A. Neves, Singularities of Lagrangian mean curvature flow: zero-Maslov class case, Invent. Math. 168(3) (2007), 449 C484.
  * [21] K. Smoczyk, Der Lagrangesche mittlere Kruemmungsfluss, Univ. Leipzig (Habil.-Schr.), 102 S. 2000.
  * [22] L. Vanhecke, Some almost Hermitian manifolds with constant holomorphic sectional curvature, J. Diff. Geom., 12 (1977), no. 4, 461-471.
  * [23] M.-T. Wang, Mean curvature flow of surfaces in Einstein four manifolds, J. Diff. Geom., 57 (2001), 301-338.
  * [24] B. White, A local regularity theoreom for mean curvature flow, Ann. Math. 161 2005, 1487-1519.