#  A fully adaptive explicit stabilized integrator for advection-diffusion-
reaction problems

Ibrahim Almuslimani*

###### Abstract

We introduce a novel second order family of explicit stabilized Runge-Kutta-
Chebyshev methods for advection-diffusion-reaction equations which outperforms
existing schemes for relatively high Peclet number due to its favorable
stability properties and explicitly available coefficients. The construction
of the new schemes is based on stabilization using second kind Chebyshev
polynomials first used in the construction of the stochastic integrator SK-
ROCK. We propose an adaptive algorithm to implement the new scheme that is
able to automatically select the suitable step size, number of stages, and
damping parameter at each integration step. Numerical experiments that
illustrate the efficiency of the new methods are presented.

Keywords:  advection-diffusion-reaction equations, explicit stabilized
methods, Runge-Kutta Chebyshev methods, RKC, SK-ROCK, ARKC.

AMS subject classification (2010):  65L04, 65L20, 65M12

11footnotetext: Univ Rennes, INRIA Rennes, IRMAR - UMR 6625, F-35000 Rennes,
France.
<EMAIL_ADDRESS>
##  1 Introduction

In this paper we use the idea of stabilization by combining first and second
kind Chebyshev polynomials introduced in [4] to derive explicit stabilized
methods for advection-diffusion problems with, possibly, costly non-stiff
reaction terms,

$\partial_{t}u(x,t)=\nabla\cdot(D\nabla
u(x,t))-\nabla\cdot(\textbf{v}u(x,t))+r(u(x,t)),\quad(x,t)\in\Omega\times[0,T],$

with initial and boundary conditions, where $\Omega\in\mathbb{R}^{d}$, $D$ is
the matrix of diffusion coefficients, and v is the velocity vector. The
function $r$ represents non-stiff, but possibly costly, reaction terms. Note
that in general, $D$ and v may also depend on $u$ leading to nonlinear
diffusion and advection terms. In the linear one dimensional setting, the
equation reduces to

$\partial_{t}u(x,t)=d\partial^{2}_{x}u(x,t)-a\partial_{x}u(x,t)+r(u(x,t)),\quad(x,t)\in\Omega\times[0,T]$
(1)

where $d$ and $a$ are positive reals, and $\Omega$ is a real interval. The
Peclet number is defined by $a/d$ and is allowed here to be quite large. When
discretizing the partial differential equation (PDE) (1) in space using
centered finite difference for example, with mesh size $\Delta x$, We obtain a
system of ordinary differential equations (ODEs) of the form

$\dot{y}(t)=F_{D}(y(t))+F_{A}(y(t)),\qquad y(0)=y_{0}\in\mathbb{R}^{d},\qquad
t\in[0,T],$ (2)

where $F_{D}$ represents the diffusion term with eigenvalues of its Jacobian
grow as $1/\Delta x^{2}$ on the negative real axis, and $F_{A}$ represents the
advection term (and possibly non-stiff reaction terms) with eigenvalues of its
Jacobian are of size $1/\Delta x$ and located close to the imaginary axis and
symmetric with respect to the origin. This means that the eigenvalues of the
Jacobian of the obtained system are approximately located in an ellipse with
the length of its minor axis proportional to the square root of the length of
the major axis.

Explicit stabilized Runge-Kutta-Chebyshev methods were originally introduced
in the context of purely diffusive or diffusion dominated advection-diffusion
problems (very small Peclet number) as a compromise between costly implicit
methods and restrictive usual explicit schemes [1, 2, 3, 6, 17]. Due to their
versatility, they were extended to many other types of problems such as
advection-diffusion-reaction equations [7, 18, 19, 20], stochastic
differential equations (SDEs) [4, 5, 8], and optimal control problems [10].

Typically, the stability domain of an explicit stabilized method contains a
long narrow strip around the negative real axis. In the context of advection-
diffusion problems, the authors of the article [19] propose the usage of the
RKC method with very large damping parameter to make the strip wider, which
means that eigenvalues with slightly larger imaginary parts coming from the
advection terms can be put in. This comes at the cost of a serious shortening
of the strip, which means that, for a fixed time step, less eigenvalues with
negative real parts can be put in. Later, a partitioned Ruge-Kutta-Chebyshev
method (PRKC) of order $2$ was designed in [20] based on the RKC method [15]
for the integration of ODEs that have a moderately stiff term (diffusion) and
non-stiff terms (advection or costly reaction terms). PRKC has a limited
stability for the advection term, and it shares with the standard RKC the same
stability domain length over the negative real axis. In [7], the authors
propose a partitioned implicit-explicit orthogonal Runge-Kutta method (called
PIROCK) for the time integration of advection-diffusion-reaction problems with
possibly severely stiff reaction terms and stiff stochastic terms. The
diffusion terms are solved by the explicit, nearly optimal, second order
orthogonal Chebyshev method (ROCK2). Applied to advection-diffusion problems,
the method has order 2 of accuracy and can handle the large Peclet number
regime but it needs very large damping that reduces a lot its stability domain
length over the negative real axis. In addition, PIROCK relies on the ROCK2
method, for which no explicit formulas are available to compute the
coefficients for a given stage number. This rich literature shows that the
domain of stabilized schemes is very active, and that the construction of an
adaptively efficient explicit stabilized integrator for such important class
of problems is quite challenging.

In [9, Sect. 3.7.2], the author of the present paper profited from the idea of
stabilization using second kind Chebyshev polynomials, first introduced in [4]
to construct a first order explicit stabilized method for advection-diffusion-
reaction equations with optimal stability domain. What made this quite
intuitive is the similarity between mean square stability for SDEs and the
stability of the test equation for ODEs of the form (2) (see Sect. 4.1). In
this work, we construct a second order integrator based on RKC that
outperforms the existing methods in the literature. We propose a fully
adaptive algorithm to implement the new second order method.

This paper is organized as follows: in Section 2, we give a fast revision on
explicit stabilized methods. In Section 3 we present an _optimal_ first order
explicit stabilized method for advection-diffusion-reaction equations that we
constructed in the thesis [9] inspired by previous work on SDEs [4]. In
Section 4 we derive and analyze our new second order adaptive scheme in terms
of stability and convergence, and we propose a local error estimator for
automatic step size selection. Section 5 is dedicated to present and analyze
some numerical experiments that illustrate the efficiency of the new schemes.

##  2 Preliminaries on explicit stabilized methods

This section is devoted to present useful standard materials.

In order to study the stability of a Runge-Kutta integrator applied to an ODE,
the following approach is widely used [3, 14]. Consider the test ODE

$\dot{y}(t)=\lambda y(t),\quad y(0)=y_{0},$ (3)

where $\lambda\in\mathbb{C}$ with negative real part. If we apply a Runge-
Kutta method with step size $h$ on (3), we get the relation
$y_{n+1}=R(h\lambda)^{n}y_{0}$, where the rational function $R(z)$ is called
the stability function. Hence, the stability domain of the method is defined
as

${\mathcal{S}}:=\\{z\in\mathbb{C};|R(z)|\leq 1\\}.$

In the particular case of explicit methods, $R(z)$ is a polynomial which means
that the stability domain is necessarily bounded. For example, the stability
domain of the explicit Euler method is just a disk of radius $1$ which
explains the time step restriction it faces for stiff ODEs, while that of the
implicit Euler method is the complementary of a disk of radius $1$. This shows
the advantage of implicit methods in terms of stability, however, for large
dimensional problems and especially the nonlinear and ill-conditioned ones,
implicit methods become very costly and difficult to implement. Here appears
the need to a compromise between classical explicit methods and implicit
integrators.

This compromise is ”explicit stabilized methods” (see the survey [3]). The
idea is to construct explicit Runge-Kutta integrators with extended stability
domain that grows quadratically with the number of stages $s$ of the method
along the negative real axis, and then allows to use large time steps
typically for problems arising from _diffusion dominant_ advection-diffusion-
reaction PDEs for which the eigenvalues are close to the negative real axis
and are very large in modulus.

Before proceeding, let us recall some useful facts on Chebyshev polynomials.
The first kind Chebyshev polynomials are defined by

$T_{0}(x)=1,\quad T_{1}(x)=x,\quad T_{j}(x)=2xT_{j-1}(x)-T_{j-2}(x),\quad
j\geq 2.$ (4)

The second kind Chebyshev polynomials are defined by

$U_{0}(x)=1,\quad U_{1}(x)=2x,\quad U_{j}(x)=2xU_{j-1}(x)-U_{j-2}(x),\quad
j\geq 2.$ (5)

Moreover, the two kinds polynomials satisfy the following

$U_{j-1}(x)=\frac{T_{j}^{\prime}(x)}{j}.$ (6)

The stabilization procedure is based on the above relations. The fact that
both kinds share the same recurrence relation will be very useful in our
analysis.

###  2.1 Optimal first order Chebyshev methods

Consider the ODE

$\dot{y}=f(y),\qquad y(0)=y_{0},\qquad t\in[0,T].$ (7)

Given $y_{0}$, in order to compute $y_{1}\approx y(h)$ using the optimal first
order Chebyshev method applied to (7) with step size $h$, the following
recurrence is applied

$\displaystyle K_{0}$ $\displaystyle=$ $\displaystyle y_{0},\quad
K_{1}=K_{0}+\mu_{1}hf(K_{0}),$ $\displaystyle K_{j}$ $\displaystyle=$
$\displaystyle\mu_{i}hf(K_{j-1})+\nu_{i}K_{j-1}+(1-\nu_{i})K_{j-2},\quad
j=2,\ldots,s$ (8) $\displaystyle y_{1}$ $\displaystyle=$ $\displaystyle
K_{s},$

where
$\omega_{0}:=1+\frac{\eta}{s^{2}},\,\omega_{1}:=\frac{T_{s}(\omega_{0})}{T^{\prime}_{s}(\omega_{0})}$,
and

$\mu_{1}:=\frac{\omega_{1}}{\omega_{0}},\quad\mu_{j}:=\frac{2\omega_{1}T_{j-1}(\omega_{0})}{T_{j}(\omega_{0})},\quad\nu_{j}:=\frac{2\omega_{0}T_{j-1}(\omega_{0})}{T_{j}(\omega_{0})},\quad
j=2,\ldots,s.$ (9)

The parameter $\eta$ is called the damping parameter and it is necessary to
avoid singularities in the stability domain which ensures the robustness of
the method. Typically for this method, $\eta$ is fixed to $0.05$. It can be
easily verified, using the recurrence (4) and proceeding by induction, that
applied to the test problem (3), the above method produces after one step
$y_{1}=K_{s}=R_{s}(h\lambda)y_{0}$ with

$R_{s}(z)=\frac{T_{s}(\omega_{0}+\omega_{1}z)}{T_{s}(\omega_{0})},$

for which the stability domain contains a narrow strip around the interval
$[-C_{\eta}s^{2},0]$ with
$C_{\eta}~{}=~{}\frac{(1-\omega_{0})}{\omega_{1}}\simeq 2-4/3\eta$ is very
close to $2$ (the optimal value for order 1). The method has low memory
requirements (only two stages have to be stored) and reasonable propagation of
round-off errors even for large values of $s$ needed in practice [17]. The
fact that the length of the stability domain on the negative real axis enjoys
a quadratic growth with respect to the number of stages $s$ is crucial to the
success of explicit stabilized Runge-Kutta methods.

###  2.2 Second order RKC methods

To design a second order method, we need the stability polynomial to
satisfy***Indeed, up to order two, the order conditions for nonlinear problems
are the same as the order conditions for linear problems [12, Chap. III].

$R(z)=1+z+\frac{z^{2}}{2}+\mathcal{O}(z^{3}).$

In [11], Bakker introduced a correction to the first order shifted Chebyshev
polynomials to get the following second order polynomial

$R_{s}(z)=a_{s}+b_{s}T_{s}(\omega_{0}+\omega_{2}z),$

where,

$a_{s}=1-b_{s}T_{s}(\omega_{0}),\quad
b_{s}=\frac{T^{\prime\prime}_{s}(\omega_{0})}{(T^{\prime}_{s}(\omega_{0})^{2})},\quad\omega_{0}=1+\frac{\eta}{s^{2}},\quad\omega_{2}=\frac{T^{\prime}_{s}(\omega_{0})}{T^{\prime\prime}_{s}(\omega_{0})},\quad\eta=0.15.$
(10)

For each $s$, $|R_{s}(z)|$ remains bounded by
$a_{s}+b_{s}=1-\eta/3+\mathcal{O}(\eta^{2})$ for $z$ in the stability interval
(except for a small interval near the origin). The stability interval along
the negative real axis is $[-\frac{1+\omega_{0}}{\omega_{2}},0]$ which is
approximately $[-0.65s^{2},0]$, and covers about $80\%$ of the optimal
stability interval for second order stability polynomials, and the formula now
for calculating $s$ for a given time step $h$ is

$s:=\left[\sqrt{\frac{h\lambda_{\max}+1.5}{0.65}}+0.5\right].$

Using the recurrence relation of the Chebyshev polynomials, the RKC method as
introduced in [17] is defined by

$\begin{split}K_{0}&=y_{0},\quad K_{1}=K_{0}+hb_{1}\omega_{2}f(K_{0}),\\\
K_{j}&=\mu_{j}h(f(K_{j-1})-a_{j-1}f(K_{0}))+\nu_{j}K_{j-1}+\kappa_{j}K_{j-2}+(1-\nu_{j}-\kappa_{j})K_{0},\\\
y_{1}&=K_{s},\end{split}$ (11)

where

$\mu_{j}=\frac{2b_{j}\omega_{2}}{b_{j-1}},~{}~{}\nu_{i}=\frac{2b_{j}\omega_{0}}{b_{j-1}},~{}~{}\kappa_{j}=-\frac{b_{j}}{b_{j-2}},~{}~{}b_{j}=\frac{T^{\prime\prime}_{j}(\omega_{0})}{T^{\prime}_{j}(\omega_{0})^{2}},~{}~{}a_{j}=1-b_{j}T_{j}(\omega_{0}),$
(12)

for $j=2,\ldots,s$. The parameters $b_{0}$ and $b_{1}$ are free ($R_{0}(z)$ is
constant and only order 1 is possible for $R_{1}$(z)) and the values
$b_{0}=b_{1}=b_{2}$ are suggested in [16].

######  Remark 2.1.

For simplicity of the presentation, we consider the case of autonomous
problems (with $f$ independent of time) but we highlight that our approach
also applies straightforwardly to non-autonomous problems
$\dot{y}(t)=f(t,u(t),y(t))$. Indeed, a standard approach is to consider the
augmented system with $z(t)=t$, i.e. $\frac{dz}{dt}=1,~{}z(0)=0$ and define
$\tilde{y}(t)=(y(t),z(t))^{T}$, see e.g. [12, Chap. III] for details.

##  3 Optimal first order scheme

Note that we have discussed the content of this section in the thesis [9], but
we recall it here since it gives insight about the construction of the second
order adaptive integrator in the next section.

Consider the linear test problem

$\dot{y}=\lambda y+i\mu y,\qquad y(0)=y_{0},$ (13)

where $\lambda\in\mathbb{R}^{-}$, $\mu\in\mathbb{R}$, and $i=\sqrt{-1}$.
Applying a Runge-Kutta method to the above equation, one gets an induction of
the form

$y_{n+1}=R(p,q)y_{n},$

with $p=h\lambda$ and $q=h\mu$. We define the stability domain of a Runge-
Kutta method applied to (13) by

$\mathcal{S}=\\{(p,q)\in\mathbb{R}^{2}~{};~{}|R(p,q)|\leq 1\\}.$

Equation (13) can be seen as the test equation for linear SDEs with the $i\mu$
replacing the noise. Hence, inspired by SK-ROCK [4], we consider the following
stability polynomial

$R(p,q)=A(p)+B(p)iq:=\frac{T_{s}(\omega_{0}+\omega_{1}p)}{T_{s}(\omega_{0})}+\frac{U_{s-1}(\omega_{0}+\omega_{1}p)}{U_{s-1}(\omega_{0})}(1+\frac{\omega_{1}}{2}p)iq,$
(14)

where $T_{s}$ and $U_{s}$ are the first and the second kind Chebyshev
polynomials of degree $s$ (the number of stages), and the coefficients
$\omega_{0}$ and $\omega_{1}$ are the same as for the Chebyshev method (2.1).
The stability condition $|R(p,q)|\leq 1$ is equivalent to
$A(p)^{2}+B(p)^{2}q^{2}\leq 1$ which is exactly the mean square stability
condition described in [4, Sect. 2]. By [4, Theorem 3.2] and [4, Remark 3.6],
for all $\eta>0$ and all $s\in\mathbb{N}$, $|R(p,q)|\leq 1$ for all
$p\in[-2\omega_{1}^{-1},0]$ and $|q|\leq\sqrt{-2p}$ (See Figure 1).

Figure 1: Stability domain of the new optimal first order method (15) in the
$p-q$ plane for $s=10$ and $\eta=0.05$. The dashed lines correspond to
$\pm\sqrt{-2p}$.

The new order one method for space discretized advection-diffusion-reaction
equations (2), is defined in the same way as the SK-ROCK method [4], by
replacing the noise term by the advection-reactions terms,

$\begin{split}K_{0}&=y_{0}\\\
K_{1}&=y_{0}+\mu_{1}hF_{D}(y_{0}+\nu_{1}hF_{A}(y_{0}))+\kappa_{1}hF_{A}(y_{0})\\\
K_{j}&=\mu_{j}hF_{D}(K_{j-1})+\nu_{j}K_{j-1}+\kappa_{j}K_{j-2},\quad
j=2,\dots,s,\\\ y_{1}&=K_{s},\end{split}$ (15)

where $\nu_{1}=s\omega_{1}/2$, $\kappa_{1}=s\omega_{1}/\omega_{0}$ and the
rest of the coefficients are identical to those defined in (9). Assuming
enough regularity on $F_{D}$ and $F_{A}$, the convergence proof is
straightforward and based on [4, Lemma 4.2].

Note that the method requires only 1 evaluation of the advection-reaction
terms $F_{A}$ per time step, while other existing methods for advection
diffusion problems such as PIROCK [7] and PRKC [20] require 3 and 4
evaluations of $F_{A}$ respectively per time step. It also outperform the
mentioned methods by far in terms of stability. However, it has only order 1
of accuracy, which motivates the construction of a second order integrator
that competes with existing schemes in terms of stability and convergence,
taking advantage of the above analysis.

Note that we can adaptively increase the stability for advection by increasing
$\eta$ in the price of loosing some length on the negative real axis. This
will be explained in details in the next section for the second order scheme,
for which this adaptive increase of damping is a key feature.

##  4 New second order scheme

In this section, we introduce the new adaptive second order integrator based
on the material presented in the first 3 sections and inspired by the SK-ROCK
method introduced in [4].

#### ARKC integrator

The ARKC integrator (for adaptive RKC) applied to (2) is defined as follows
for $s\geq 2$:

$\begin{split}G&=hF_{A}(y_{0}+\frac{h}{2}F_{A}(y_{0}+\frac{\omega_{2}}{2}hF_{D}(y_{0}))+\frac{h}{2}F_{D}(y_{0}))+hF_{D}(y_{0}+\frac{\omega_{2}-1}{2}hF_{A}(y_{0}))-hF_{D}(y_{0}),\\\
K_{-1}&=y_{0},\\\ K_{0}&=K_{-1}+\frac{\omega_{2}}{2}G,\\\
K_{1}&=K_{0}+b_{1}\omega_{2}hF_{D}(y_{0})+\alpha G,\quad\text{and for
}j=2,\dots,s,\\\
K_{j}&=\mu_{j}h(F_{D}(K_{j-1})-F_{D}(K_{0})+(1-a_{j-1})F_{D}(y_{0}))+\nu_{j}K_{j-1}+\kappa_{j}K_{j-2}+(1-\nu_{j}-\kappa_{j})K_{0},\\\
y_{1}&=K_{s},\end{split}$ (16)

where $\alpha=\left(1-\frac{\omega_{2}}{2}\right)b_{1}s\omega_{2},$ and the
other coefficients are identical to those defined in (10) and (12). Note that
for purely diffusive equations ($F_{A}\equiv 0$), the method reduces to
standard RKC (11).

#### Complexity

The method requires $s+2$ evaluations of $F_{D}$ and $3$ evaluations of
$F_{A}$ per time step. The standard RKC and the PRKC methods both require $s$
evaluations of $F_{D}$ per time step, but RKC requires $s$ evaluations of
$F_{A}$ while PRKC requires only $4$. Finally, PIROCK merhod needs $s+2+l$
evaluations of $F_{D}$ ($l=1$ or $2$) and $3$ evaluations of $F_{A}$ per time
step. The advantage of our new scheme comes from having larger stability
region that changes adaptively, which allows to use bigger time steps.

#### Construction

In order to construct a second order scheme, we need our stability function to
satisfy the following equality

$\begin{split}R_{2}(p,q)&=1+p+iq+\frac{1}{2}(p+iq)^{2}+\mathcal{O}(((p+iq)^{3})\\\
&=1+p+iq+\frac{p^{2}}{2}+ipq-\frac{q^{2}}{2}+\mathcal{O}((p+iq)^{3}).\end{split}$
(17)

Inspired by the previous section, we consider the following polynomial

$\begin{split}R_{2}(p,q)&=A_{2}(p)+B_{2}(p)(iq-\frac{q^{2}}{2})\\\
&:=a_{s}+b_{s}T_{s}(\omega_{0}+\omega_{2}p)+\left(\frac{\omega_{2}}{2}+\left(1-\frac{\omega_{2}}{2}\right)\frac{U_{s-1}(\omega_{0}+\omega_{2}p)}{U_{s-1}(\omega_{0})}\right)(1+\frac{\omega_{2}}{2}p)(iq-\frac{q^{2}}{2}),\end{split}$
(18)

where all the coefficients are defined in (10), from which we can derive the
new adaptive second order ARKC method (16) for advection-diffusion-reaction
problems, using the relations (4), (5), and (6).

Unlike the case of standard ODEs where the term $z^{2}/2$ is enough to have
order two for nonlinear problems, some additional _coupling_ conditions need
to be satisfied here due to the partitioned nature of the scheme (see [12,
Sect. III.2]). These conditions result from the fact that the term $ipq$ in
(17) is in fact the sum of two terms: $\frac{1}{2}ipq$ and $\frac{1}{2}iqp$,
which do not necessarily commute for multidimensional and nonlinear problems.
The quantity $G$ is constructed carefully to get order $2$ of convergence for
general nonlinear problems. Analogously to some existing second order methods
[7, 19], our damping parameter is not fixed to a small value, but it is an
increasing function of the stage number $s$. However, for relatively large
Peclet numbers, no huge damping is needed, and the method still perform very
well as will be shown in Section 5.

###  4.1 Stability analysis

######  Lemma 4.1.

The scheme (16) applied to the linear test problem (13) with step size $h$,
produces the following recurrence

$y_{n+1}=R_{s}(h\lambda,h\mu)y_{n},$

where the stability polynomial $R_{s}(p,q)$ is defined in (18).

###### Proof.

By induction on $j$, it can be shown that, for every $j\geq 1$, the internal
stages satisfy

$K_{j}=R_{j}(p,q)y_{n},$

where,

$R_{j}(p,q)\coloneqq
a_{j}+b_{j}T_{j}(\omega_{0}+\omega_{2}p)+\left(\frac{\omega_{2}}{2}+\left(1-\frac{\omega_{2}}{2}\right)U_{j-1}(\omega_{0}+\omega_{2}p)\right)b_{j}s\omega_{2}(1+\frac{\omega_{2}}{2}p)(iq-\frac{q^{2}}{2}).$

For $j=s$, we have that $b_{s}s\omega_{2}=1/U_{s-1}(\omega_{0})$ because
$s^{-1}T_{s}^{\prime}(\omega_{0})=U_{s-1}(\omega_{0})$ and the proof is done.
∎

In contrast to the first order scheme (15) and the stochastic integrator SK-
ROCK [4], the damping parameter for the ARKC method is not fixed. Indeed, it
is an increasing function of the Peclet number and the number of stages $s$.
We provide numerical stability analysis that illustrates the nice features of
the new scheme (see Figures 2 and 3).

(a) $\eta=0.15$, $d_{s}=0.65s^{2}$, $a_{s}=0.17s$.

(b) $\eta=1.5$, $d_{s}=0.56s^{2}$, $a_{s}=0.35s$.

(c) $\eta=3$, $d_{s}=0.5s^{2}$, $a_{s}=0.5s$.

(d) $\eta=10$, $d_{s}=0.35s^{2}$, $a_{s}=0.9s$.

Figure 2: Illustration of the stability regions of the ARKC method (16) in the
$p-q$ plane for $s=20$ and different values of the damping parameter $\eta$.

In what follows, we will respectively denote by $d_{s}$ and $a_{s}$ the half
width and the half height of the largest ellipse that can be put in the
stability region of the corresponding method.

The method (16) is designed to handle quite large Peclet numbers, this means
that we need the width of the stability region in the imaginary direction to
be as large as possible. While the length over the negative real axis grows
quadratically with $s$, we cannot achieve more than linear growth on the
imaginary direction (this fact is given as an exercise in [14, Sect. IV]). In
[19], the RKC method is proposed with large damping value which reduces the
length of the stability domain over the negative real axis to
$d_{s}=0.34s^{2}$, and leads to a growth of only $\mathcal{O}(\sqrt{s})$ for
the ellipse half-height $a_{s}$. The PRKC method proposed in [20] uses the
standard small damping $\eta=0.15$ that keeps $d_{s}\approx 0.65s^{2}$, but
$a_{s}$ is fixed to $1.7$ no matter how large is the number of stages $s$,
which does not add much to the standard RKC method since this might be useful
only in the small Peclet number regime. However, the main feature of the
method is to reduce the number of evaluations of possible non stiff terms such
as non-stiff advection or reaction terms. The PIROCK integrator [7], is
notably better than the two mentioned methods RKC and PRKC. Two kinds of
damping are proposed for PIROCK, the first lead to a nearly optimal length
over the negative real axis $d_{s}\approx 0.81s^{2}$ but $a_{s}$ is limited to
$0.07696s+1.878$. For the second damping, $a_{s}=0.5321s+0.4996$ but $d_{s}$
reduces to $0.43s^{2}$. It shares as well the feature of reducing the
evaluation of non diffusive terms. Nevertheless, PIROCK still have two
drawbacks: the first one is the non-availability of explicit closed form
formulas to compute its coefficients for a given number of stages $s$, the
other is the limited choice of damping.

The stability region of the ARKC method changes adaptively with the spectrum
of the Jacobian of the vector fields at each time step. This change is due to
the damping parameter $\eta$ which is not necessarily constant. In fact, the
ellipse half width $d_{s}$ still grows quadratically with $s$, while the
ellipse half height $a_{s}$ grows linearly with $s$, in contrast to RKC and
PRKC methods. For example, for the standard value of $\eta=2/13$ used for RKC
in [15] and PRKC in [20], $d_{s}$ is still equal to $0.65s^{2}$ and
$a_{s}\simeq 0.17s$ (see Figure 2(a)), whereas for PRKC $a_{s}$ is fixed to
$1.7$ . In addition, increasing the value of the damping parameter $\eta$ adds
more space in the imaginary direction which is very favorable for the
advection dominated problems. Figure 2(c) illustrates the stability region in
the p-q plane for $s=20$ and $\eta=3$. We can see that $d_{s}=0.5s^{2}$ and
$a_{s}=0.5s$, while the best that PIROCK [7] could achieve for almost the same
$a_{s}$ is $d_{s}=0.43s^{2}$ which makes difference for large values of $s$.
For $\eta=10$ we have $d_{s}=0.35s^{2}$ and $a_{s}=0.9s$, compare that with
the case of standard RKC with infinite damping, considered in [19], where
$d_{s}=0.34s^{2}$ while $a_{s}=\mathcal{O}(\sqrt{s})$.

###  4.2 Choice of damping

In fact, $a_{s}$ is a function of $s$ and $\eta$, this introduces additional
difficulty in the estimation of the value of $a_{s}$. Therefore, we will
introduce many choices of the range of Peclet number in order to simplify the
implementation. Let $\rho_{D}$ and $\rho_{A}$ be the spectral radii of the
Jacobians of $F_{D}$ and $F_{A}$ respectively. For $P_{e}\leq 0.1$, i.e,
$q\leq\sqrt{-p}/20$ (which corresponds to $\rho_{A}/\sqrt{\rho_{D}}\simeq
1/20$), we fix $\eta=0.15$ up to $s=200$, and for $s$ between $200$ and $500$
we set $\eta=0.6$ (we do not allow $s$ to be more than $500$). For
$P_{e}\simeq 1/2$ which means that $q\simeq\sqrt{-p}/4$
($\rho_{A}/\sqrt{\rho_{D}}\simeq 1/4$), we consider the following choices for
$s$ and $\eta$:

$2\leq s\leq 30$ | $31\leq s\leq 60$ | $61\leq s\leq 110$ | $111\leq s\leq 160$ | $161\leq s\leq 260$ | $261\leq s\leq 360$ | $361\leq s\leq 500$
---|---|---|---|---|---|---
$\eta=0.2$ | $\eta=0.45$ | $\eta=1$ | $\eta=1.5$ | $\eta=2.4$ | $\eta=3$ | $\eta=4$

Table 1: $q\simeq\sqrt{-p}/4$

Tables for other Peclet number regimes are given in appendix A. Infinitely
many choices could be done, but we will present and use only a few choices
since they are enough for the method to perform very well. The methodology to
compute the above values is easy, it is enough to plot $R_{s}(p,q)$ for the
corresponding choice of $q$ and to vary $\eta$ in a way that we stay stable
for the given value of $s$ (see Figure 3).

(a) $\eta=0.15$, $d_{s}=0.65s^{2}$.

(b) $\eta=1$, $d_{s}=0.59s^{2}$.

(c) $\eta=5$, $d_{s}=0.4345s^{2}$.

Figure 3: Stability polynomial (18) for $s=200$ (left) and $s=30$ (middle and
right), and different values of $\eta$ and $q$.

###  4.3 Convergence analysis

In this section we will prove the second order of convergence of the scheme
(16) when applied to ODEs of the form (2) arising from the discretization of
advection-diffusion-reaction problems.

######  Theorem 4.2.

Let $T>0$ and consider the system of ODEs (2) on the time interval $[0,T]$,
where $F_{D}$ and $F_{A}$ are of class $\mathcal{C}^{2}$ and are Lipschitz
continuous. Suppose in addition that the first and second derivatives of
$F_{D}$ and $F_{A}$ are bounded. Let $N\in\mathbb{N}$, $h=T/N$, and
$t_{n}=nh,~{}n=1,\dots,N$, then the method (16) applied to (2) with step size
$h$ has order $2$ of convergence. In other words, we have for all
$n=1\dots,N$,

$\|y(t_{n})-y_{n}\|\leq Ch^{2},$ (19)

where $C$ is independent of $h$ and $n$.

###### Proof.

Let us prove first that the local error (the error after one step) satisfies

$\|y(h)-y_{1}\|=\mathcal{O}(h^{3}).$ (20)

Throughout the proof, $F_{i}(y_{0})$ will be simply denoted by $F_{i}$, and
$F_{i}^{\prime}(y_{0})F_{j}(y_{0})$ will be denoted by $F_{i}^{\prime}F_{j}$
with $i,j\in\\{D,A\\}$.

Using Taylor expansion, we can easily see that

$G=hF_{A}+\frac{h^{2}}{2}F_{A}^{\prime}F_{A}+\frac{h^{2}}{2}F_{A}^{\prime}F_{D}+(\omega_{2}-1)\frac{h^{2}}{2}F_{D}^{\prime}F_{A}+\mathcal{O}(h^{3}).$

Now, let us suppose that

$K_{j}=y_{0}+\gamma_{1}^{j}hF_{D}+\gamma_{2}^{j}hF_{A}+\gamma_{3}^{j}h^{2}F_{D}^{\prime}F_{D}+\gamma_{4}^{j}h^{2}F_{D}^{\prime}F_{A}+\gamma_{5}^{j}h^{2}F_{A}^{\prime}F_{D}+\gamma_{6}^{j}h^{2}F_{A}^{\prime}F_{A}+\mathcal{O}(h^{3}).$
(21)

The first two stages of the method satisfy:

$K_{0}=y_{0}+\frac{\omega_{2}}{2}G,\quad
K_{1}=y_{0}+b_{1}\omega_{2}hF_{D}(y_{0})+(\alpha+\frac{\omega_{2}}{2})G,$

hence, we have, $\gamma_{1}^{0}=0$, $\gamma_{1}^{1}=b_{1}\omega_{2}$,
$\gamma_{2}^{0}=\frac{\omega_{2}}{2}$,
$\gamma_{2}^{1}=\alpha+\frac{\omega_{2}}{2}$,
$\gamma_{3}^{0}=\gamma_{3}^{1}=0$,
$\gamma_{4}^{0}=\frac{\omega_{2}(\omega_{2}-1)}{4}$,
$\gamma_{4}^{1}=(\alpha+\frac{\omega_{2}}{2})\frac{\omega_{2}-1}{2}$,
$\gamma_{5}^{0}=\frac{\omega_{2}}{4}$,
$\gamma_{5}^{1}=\frac{1}{2}(\alpha+\frac{\omega_{2}}{2})$,
$\gamma_{6}^{0}=\frac{\omega_{2}}{4}$, and finally,
$\gamma_{6}^{1}=\frac{1}{2}(\alpha+\frac{\omega_{2}}{2})$. By performing a
Taylor expansion of the stages $K_{j}$ defined in (16) and replacing $K_{j-1}$
and $K_{j-2}$ by the expansion defined in(21), and finally identifying the
coefficients, we get the following relations

$\begin{split}\gamma_{1}^{j}&=\mu_{j}(1-a_{j-1})+\nu_{j}\gamma_{1}^{j-1}+\kappa_{j}\gamma_{1}^{j-2},\\\
\gamma_{2}^{j}&=\nu_{j}\gamma_{2}^{j-1}+\kappa_{j}\gamma_{2}^{j-2}+(1-\nu_{j}-\kappa_{j})\frac{\omega_{2}}{2},\\\
\gamma_{3}^{j}&=\mu_{j}\gamma_{1}^{j-1}+\nu_{j}\gamma_{3}^{j-1}+\kappa_{j}\gamma_{3}^{j-2},\\\
\gamma_{4}^{j}&=\mu_{j}(\gamma_{2}^{j-1}-\frac{\omega_{2}}{2})+\nu_{j}\gamma_{4}^{j-1}+\kappa_{j}\gamma_{4}^{j-2}+(1-\nu_{j}-\kappa_{j})\frac{\omega_{2}(\omega_{2}-1)}{4},\\\
\gamma_{5}^{j}&=\nu_{j}\gamma_{5}^{j-1}+\kappa_{j}\gamma_{5}^{j-2}+(1-\nu_{j}-\kappa_{j})\frac{\omega_{2}}{4},\\\
\gamma_{6}^{j}&=\nu_{j}\gamma_{6}^{j-1}+\kappa_{j}\gamma_{6}^{j-2}+(1-\nu_{j}-\kappa_{j})\frac{\omega_{2}}{4}.\end{split}$

In order to prove (20), it is sufficient to show that
$\gamma_{1}^{s}=\gamma_{2}^{s}=1$ and
$\gamma_{3}^{s}=\gamma_{4}^{s}=\gamma_{5}^{s}=\frac{1}{2}$. Obviously,
$\gamma_{6}^{j}=\gamma_{5}^{j}$ for all $j\geq 0$. We will provide proofs for
the first 3 coefficients, the other two can be done using similar arguments.

  * •

For $\gamma_{1}^{j}$ we have: $\gamma_{1}^{0}=0$,
$\gamma_{1}^{1}=b_{1}\omega_{2}$, and
$\gamma_{1}^{j}=\mu_{j}(1-a_{j-1})+\nu_{j}\gamma_{1}^{j-1}+\kappa_{j}\gamma_{1}^{j-2}$,
thus $\gamma_{1}^{j}$ are the internal stages of the RKC method
$\eqref{eq:rkc}$ applied to the problem $\dot{y}=1,~{}y(0)=0$ with step size
$h=1$. This means that for all $0\leq j\leq s$, $\gamma_{1}^{j}=c_{j}$, where
$c_{j}$ is the $j^{th}$ node of the RKC method. Therefore,
$\gamma_{1}^{s}=c_{s}=1$. Moreover, according to [15, Sect. 2], for all
$j=2\dots,s$ we have

$\gamma_{1}^{j}=c_{j}=\omega_{2}\frac{T_{j}^{\prime\prime}(\omega_{0})}{T_{j}^{\prime}(\omega_{0})},\qquad
c_{1}=\frac{c_{2}}{T_{2}^{\prime}(\omega_{0})}=b_{1}\omega_{2},\quad c_{0}=0.$

  * •

It can be proved that $\gamma_{2}^{j},~{}j=2,\dots,s$ are given by
$\gamma_{2}^{j}=\frac{\omega_{2}}{2}+\alpha\frac{b_{j}}{b_{1}}P_{j}(\omega_{0}),$
where $P_{j}(x)$ are polynomials that satisfy the following two term
recurrence relation

$P_{0}(x)=0,\quad P_{1}(x)=1,\quad P_{j}(x)=2xP_{j-1}(x)-P_{j-2}(x),~{}j\geq
2.$

Comparing with the relation (5), it can be easily seen that for all $j\geq 1$,
$P_{j}(x)=U_{j-1}(x)=T_{j}^{\prime}(x)/j$. Hence,
$\gamma_{2}^{s}=\frac{\omega_{2}}{2}+\alpha\frac{b_{s}}{sb_{1}}T_{s}^{\prime}(\omega_{0})=\frac{\omega_{2}}{2}+(1-\frac{\omega_{2}}{2})\omega_{2}b_{s}T_{s}^{\prime}(\omega_{0})=1.$

  * •

The proof for $\gamma_{3}^{j}$ is very similar to that of $\gamma_{2}^{j}$.
Indeed, we can prove for each $j=0,\dots,s$ the equality
$\gamma_{3}^{j}=2\omega_{2}^{2}b_{j}Q_{j}(\omega_{0})$, where

$Q_{0}(x)=0,\quad Q_{1}(x)=0,\quad
Q_{j}(x)=T_{j-1}^{\prime}+2xQ_{j-1}(x)-Q_{j-2}(x),~{}j\geq 2.$

Using the relation (4), we observe that for all $j\geq 0$,
$Q_{j}(x)=\frac{T_{j}^{\prime\prime}(x)}{4}$, which implies that
$\gamma_{2}^{s}=2\omega_{2}^{2}b_{s}T_{s}^{\prime\prime}(\omega_{0})/4=1.$

Thus, (20) is proved, and using regularity assumptions made on the vector
fields, [13, Theorem 3.6] implies the global convergence estimate (19).

∎

###  4.4 Variable step size control and the fully adaptive algorithm

We introduce the following local error estimator that allows us to adaptively
select the time step size in order to reach a given accuracy,

$Est_{n+1}=C(12(y_{n}-y_{n+1})+6h(F_{D}(y_{n})+F_{A}(y_{n})+F_{D}(y_{n+1})+F_{A}(y_{n+1}))),$

where $C=1/6-c_{2}+(1/2-c_{1})\zeta-1/6\zeta$, with $\zeta=0$ if $F_{A}\equiv
0$ and $\zeta=1$ otherwise, and

$c_{1}=\frac{\omega_{2}}{2}\left(1-\frac{\omega_{2}}{2}\right)\left(1+\omega_{2}\frac{U^{\prime\prime}_{s-1}(\omega_{0})}{U_{s-1}(\omega_{0})}\right),\quad
c_{2}=sb_{s}U^{\prime\prime}_{s-1}(\omega_{0})\frac{\omega_{2}^{3}}{6}.$

To get an intuition about the above coefficients, compare the third order term
in the exact polynomial which is
$(p+iq)^{3}/6=p^{3}/6-iq^{3}/6+iqp^{2}/2-pq^{2}/2$ and the third order term of
stability polynomial (14) that is equal to $c_{2}p^{3}+c_{1}iqp^{2}-pq^{2}/2$.
Note that the above estimator is inspired by the one considered for RKC in the
paper [15], and they coincide for $\zeta=0$. In contrast to the estimators
introduced for PRKC [20] and PIROCK [7], we do not consider two separate
estimators for $F_{D}$ and $F_{A}$, since our method is defined in a different
way. Indeed, for improved stabilization, the advection-reaction terms are
computed at the beginning and not separately at the end, which makes the
method more similar to RKC. We adopt the standard step size selection strategy
proposed in [15] for RKC (see also [13, page 167]). Now, we are ready to
present our fully adaptive algorithm,

######  Algorithm 4.3 ($y_{0}\mapsto y_{1}$).

Given a time step size $h$:

  * •

Calculate $\rho_{A}/\sqrt{\rho_{D}}$ and choose the best table.

  * •

Search the minimum $s$ (and the corresponding $\eta$) in the chosen table such
that $\frac{1+\omega_{0}}{\omega_{1}}>h\rho_{D}$.

  * •

Generate the coefficients (10) and (12) and apply the recurrence (16).

  * •

Update the step size according to the automatic step size selection procedure.

The code of the ARKC integrator as well as the drivers that reproduce the
numerical experiments will be made publicly available on the page
https://sites.google.com/view/ibrahim-almuslimani.

##  5 Numerical experiments

We will compare ARKC with PIROCK and PRKC as they were shown to outperform the
other stabilized methods for advection-diffusion-reaction problems where the
reaction term is not stiff.

###  5.1 Linear 1D advection-diffusion problem

An excellent example to compare the performance of ARKC with PRKC and PIROCK
is the following 1-dimensional advection-diffusion equation with periodic
boundary conditions

$\begin{split}\partial_{t}u+a\partial_{x}u&=\partial_{x}^{2}u,\\\
u(x,0)&=\sin(2\pi x),\\\ u(0,t)&=u(1,t),\end{split}$ (22)

where $a$ is a positive constant ($P_{e}=a$), $x\in[0,1]$, and $t\geq 0$. We
discretize the space interval $[0,1]$ to a uniform grid
$\\{x_{k}\\}_{k=0}^{N}$ with $x_{k}=kh$, and $h=1/N$. We use second order
central differences for the advection and the diffusion terms. We denote by
$u_{k}(t)$ the approximation of $u(x_{k},t)$, and the periodic boundary
conditions propose that $u_{0}(t)=u_{N}(t)$. The eigenvalues of the obtained
matrix are

$\lambda_{k}=\frac{2}{h^{2}}(\cos(2k\pi h)-1)-\frac{ia}{h}\sin(2k\pi h),\quad
k=1,\dots,N,$

and are located in an ellipse in the left half-plane $\mathbb{C}^{-}$, which
makes the problem typical for the comparison of the three schemes.

Method | $a$ | $Steps$ |  $F_{D}$ evals |  $F_{A}$ evals | $s\\_{max}$ |  $L_{\infty}$ error at $t=1/2$
---|---|---|---|---|---|---
PRKC | $0.1$ | $14|71$ | $862|1860$ | $56|284$ | $126|106$ | $4.8\times 10^{-4}|3\times 10^{-7}$
PIROCK | $0.1$ | $13|237$ | $789|3654$ | $39|711$ | $150|104$ | $1.8\times 10^{-3}|7.5\times 10^{-7}$
ARKC | $0.1$ | $14|79$ | $886|2098$ | $42|237$ | $145|97$ | $4.3\times 10^{-4}|3.3\times 10^{-7}$
PRKC | $0.5$ | $27|77$ | $1341|2067$ | $108|308$ | $57|57$ | $8.7\times 10^{-9}|3.8\times 10^{-9}$
PIROCK | $0.5$ | $13|237$ | $789|3648$ | $39|711$ | $150|104$ | $1.8\times 10^{-3}|7.3\times 10^{-7}$
ARKC | $0.5$ | $13|79$ | $909|2132$ | $39|237$ | $142|117$ | $2.5\times 10^{-4}|2.2\times 10^{-7}$
PRKC | $1$ | $47|89$ | $1803|2342$ | $188|356$ | $40|40$ | $9\times 10^{-10}|6.1\times 10^{-10}$
PIROCK | $1$ | $13|237$ | $849|3648$ | $39|711$ | $200|104$ | $1.7\times 10^{-4}|6.6\times 10^{-7}$
ARKC | $1$ | $11|74$ | $896|2104$ | $33|222$ | $194|153$ | $2\times 10^{-4}|3.6\times 10^{-7}$
PRKC | $2$ | $90|120$ | $2575|2893$ | $360|480$ | $29|29$ | $1.2\times 10^{-10}|1\times 10^{-10}$
PIROCK | $2$ | $13|237$ | $934|3720$ | $39|711$ | $200|150$ | $2\times 10^{-6}2.5\times 10^{-8}$
ARKC | $2$ | $10|56$ | $995|2267$ | $30|168$ | $228|172$ | $4.8\times 10^{-5}|1.8\times 10^{-7}$
PRKC | $5$ | $222|229$ | $3980|4049$ | $888|916$ | $18|18$ | $4.1\times 10^{-11}|4\times 10^{-11}$
PIROCK | $5$ | $14|238$ | $1043|3826$ | $42|714$ | $182|150$ | $1.8\times 10^{-5}1.7\times 10^{-7}$
ARKC | $5$ | $12|59$ | $1272|2764$ | $36|177$ | $237|184$ | $1.9\times 10^{-6}|2.9\times 10^{-8}$
PRKC | $10$ | $442|453$ | $5738|5818$ | $1768|1812$ | $13|13$ | $6.2\times 10^{-11}|6\times 10^{-11}$
PIROCK | $10$ | $24|246$ | $1574|4086$ | $72|738$ | $125|125$ | $1.3\times 10^{-3}|1.2\times 10^{-5}$
ARKC | $10$ | $15|84$ | $1359|3207$ | $45|252$ | $234|160$ | $5.4\times 10^{-6}|7.3\times 10^{-8}$
PRKC | $12$ | $530|543$ | $6355|6436$ | $2120|2172$ | $12|12$ | $6.4\times 10^{-11}|6.2\times 10^{-11}$
PIROCK | $12$ | $28|255$ | $1750|4306$ | $84|765$ | $104|104$ | $1.1\times 10^{-2}|4.3\times 10^{-6}$
ARKC | $12$ | $18|104$ | $1557|3593$ | $54|312$ | $196|150$ | $3.5\times 10^{-5}|4.3\times 10^{-7}$
Table 2: Comparison of ARKC, PRKC, and PIROCK for the linear advection-
diffusion problem (22). The numbers on the left correspond to $tol=10^{-2}$,
while those on the right correspond to $tol=10^{-5}$.

For the numerical experiments, we take $N=150$, $t\in[0,1/2]$, and we fix
$Atol=Rtol=tol$. The initial step is fixed to $10^{-3}$. The number of
rejected steps is always very small and thus neglected. We can clearly see in
Table 2 that PRKC can compete with our scheme ARKC only in the very small
Peclet number regime, while For moderate and large Peclet number, ARKC is
notably better. On the other hand, the cost of ARKC is close to that of PIROCK
for moderate $P_{e}$ with small advantage for the latter when using large
tolerance $10^{-2}$ . However, ARKC becomes cheaper and more accurate at the
same time for large values of $P_{e}$ (in our experiment for $P_{e}=10$ and
$12$). For small tolerance $10^{-5}$ ARKC outperforms PIROCK in all Peclet
number regimes. This is expected because the damping (and so the vertical with
of the stability region) is fixed for PIROCK (when $F_{A}\neq 0$, PIROCK has
$d_{s}=0.43s^{2}$ and $a_{s}\simeq 0.53s+0.5$). The adaptivity of ARKC allows
to continuously increase the vertical width of its stability region which lets
it to be more flexible with respect to the change in $P_{e}$. One should not
forget that the explicit availability of ARKC coefficients helps a lot in
making the scheme adaptive with respect to the change in $P_{e}$, the feature
that is missing in PIROCK. The smaller number of steps needed in ARKC compared
to the two other methods, leads to a significantly lower number of $F_{A}$
evaluations no matter how big is the number of stages.

###  5.2 Burgers equation with a nonlinear reaction term

As an example of a PDE with variable Peclet number, we consider the following
Burgers equation with a nonlinear reaction term and periodic boundary
conditions

$\begin{split}\partial_{t}u+10u\partial_{x}u&=\partial^{2}_{x}u+\sin(u^{2}),\quad(x,t)\in[0,1]\times[0,1/2],\\\
u(x,0)&=1+\sin(2\pi x),\\\
u(0,t)&=u(1,t)\quad\forall\,t\in[0,1/2].\end{split}$ (23)

We discretize the above equation in space using second order central
differences with $\Delta x=10^{-2}$, into $N+1$ grid points with $N=100$. We
calculate a reference solution using the Radau IIA method of order 5 [14].

Figure 4: Change in Peclet number with respect to time in problem (23).

In Figure 4, we plot the Peclet number of equation (23) as a function of time.
we see that it is variable and of quite large magnitude.

Figure 5: Solution of equation (23) obtained using ARKC method (16) at
different time moments.

Figure 5 shows the solution $u(x,t)$ of equation (23) obtained using ARKC
method (16) at different time moments.

(a) $F_{D}$ evaluations

(b) $F_{A}$ evaluations

Figure 6: Comparison between ARKC, PRKC, and PIROCK in terms of cost for
problem (23).

In Figure 6, we compare the number of functions evaluations needed to obtain a
given accuracy for the solution of the Burgers equation (23) using 3 different
numerical methods: ARKC, PRKC, and PIROCK. The results are obtained for
$tol=10^{-r},\,r=1,\dots,6$. The advantage of our scheme ARKC (16) is very
clear. For PRKC, we get very close results for any value of $tol$ because of
the very high $P_{e}$ which oblige the method to severely restrict the time
step in order to stay stable. For PIROCK, the fixed damping increases the
number of steps and of functions evaluations with respect to ARKC. The
flexibility of ARKC gives it remarkable advantage over the other schemes. We
can also see that even when the number of $F_{D}$ evaluations for ARKC is
close to that of the other two schemes, the number of evaluations of $F_{A}$
containing the nonlinear advection and reaction terms stays much lower, that
is because of the low number of time steps needed (The number of evaluations
of $F_{A}$ is independent of the stage number for the 3 schemes).

#### Acknowledgment.

The author is grateful to Gilles Vilmart for helpful discussions and comments.
This work is supported by the Swiss National Science Foundation, project No:
P2GEP2_195212.

## References

  * [1] A. Abdulle.  Chebyshev methods based on orthogonal polynomials.  PhD Thesis, University of Geneva, Department of Mathematics. University of Geneva, 2001.
  * [2] A. Abdulle.  Fourth order Chebyshev methods with recurrence relation.  SIAM J. Sci. Comput., 23(6):2041–2054, 2002.
  * [3] A. Abdulle.  Explicit Stabilized Runge–Kutta Methods, pages 460–468.  Encyclopedia of Applied and Computational Mathematics, Springer Berlin Heidelberg, 2015.
  * [4] A. Abdulle, I. Almuslimani, and G. Vilmart.  Optimal explicit stabilized integrator of weak order 1 for stiff and ergodic stochastic differential equations.  SIAM/ASA J. Uncertain. Quantif., 6(2):937–964, 2018.
  * [5] A. Abdulle and T. Li.  S-ROCK methods for stiff Ito SDEs.  Commun. Math. Sci., 6(4):845–868, 2008.
  * [6] A. Abdulle and A. Medovikov.  Second order chebyshev methods based on orthogonal polynomials.  Numer. Math., 90(1):1–18, 2001.
  * [7] A. Abdulle and G. Vilmart.  PIROCK: a swiss-knife partitioned implicit-explicit orthogonal Runge-Kutta Chebyshev integrator for stiff diffusion-advection-reaction problems with or without noise.  J. Comput. Phys., 242:869–888, 2013.
  * [8] A. Abdulle, G. Vilmart, and K. C. Zygalakis.  Weak second order explicit stabilized methods for stiff stochastic differential equations.  SIAM J. Sci. Comput., 35(4):A1792–A1814, 2013.
  * [9] I. Almuslimani.  Explicit Stabilized Methods for Stiff Stochastic Differential Equations and Stiff Optimal Control Problems.  University of Geneva. PhD thesis, 2020.
  * [10] I. Almuslimani and G. Vilmart.  Explicit stabilized integrators for stiff optimal control problems.  SIAM J. Sci. Comput., 43(2):A721–A743, 2021.
  * [11] M. Bakker.  Analytical aspects of a minimax problem.  1971\.  Technical Note TN 62 (in Dutch), Mathematical centre, Amsterdam.
  * [12] E. Hairer, C. Lubich, and G. Wanner.  Geometric numerical integration, volume 31 of Springer Series in Computational Mathematics.  Springer-Verlag, Berlin, second edition, 2006.  Structure-preserving algorithms for ordinary differential equations.
  * [13] E. Hairer, S. Nørsett, and G. Wanner.  Solving Ordinary Differential Equations I. Nonstiff Problems, volume 8.  Springer Verlag Series in Comput. Math., Berlin, 1993.
  * [14] E. Hairer and G. Wanner.  Solving ordinary differential equations II. Stiff and differential-algebraic problems.  Springer-Verlag, Berlin and Heidelberg, 1996.
  * [15] B. Sommeijer, L. Shampine, and J. Verwer.  RKC: an explicit solver for parabolic PDEs.  J. Comput. Appl. Math., 88:316–326, 1998.
  * [16] B. P. Sommeijer and J. G. Verwer.  A performance evaluation of a class of Runge-Kutta-Chebyshev methods for solving semidiscrete parabolic differential equations.  Afdeling Numerieke Wiskunde [Department of Numerical Mathematics], 91\. Mathematisch Centrum, Amsterdam, 1980.
  * [17] P. J. van der Houwen and B. P. Sommeijer.  On the internal stability of explicit, $m$-stage Runge-Kutta methods for large $m$-values.  Z. Angew. Math. Mech., 60(10):479–485, 1980.
  * [18] J. G. Verwer and B. P. Sommeijer.  An implicit-explicit Runge-Kutta-Chebyshev scheme for diffusion-reaction equations.  SIAM J. Sci. Comput., 25(5):1824–1835, 2004.
  * [19] J. G. Verwer, B. P. Sommeijer, and W. Hundsdorfer.  RKC time-stepping for advection-diffusion-reaction problems.  J. Comput. Phys., 201(1):61–79, 2004.
  * [20] C. J. Zbinden.  Partitioned Runge-Kutta-Chebyshev methods for diffusion-advection-reaction problems.  SIAM J. Sci. Comput., 33(4):1707–1725, 2011.

##  Appendix A Damping and number of stages for some choices of Peclet number

For $P_{e}\simeq 1$ which means that $q\simeq\sqrt{-p}/2$
($\rho_{A}/\sqrt{\rho_{D}}\simeq 1/2$), we consider the following choices for
$s$ and $\eta$:

$s$ | $<11$ | $<21$ | $<31$ | $<41$ | $<51$ | $<61$ | $<71$ | $<81$ | $<91$ | $<101$
---|---|---|---|---|---|---|---|---|---|---
$\eta$ | $0.15$ | $0.6$ | $1$ | $1.4$ | $1.7$ | $2.1$ | $2.4$ | $2.7$ | $3$ | $3.3$
$s$ | $<121$ | $<141$ | $<161$ | $<181$ | $<201$ | $<251$ | $<301$ | $<401$ | $<501$ |
$\eta$ | $3.7$ | $4.1$ | $4.5$ | $4.9$ | $5.3$ | $6$ | $6.6$ | $7.7$ | $8.8$ |

Table 3: $q\simeq\sqrt{-p}/2$

For $P_{e}\simeq 3/2$ which means that $q\simeq 3\sqrt{-p}/4$
($\rho_{A}/\sqrt{\rho_{D}}\simeq 3/4$), we consider the following choices for
$s$ and $\eta$:

$s$ | $<11$ | $<21$ | $<31$ | $<41$ | $<51$ | $<61$ | $<71$ | $<81$
---|---|---|---|---|---|---|---|---
$\eta$ | $0.7$ | $1.5$ | $2.3$ | $2.9$ | $3.5$ | $4$ | $4.5$ | $4.9$
$s$ | $<91$ | $<101$ | $<141$ | $<181$ | $<251$ | $<301$ | $<401$ | $<501$
$\eta$ | $5.2$ | $5.5$ | $6.7$ | $7.7$ | $8.8$ | $9.8$ | $11$ | $12$

Table 4: $q\simeq 3\sqrt{-p}/4$

For $P_{e}\simeq 2$ which means that $q\simeq\sqrt{-p}$
($\rho_{A}/\sqrt{\rho_{D}}\simeq 1$), we consider the following choices for
$s$ and $\eta$:

$s$ | $<11$ | $<21$ | $<31$ | $51$ | $<71$ | $<111$ | $<151$ | $<311$ | $501$
---|---|---|---|---|---|---|---|---|---
$\eta$ | $1$ | $2.5$ | $3.5$ | $4.8$ | $6$ | $7.8$ | $9$ | $12.5$ | $15$

Table 5: $q\simeq\sqrt{-p}$

For $P_{e}\simeq 2\sqrt{2}$ which means that $q\simeq\sqrt{-2p}$
($\rho_{A}/\sqrt{\rho_{D}}\simeq\sqrt{2}$), we consider the following choices
for $s$ and $\eta$:

$s$ | $<11$ | $<21$ | $<31$ | $51$ | $<71$ | $<111$ | $<151$ | $<311$ | $501$
---|---|---|---|---|---|---|---|---|---
$\eta$ | $2$ | $3.8$ | $5$ | $6.8$ | $8$ | $10.4$ | $12$ | $16$ | $19$

Table 6: $q\simeq\sqrt{-p}$

For $P_{e}\geq 4$ which means that $q\geq 2\sqrt{-p}$
($\rho_{A}/\sqrt{\rho_{D}}\geq 2$), we consider the following choices for $s$
and $\eta$:

$s$ | $<11$ | $<31$ | $<71$ | $<151$ | $<311$ | $<501$
---|---|---|---|---|---|---
$\eta$ | $4$ | $9$ | $13.5$ | $18$ | $23$ | $27$

Table 7: $q\geq 2\sqrt{-p}$


# Identification of diffracted vortex beams at different propagation distances
using deep learning

Heng Lv1,a Yan Guo1,a Zi-Xiang Yang1 Chunling Ding1 Wu-Hao Cai1 Chenglong You2
<EMAIL_ADDRESS>Rui-Bo Jin1,3<EMAIL_ADDRESS>1Hubei Key Laboratory of Optical
Information and Pattern Recognition, Wuhan Institute of Technology, Wuhan
430205, China  2Quantum Photonics Laboratory, Department of Physics and
Astronomy, Louisiana State University, Baton Rouge, Louisiana, USA  3Guangdong
Provincial Key Laboratory of Quantum Science and Engineering, Southern
University of Science and Technology, Shenzhen, China  aThese authors
contributed equally to this work.

###### Abstract

Orbital angular momentum of light is regarded as a valuable resource in
quantum technology, especially in quantum communication and quantum sensing
and ranging. However, the OAM state of light is susceptible to undesirable
experimental conditions such as propagation distance and phase distortions,
which hinders the potential for the realistic implementation of relevant
technologies. In this article, we exploit an enhanced deep learning neural
network to identify different OAM modes of light at multiple propagation
distances with phase distortions. Specifically, our trained deep learning
neural network can efficiently identify the vortex beam’s topological charge
and propagation distance with 97% accuracy. Our technique has important
implications for OAM based communication and sensing protocols.

##  I Introduction

Vortex beam generally refers to the phase vortex beam, which has a spiral
wavefront, a phase singularity in the center of the beam, and ring-shaped
intensity distribution Franke-Arnold  _et al._ (2008); Rubinsztein-Dunlop  _et
al._ (2016). The beam with orbital angular momentum (OAM) has the phase term
$e^{i\ell\phi}$ in the complex amplitude equation, where $\phi$ is the
azimuthal angle and $\ell$ is the angular quantum number or topological
charge. OAM is an inherent characteristic of vortex beam photons, and each
photon carries OAM, which is $\ell\hbar$ Fickler  _et al._ (2012); Bai  _et
al._ (2022). Due to the high-dimensional characteristics of the photon OAM, it
is utilized in applications such as optical tweezers Friese  _et al._ (1998);
Padgett and Bowman (2011), micromanipulation Grier (2003); Curtis and Grier
(2003), angular velocity sensing Lavery  _et al._ (2013), quantum information
Mair  _et al._ (2001); Molina-Terriza  _et al._ (2007); Magaña-Loaiza and Boyd
(2019); Ding  _et al._ (2015), quantum computing Langford  _et al._ (2004);
Zhang  _et al._ (2007); Nagali  _et al._ (2009); Zhang  _et al._ (2012),
optical communications Bozinovic  _et al._ (2013); Wang  _et al._ (2012);
Ouyang  _et al._ (2021); Zhang  _et al._ (2021); Wang  _et al._ (2018) and
quantum cryptography Molina-Terriza  _et al._ (2001). Once the value of $\ell$
is identified, the orbital angular momentum can be calculated, allowing the
features of the vortex beam to be determined. Unfortunately, the vortex beam
will diffract during propagation, and its spatial profile will be easily
distorted in a real-world environment Rodenburg  _et al._ (2012).
Detrimentally, the information encoded in the structured beam can be destroyed
by random phase distortions Rodenburg  _et al._ (2014); Paterson (2005); Tyler
and Boyd (2009) and diffraction effects, resulting in mode loss and mode
cross-talk Ndagano  _et al._ (2017); Krenn  _et al._ (2014). As a result,
capturing the vortex beam and identifying its information using equipment such
as a Charge Coupled Device (CCD) camera or a Complementary Metal Oxide
Semiconductor (CMOS) camera is difficult Cox  _et al._ (2019). Hitherto, the
traditional methods of identifying vortex beams have included methods such as
the interferometer method, plane wave interferometry, and triangular aperture
diffraction measurement, to name a few Courtial  _et al._ (1998); Zhang  _et
al._ (2014); Padgett  _et al._ (1996); Yongxin  _et al._ (2011). These
traditional methods are much more challenging due to the need for more
equipment, as well as complicated data analysis process. In addition, some of
these methods can only identify specific vortex beams Karimi  _et al._ (2009).
Moreover, the accuracy of these methods will be greatly reduced when
turbulence is considered. These aforementioned factors have significantly
hampered the performance of communication, cryptography, and remote sensing.
As a result, identifying OAM efficiently and correctly while accounting for
diffraction and turbulence is a critical and unresolved challenge.

In recent years, methods such as deep learning algorithm Shin  _et al._ (2016)
and transfer learning Quattoni  _et al._ (2008) have considerably increased
the accuracy of automatic image recognition Taigman  _et al._ (2014); Melnikov
_et al._ (2018). A significant number of recent articles have proved the
potential of artificial neural networks for efficient pattern recognition and
spatial mode identification Doster and Watnik (2017); You  _et al._ (2020);
Bhusal  _et al._ (2021a), and its accuracy is far superior to some traditional
identification detection methods Gibson  _et al._ (2004); Liao  _et al._
(2021); Bhusal  _et al._ (2021b). However, due to the complex diffraction
effect in the OAM propagation process, there is little relevant work in the
identification of the propagation distance value. In the related research, the
propagation distance of the vortex beam ranges from the order of centimeters
to the order of kilometers, and it is used as a known parameter Zhang and Zhao
(2019); Krenn  _et al._ (2016). Different propagation distance $z$ will
drastically change the size of the central aperture of the vortex beam. As a
result, it remains difficult to identify the $z$ value using only the
intensity pattern. Additionally, changing the value of topological charge
$\ell$ also changes the size of the central aperture, making the
identification task more challenging. Finally, turbulence in real-world
applications exacerbates the difficulty of such an identification task Zhao
_et al._ (2011).

In this report, we take advantage of the deep learning algorithm to identify
vortex beams and their propagation distances while considering the effects of
undesired turbulence. Through theoretical simulations and experiments, we
generated vortex beams with different propagation distances and topological
charges. In addition, using the transfer learning method, we designed a deep
learning model to classify vortex beams. For the first time, our approach
utilizes artificial intelligence to simultaneously identify the propagation
distance and topological charge of a vortex beam under turbulence’s effects.
Our research enables the encoding of vortex beams with different propagation
distances. As a result, the vortex beam propagation distance may become a new
encoding variable. With the improvement of the accuracy of distance
recognition, it is even possible to realize precise distance measurement based
on the intensity of vortex beams. Our research opens up a new direction for
OAM communication and has great significance in OAM based sensing.

Figure 1: The experimental setup (upper panel) and our customized deep
learning algorithm (lower panel). We perform the experiment using a collimated
He-Ne laser beam. The vortex beam is generated using an SLM with the computer-
generated hologram. Finally, the intensity images are collected by a CCD and
used for training and testing. Our deep learning network consists of the
unaltered ResNet-101 bottom layer and our redesigned top layer.  Figure 2: The
spatial profiles of different vortex beams. We show the experimentally
measured images and simulated images for different topological charges $\ell$
and different propagation distances $z$, without turbulence in (a) and with
turbulence in (b). The first and third rows are the images acquired from the
experiment, and the second and fourth rows represent the theoretically
simulated images.

##  II Theory and Methods

###  II.1 Generation of the vortex beam

The fundamental beam used to produce the vortex beam is a Gaussian beam. By
applying a phase mask on the spatial light modulator (SLM), the beam amplitude
on the plane of SLM becomes Zhou  _et al._ (2016):

$E_{1}(r,\theta)=\sqrt{\frac{2}{\pi}}\exp(-\frac{r^{2}}{\omega^{2}_{0}})\exp(-i\ell\theta),$
(1)

where $\ell$ is the topological charge, $\omega_{0}$ is the Gaussian beam
waist, $r$ and $\theta$ are radial and azimuthal coordinates, respectively.

Within the framework of paraxial approximation, the field distribution of
$E_{1}(r,\theta)$ after propagation can be calculated using the Collins
integral equation Collins (1970):

$\displaystyle E_{2}\left(r_{1},\theta_{1},z\right)=\frac{i}{\lambda
B}\exp(-ikz)\int_{0}^{2\pi}\int_{0}^{\infty}E_{1}(r,\theta)$ (2)
$\displaystyle\times\exp\left[-\frac{ik}{2B}\left(Ar^{2}-2rr_{1}\cos\left(\theta_{1}-\theta\right)+Dr_{1}^{2}\right)\right]rdrd\theta,$

where $r_{1}$ and $\theta_{1}$ are radial and azimuthal coordinates in the
output plane, $z$ is the propagation distance, and $k=2\pi/\lambda$ is the
wave number with $\lambda$ being the wavelength. The ABCD transfer matrix of
light propagation in free space of distant $z$ is

$\left({\begin{array}[]{*{20}c}A&B\\\ C&D\\\
\end{array}}\right)=\left({\begin{array}[]{*{20}c}1&{z}\\\ 0&1\\\
\end{array}}\right).$ (3)

By inserting Eq. (1) and Eq. (3) into Eq. (2), we can obtain the beam
amplitude as

$\displaystyle{E_{2}}({r_{1}},{\theta_{1}},z)=\frac{{{i^{\ell+1}}\pi}}{{\lambda
z}}\exp\left({-ikz}\right)\exp\left({-\frac{{ikr_{1}^{2}}}{{2z}}}\right)$ (4)
$\displaystyle\times\exp\left({-i\ell{\theta_{1}}}\right)\frac{{b_{1}^{\ell}}}{{\varepsilon_{1}^{1+\frac{\ell}{2}}}}{\frac{{\Gamma\left({\frac{\ell}{2}+1}\right)}}{{\Gamma\left({\ell+1}\right)}}_{1}}{F_{1}}\left({\frac{{\ell+2}}{2},\ell+1,-\frac{{{b^{2}}}}{{{\varepsilon_{1}}}}}\right).$

Equation (4) represents the hypergeometric Gaussian mode.
${}_{1}F_{1}(\alpha,\beta,z)$ is a confluent hypergeometric function,
$\Gamma(n)$ is the Gamma function, $b_{1}$ and $\varepsilon_{1}$ are defined
as:

${b}_{1}=\frac{{kr_{1}}}{{2z}},\quad\varepsilon_{1}=\frac{1}{{\omega_{0}^{2}}}+\frac{{ik}}{{2z}}.$
(5)

Based on the above calculations, we can obtain transverse intensity images of
vortex beam with different values of $\ell$ after propagating different
distances $z$.

In actual communication, turbulence can lead to phase distortion of optical
mode spatial distribution. Therefore, in our experiment, we use the Kolmogorov
model with Von Karman spectrum of turbulence to simulate the atmospheric
turbulence in SLM to achieve a distorted communication mode Bos  _et al._
(2015); Glindemann  _et al._ (1993). The degree of distortion is quantified by
the Fried parameter $r_{0}$. The expression of the turbulence phase mask we
added on the SLM is Bhusal  _et al._ (2021a):

$\Phi(x,y)=\mathbb{R}\left\\{\mathcal{F}^{-1}\left(\mathbb{M}_{NN}\sqrt{\phi_{NN}(\kappa)}\right)\right\\},$
(6)

with
$\phi_{NN}(\kappa)=0.023r_{0}^{-5/3}(\kappa^{2}+\kappa_{0}^{2})^{-11/6}e^{-\kappa^{2}/\kappa_{m}^{2}}$
and the Fried’s parameter $r_{0}=(0.423k^{2}C_{n}^{2}z)^{-3/5}$. The symbol
$\mathbb{R}$ represents the real part of the complex field, and
$\mathcal{F}^{-1}$ indicates the inverse Fourier transform operation. In
addition, $\kappa$, $\kappa_{0}$, and $\mathbb{M}_{NN}$ represent the spatial
frequency, the central spatial frequency, and encoded random matrix,
respectively. $C_{n}^{2}$ is the standard refractive index, which is a
constant representing the turbulence intensity.

After the turbulence term is added to the original phase mask and loaded on
the SLM, beam amplitude on the plane of SLM becomes

$E_{1}^{\prime}(r,\theta)=E_{1}(r,\theta)\exp(i\Phi(x,y)).$ (7)

By substituting Eq.(7) into Eq.(2), we can numerically obtain the field
distribution of the turbulence distorted $E_{2}^{\prime}(r_{1},\theta_{1},z)$
after propagation.

The experimental setup and the deep learning model are shown in Fig. 1. In our
experiment, the vortex beams are generated by using an SLM and computer-
generated holograms. Utilizing the first-diffraction order of SLM, we can
obtain the vortex beam of arbitrary topological charge. A laser beam from a
He-Ne laser (wavelength of 632.8 nm) is coupled into a single-mode fiber for
spatial mode cleaning. A half-wave plate (HWP) and a quarter-wave plate (QWP)
are employed to adjust the polarization of the laser beam at the output port
of the fiber. An objective lens (magnification of 10$\times$ and an effective
focal length of 17 mm) is used to collimate the light from the fiber, and the
beam waist after collimation is around 2 mm. By loading a computer-generated
phase hologram onto the SLM, a Gaussian beam is converted into a vortex beam.
In order to simulate the turbulence in an atmosphere transmission process, we
can add an additional turbulence phase to the hologram. Finally, a CCD camera
is used to collect the intensity images of the vortex beam, and the
transmission distance is controlled by changing the distance between the CCD
and the SLM. The images collected by the CCD are sent to a computer for
training. Each training set, validation set, and test set contain 86, 10, and
10 images (360 $\times$ 360 pixels), in which the value of $\ell$ ranges from
1 to 5, and the propagation distance $z$ ranges from 40 cm to 100 cm with a
step of 5 cm. Totally, there are $86\times 5\times 13=5590$ images for the
training set and $10\times 5\times 13=650$ images for the validation set and
test set.

###  II.2 The deep learning algorithm

The lower panel of Fig. 1 shows our customized deep learning algorithm model.
Our model is a transfer learning network based on the ResNet-101 network
design He  _et al._ (2016). Since our obtained images have a high degree of
similarity, the neural network must have enough depth to extract image
features. Therefore, we adopt the CNN architecture and retrain the ResNet-101
deep learning model rather than the shallow neural networks model. More
specifically, the top layer is removed from the original ResNet-101. Moreover,
a global max pooling layer with a node count of 2048 is used to reduce the
parameters to increase the calculation speed. Following that, a dropout layer
is added to remove some parameters randomly to minimize over-fitting, and then
we use a fully connected (FC) layer to connect the local features. Another
dropout layer and an FC layer are added to lower the number of nodes from 1024
to 65. Finally, a softmax layer is applied for a 65 classification probability
output.

To train and test the deep learning model, we utilize a computer with an
Intel(R) Core(TM) i5-7300HQ CPU @2.5GHz and an Nvidia GeForce GTX 1050 Ti GPU
with 4GB of video memory. We use an adaptive moment estimation (Adam)
optimizer throughout the algorithm Kingma and Ba (2014). In our deep learning
model, we also used the transfer learning technique (TLT) Fernando  _et al._
(2014), which has two benefits. Firstly, it is highly efficient; for example,
tasks that originally required months of training without TLT can be reduced
to a few hours. The second merit of TLT is that less data is needed. This is
because transfer learning requires the use of a pre-trained model, which
allows us to achieve accurate recognition results with fewer datasets Patel
_et al._ (2015). Generally speaking, only hundreds or thousands of training
images (instead of tens of thousands or even millions of images) are enough to
achieve good training results Krizhevsky  _et al._ (2017); Liu  _et al._
(2019). Finally, the training results in each epoch are evaluated by the
categorical cross-entropy loss function Zhang and Sabuncu (2018) which is
given by

$\text{ Loss }=-\sum_{i=1}^{n}\left(\hat{y}_{i1}\ln y_{i1}+\hat{y}_{i2}\ln
y_{i2}+\ldots+\hat{y}_{im}\ln y_{im}\right),$ (8)

where $n$ is the number of samples, $m$ is the number of classifications,
$\hat{y}_{im}$ indicates that the true label (with the value of 0 or 1), and
$y_{im}$ is the predicted value of the $m$-th class given by the neural
network.

##  III Results and discussion

Figure 3: The cross-sectional view of the intensity images at $y=0$ for (a)
$\ell=3$, $z=70$ cm and (b) $\ell=4$, $z=50$ cm. The size of the central
apertures is comparable in both images. However, the number of side lobes in
the diffractive apertures is different. This subtle feature enables us to
distinguish these two cases. Figure 4: The accuracy and confusion matrix of
our trained deep learning algorithm. (a) The accuracy of the training set and
the validation set versus the epochs. The accuracy of up to 97% is achieved in
identifying vortex beams with different $\ell$ values and different $z$ values
after 90 epochs. All shaded areas correspond to the standard deviation of
accuracy. All classification results have been tested and verified, and some
test results are shown here. (b) The normalized confusion matrix between the
predicted propagation distance and the true propagation distance for $\ell=3$.
(c) Normalized confusion matrix between predicted $\ell$ values, predicted
propagation distance and true $\ell$ values, true propagation distance.

Figure 2 shows the spatial profiles of vortex beams with different propagation
distances $z$ and topological charges $\ell$ obtained from experiments and
simulations. Figure 2(a) shows the vortex beams without turbulence, and Fig.
2(b) depicts vortex beams affected by turbulence. The first and third rows
depict the spatial profiles of the vortex beams acquired in the experiment,
while the second and fourth rows depict the simulated ones under identical
conditions. For simplicity, we define the first aperture in the center as the
“center aperture” and other outer apertures as the “diffraction apertures”.
For a fixed value of $\ell$, the size of the central aperture becomes larger
as the propagation distance $z$ value increases. At the same time, for a fixed
propagation distance $z$, the size of the central aperture also increases as
the value of $\ell$ increases. From Fig. 2, we can observe that the spatial
profiles obtained from the experiment match well with our simulations,
therefore validating our theoretical model of the experiment. By comparing the
experimental and theoretical images of $\ell=4$ and $z=100$ cm, we can notice
that the central aperture in the experimental image is not distributed
uniformly. This effect is due to the slight misalignment of the collimated
beam to the center of the SLM. As a result, the brightness and shape of the
center aperture can vary slightly. This kind of deviation is also included, in
order to increase the diversities of the training data. We will show later
that, even with such deviations, the training results remain excellent.

There are many situations where the size of the central aperture is comparable
for vortex beams with different $\ell$ and propagation distance $z$. This
particular effect makes the simultaneous identification of $\ell$ and $z$
difficult. For example, by comparing the theoretical image of $\ell=3$ and
$z=70$ cm (second row, second column), and the image of $\ell=4$ and $z=50$ cm
(fourth row, first column) in Fig. 2(a), we can notice that the sizes of the
center aperture are similar, making it difficult to distinguish between these
two modes. In this case, the difference of the diffraction apertures provides
the best characteristic value to distinguish them. More intuitively, Fig. 3
shows the cross-sectional view of the intensity for these two beams at $y=0$.
It is clear that the distance between the two main peaks in Fig. 3(a) and Fig.
3(b) is almost the same. However, Fig. 3(a) has 4 side lobes in the
diffraction apertures, while Fig. 3(b) has 6 side lobes in the diffraction
apertures. This subtle difference makes it possible to distinguish these two
cases. Finally, we note that in Fig. 3 the diffraction apertures in the region
of $|x|>2.2$ mm (in black) usually cannot be well captured by a CCD in the
experiment, due to low light intensity and limited resolution of the CCD.

In practical communication applications, the spatial profile of the vortex
beam might be distorted due to atmospheric turbulence, underwater turbulence,
or other adverse circumstances. Therefore, the turbulence should be taken into
account for the propagation of the vortex beam. Figure 2(b) shows the typical
simulated and experimental diagrams with turbulence. In these theoretical (the
second and fourth rows) and experimental (the first and third rows) data, a
turbulence intensity parameter of $C_{n}^{2}=5\times 10^{-10}$ mm-2/3 is
utilized. In these data sets, we also take into account the fact that light
might not be precisely incident on the center of the SLM plane. It can be
noticed that the turbulence created huge distortions on the vortex beam’s
center aperture and diffraction apertures. For deep learning training, we
gathered 1040 distorted light intensity images.

To show the performance of our deep learning network, we plot the accuracy as
a function of the training epochs in Fig. 4(a). We can see that after 90
training epochs, the accuracy is higher than 92% for the training set, while
the accuracy of the validation set is greater than 97%. Since we added
regularization and dropout operations during the training process. These
operations will be automatically closed during the verification process,
causing the accuracy of the validation set to be higher than the accuracy of
the training set. A high validation accuracy of 97% indicates that our
approach provides a powerful way to identify the vortex beams with different
$\ell$ values and different $z$ values, even under a turbulent environment.
Finally, we note that the number of epochs required for convergence depends on
multiple factors, including the number of cases of different vortex beams
propagated and the degree of turbulence.

To show our results more comprehensibly, we calculated the normalized
confusion matrix for different $\ell$ and $z$. Figure 4(b) shows a typical
normalized confusion matrix from $\ell=3$, $z=40$ cm to $\ell=3$, $z=100$ cm.
Figure 4(c) shows the normalized confusion matrix from $\ell=1$, $z=40$ cm to
$\ell=5$, $z=45$ cm. The true propagation distance and the predicted
propagation distance given by our deep learning algorithm are basically on a
diagonal line. The result means almost all OAM modes and propagation distances
tested are correctly identified, only two images with $\ell=3$, $z=80$ cm are
predicted to be $\ell=3$, $z=75$ cm in Fig. 4(b), and only two images with
$\ell=2$, $z=40$ cm are predicted to be $\ell=2$, $z=45$ cm in Fig. 4(c).

Finally, we want to emphasize that using classical methods (e.g.,
interferometer) to analyze the distorted intensity images in Fig. 2(b) is
quite challenging. However, according to training and test results, our deep
learning model can accurately identify vortex beams with varying topological
charges and propagation distances even under the influence of severe
turbulence. This demonstrates that our approach has a high level of robustness
and is very useful for practical applications. We note that our approach can
be adapted to identify larger $\ell$ value with longer and more accurate
transmission distance. However, due to the limitation of our equipment, such
as the resolution of SLM and CCD, as well as the experimental error caused by
the laboratory environment, we limit the size of our topological charges and
the length of the propagation distance. We believe our designed deep learning
neural network does not fundamentally limit the recognition accuracy, and its
potential is far from being reached. Moreover, our scheme can be adapted to
many vortex beam related applications. For instance, we can adapt our work to
consider multiple types of vortex beams, and even the combination of them.
Furthermore, the accurate identification of the propagation distance might be
a novel technique for sensing related applications. Last but not least, our
approach can be applied to free-space OAM communication, especially the
demodulation system, to increase the robustness of the communication. We
expect that by combining the unique characteristics of vortex beams with the
advantages of the deep learning algorithm, more breakthroughs in vortex beams
research can be made in the future.

##  IV Conclusion

Vortex beams have enormous potential due to their versatility and virtually
unlimited quantum information resources. However, these beams are highly
susceptible to undesirable experimental conditions such as propagation
distance and phase distortions. In our work, we exploit the deep learning
algorithm to identify a vortex beam’s topological charge and propagation
distance. Specifically, we focus on vortex beam with topological charge $\ell$
from 1 to 5, and the propagation distance $z$ ranges from 40 cm to 100 cm.
Additionally, we consider the effect of turbulence-induced in the propagation
of the beam. We experimentally demonstrated that our customized deep learning
algorithm could accurately identify the propagation distance and topological
charge. Our work has important implications for the realistic implementation
of OAM-based optical communications and sensing protocols in a turbulent
environment.

## Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any
commercial or financial relationships that could be construed as a potential
conflict of interest.

## Author Contributions

The idea was conceived by R.-B.J., C.Y., and H.L. The experiment was designed
by H.L., Y.G., C.Y., and R.-B.J. The experiment was performed by H.L. with
help of Y.G., Z.X.Y, C.D., W.-H.C., C.Y. and R.-B.J. The theoretical
description and numerical simulation were developed by Y.G., C.Y., and R.-B.J.
The deep learning process was carried out by H.L. The data was analyzed by
H.L., Y.G., C.Y., and R.-B. J. The project was supervised by R.-B.J. and C.Y.
All authors contributed to the preparation of the manuscript.

## Funding

This work is supported by the National Natural Science Foundations of China
(Grant Nos.12074299, 91836102,11704290) and by the Guangdong Provincial Key
Laboratory (Grant No. GKLQSE202102).

## Acknowledgments

We thank Prof. Zhi-Yuan Zhou for the helpful discussion.

## References

  * Franke-Arnold  _et al._ (2008) S. Franke-Arnold, L. Allen, and M. Padgett, “Advances in optical angular momentum,” Laser & Photonics Review 2, 299–313 (2008).
  * Rubinsztein-Dunlop  _et al._ (2016) Halina Rubinsztein-Dunlop, Andrew Forbes, M V Berry, M R Dennis, David L Andrews, Masud Mansuripur, Cornelia Denz, Christina Alpmann, Peter Banzer, Thomas Bauer, Ebrahim Karimi, Lorenzo Marrucci, Miles Padgett, Monika Ritsch-Marte, Natalia M Litchinitser, Nicholas P Bigelow, C Rosales-Guzmán, A Belmonte, J P Torres, Tyler W Neely, Mark Baker, Reuven Gordon, Alexander B Stilgoe, Jacquiline Romero, Andrew G White, Robert Fickler, Alan E Willner, Guodong Xie, Benjamin McMorran, and Andrew M Weiner, “Roadmap on structured light,” Journal of Optics 19, 013001 (2016).
  * Fickler  _et al._ (2012) Robert Fickler, Radek Lapkiewicz, William N. Plick, Mario Krenn, Christoph Schaeff, Sven Ramelow, and Anton Zeilinger, “Quantum entanglement of high angular momenta,” Science 338, 640–643 (2012).
  * Bai  _et al._ (2022) Yihua Bai, Haoran Lv, Xin Fu, and Yuanjie Yang, “Vortex beam: generation and detection of orbital angular momentum [Invited],” Chin. Opt. Lett. 20, 012601 (2022).
  * Friese  _et al._ (1998) M. E. J. Friese, T. A. Nieminen, N. R. Heckenberg, and H. Rubinsztein-Dunlop, “Optical alignment and spinning of laser-trapped microscopic particles,” Nature 394, 348–350 (1998).
  * Padgett and Bowman (2011) Miles Padgett and Richard Bowman, “Tweezers with a twist,” Nature Photonics 5, 343–348 (2011).
  * Grier (2003) David G. Grier, “A revolution in optical manipulation,” Nature 424, 810–816 (2003).
  * Curtis and Grier (2003) Jennifer E. Curtis and David G. Grier, “Structure of optical vortices,” Physical Review Letters 90, 133901 (2003).
  * Lavery  _et al._ (2013) Martin P. J. Lavery, Fiona C. Speirits, Stephen M. Barnett, and Miles J. Padgett, “Detection of a spinning object using light’s orbital angular momentum,” Science 341, 537–540 (2013).
  * Mair  _et al._ (2001) Alois Mair, Alipasha Vaziri, Gregor Weihs, and Anton Zeilinger, “Entanglement of the orbital angular momentum states of photons,” Nature 412, 313–316 (2001).
  * Molina-Terriza  _et al._ (2007) Gabriel Molina-Terriza, Juan P. Torres, and Lluis Torner, “Twisted photons,” Nature Physics 3, 305–310 (2007).
  * Magaña-Loaiza and Boyd (2019) Omar S Magaña-Loaiza and Robert W Boyd, “Quantum imaging and information,” Reports on Progress in Physics 82, 124401 (2019).
  * Ding  _et al._ (2015) Dong-Sheng Ding, Wei Zhang, Zhi-Yuan Zhou, Shuai Shi, Guo-Yong Xiang, Xi-Shi Wang, Yun-Kun Jiang, Bao-Sen Shi, and Guang-Can Guo, “Quantum storage of orbital angular momentum entanglement in an atomic ensemble,” Physical Review Letters 114, 050502 (2015).
  * Langford  _et al._ (2004) N. K. Langford, R. B. Dalton, M. D. Harvey, J. L. O’Brien, G. J. Pryde, A. Gilchrist, S. D. Bartlett, and A. G. White, “Measuring entangled qutrits and their use for quantum bit commitment,” Physical Review Letters 93, 053601 (2004).
  * Zhang  _et al._ (2007) Pei Zhang, Xi-Feng Ren, Xu-Bo Zou, Bi-Heng Liu, Yun-Feng Huang, and Guang-Can Guo, “Demonstration of one-dimensional quantum random walks using orbital angular momentum of photons,” Physical Review A 75, 052310 (2007).
  * Nagali  _et al._ (2009) Eleonora Nagali, Fabio Sciarrino, Francesco De Martini, Lorenzo Marrucci, Bruno Piccirillo, Ebrahim Karimi, and Enrico Santamato, “Quantum information transfer from spin to orbital angular momentum of photons,” Physical Review Letters 103, 013601 (2009).
  * Zhang  _et al._ (2012) Pei Zhang, Yan Jiang, Rui-Feng Liu, Hong Gao, Hong-Rong Li, and Fu-Li Li, “Implementing the deutsch's algorithm with spin-orbital angular momentum of photon without interferometer,” Optics Communications 285, 838–841 (2012).
  * Bozinovic  _et al._ (2013) Nenad Bozinovic, Yang Yue, Yongxiong Ren, Moshe Tur, Poul Kristensen, Hao Huang, Alan E. Willner, and Siddharth Ramachandran, “Terabit-scale orbital angular momentum mode division multiplexing in fibers,” Science 340, 1545–1548 (2013).
  * Wang  _et al._ (2012) Jian Wang, Jeng-Yuan Yang, Irfan M. Fazal, Nisar Ahmed, Yan Yan, Hao Huang, Yongxiong Ren, Yang Yue, Samuel Dolinar, Moshe Tur, and Alan E. Willner, “Terabit free-space data transmission employing orbital angular momentum multiplexing,” Nature Photonics 6, 488–496 (2012).
  * Ouyang  _et al._ (2021) Xu Ouyang, Yi Xu, Mincong Xian, Ziwei Feng, Linwei Zhu, Yaoyu Cao, Sheng Lan, Bai-Ou Guan, Cheng-Wei Qiu, Min Gu, and Xiangping Li, “Synthetic helical dichroism for six-dimensional optical orbital angular momentum multiplexing,” Nature Photonics 15, 901–907 (2021).
  * Zhang  _et al._ (2021) Mingsi Zhang, Haoran Ren, Xu Ouyang, Meiling Jiang, Yudong Lu, Yanwen Hu, Shenhe Fu, Zhen Li, Zhenqiang Chen, Bai-Ou Guan, Yaoyu Cao, and Xiangping Li, “Nanointerferometric discrimination of the spin–orbit hall effect,” ACS Photonics 8, 1169–1174 (2021).
  * Wang  _et al._ (2018) Le Wang, Xincheng Jiang, Li Zou, and Shengmei Zhao, “Two-dimensional multiplexing scheme both with ring radius and topological charge of perfect optical vortex beam,” Journal of Modern Optics 66, 87–92 (2018).
  * Molina-Terriza  _et al._ (2001) Gabriel Molina-Terriza, Juan P. Torres, and Lluis Torner, “Management of the angular momentum of light: Preparation of photons in multidimensional vector states of angular momentum,” Physical Review Letters 88, 013601 (2001).
  * Rodenburg  _et al._ (2012) Brandon Rodenburg, Martin P. J. Lavery, Mehul Malik, Malcolm N. O’Sullivan, Mohammad Mirhosseini, David J. Robertson, Miles Padgett, and Robert W. Boyd, “Influence of atmospheric turbulence on states of light carrying orbital angular momentum,” Optics Letters 37, 3735 (2012).
  * Rodenburg  _et al._ (2014) Brandon Rodenburg, Mohammad Mirhosseini, Mehul Malik, Omar S Magaña-Loaiza, Michael Yanakas, Laura Maher, Nicholas K Steinhoff, Glenn A Tyler, and Robert W Boyd, “Simulating thick atmospheric turbulence in the lab with application to orbital angular momentum communication,” New Journal of Physics 16, 033020 (2014).
  * Paterson (2005) C. Paterson, “Atmospheric turbulence and orbital angular momentum of single photons for optical communication,” Physical Review Letters 94, 153901 (2005).
  * Tyler and Boyd (2009) Glenn A. Tyler and Robert W. Boyd, “Influence of atmospheric turbulence on the propagation of quantum states of light carrying orbital angular momentum,” Optics Letters 34, 142 (2009).
  * Ndagano  _et al._ (2017) Bienvenu Ndagano, Nokwazi Mphuthi, Giovanni Milione, and Andrew Forbes, “Comparing mode-crosstalk and mode-dependent loss of laterally displaced orbital angular momentum and hermite–gaussian modes for free-space optical communication,” Optics Letters 42, 4175 (2017).
  * Krenn  _et al._ (2014) Mario Krenn, Robert Fickler, Matthias Fink, Johannes Handsteiner, Mehul Malik, Thomas Scheidl, Rupert Ursin, and Anton Zeilinger, “Communication with spatially modulated light through turbulent air across vienna,” New Journal of Physics 16, 113028 (2014).
  * Cox  _et al._ (2019) Mitchell A. Cox, Luthando Maqondo, Ravin Kara, Giovanni Milione, Ling Cheng, and Andrew Forbes, “The resilience of hermite– and laguerre–gaussian modes in turbulence,” J. Lightwave Technol. 37, 3911–3917 (2019).
  * Courtial  _et al._ (1998) J. Courtial, K. Dholakia, D. A. Robertson, L. Allen, and M. J. Padgett, “Measurement of the rotational frequency shift imparted to a rotating light beam possessing orbital angular momentum,” Physical Review Letters 80, 3217–3219 (1998).
  * Zhang  _et al._ (2014) Wuhong Zhang, Qianqian Qi, Jie Zhou, and Lixiang Chen, “Mimicking faraday rotation to sort the orbital angular momentum of light,” Physical Review Letters 112, 153601 (2014).
  * Padgett  _et al._ (1996) M. Padgett, J. Arlt, N. Simpson, and L. Allen, “An experiment to observe the intensity and phase structure of laguerre–gaussian laser modes,” American Journal of Physics 64, 77–82 (1996).
  * Yongxin  _et al._ (2011) Liu Yongxin, Tao Hua, Pu Jixiong, and L Baida, “Detecting the topological charge of vortex beams using an annular triangle aperture,” Optics & Laser Technology 43, 1233–1236 (2011).
  * Karimi  _et al._ (2009) Ebrahim Karimi, Bruno Piccirillo, Eleonora Nagali, Lorenzo Marrucci, and Enrico Santamato, “Efficient generation and sorting of orbital angular momentum eigenmodes of light by thermally tuned q-plates,” Applied Physics Letters 94, 231124 (2009).
  * Shin  _et al._ (2016) Hoo-Chang Shin, Holger R. Roth, Mingchen Gao, Le Lu, Ziyue Xu, Isabella Nogues, Jianhua Yao, Daniel Mollura, and Ronald M. Summers, “Deep convolutional neural networks for computer-aided detection: CNN architectures, dataset characteristics and transfer learning,” IEEE Transactions on Medical Imaging 35, 1285–1298 (2016).
  * Quattoni  _et al._ (2008) Ariadna Quattoni, Michael Collins, and Trevor Darrell, “Transfer learning for image classification with sparse prototype representations,” in  _2008 IEEE Conference on Computer Vision and Pattern Recognition_ (IEEE, 2008).
  * Taigman  _et al._ (2014) Yaniv Taigman, Ming Yang, Marc’Aurelio Ranzato, and Lior Wolf, “Deepface: Closing the gap to human-level performance in face verification,” Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition (CVPR) (2014).
  * Melnikov  _et al._ (2018) Alexey A. Melnikov, Hendrik Poulsen Nautrup, Mario Krenn, Vedran Dunjko, Markus Tiersch, Anton Zeilinger, and Hans J. Briegel, “Active learning machine learns to create new quantum experiments,” Proceedings of the National Academy of Sciences 115, 1221–1226 (2018).
  * Doster and Watnik (2017) Timothy Doster and Abbie T. Watnik, “Machine learning approach to OAM beam demultiplexing via convolutional neural networks,” Applied Optics 56, 3386 (2017).
  * You  _et al._ (2020) Chenglong You, Mario A. Quiroz-Juárez, Aidan Lambert, Narayan Bhusal, Chao Dong, Armando Perez-Leija, Amir Javaid, Roberto de J. León-Montiel, and Omar S. Magaña-Loaiza, “Identification of light sources using machine learning,” Applied Physics Reviews 7, 021404 (2020).
  * Bhusal  _et al._ (2021a) Narayan Bhusal, Sanjaya Lohani, Chenglong You, Mingyuan Hong, Joshua Fabre, Pengcheng Zhao, Erin M. Knutson, Ryan T. Glasser, and Omar S. Magaña-Loaiza, “Spatial mode correction of single photons using machine learning,” Advanced Quantum Technologies 4, 2000103 (2021a).
  * Gibson  _et al._ (2004) Graham Gibson, Johannes Courtial, Miles J. Padgett, Mikhail Vasnetsov, Valeriy Pas'ko, Stephen M. Barnett, and Sonja Franke-Arnold, “Free-space information transfer using light beams carrying orbital angular momentum,” Optics Express 12, 5448 (2004).
  * Liao  _et al._ (2021) Kun Liao, Ye Chen, Zhongcheng Yu, Xiaoyong Hu, Xingyuan Wang, Cuicui Lu, Hongtao Lin, Qingyang Du, Juejun Hu, and Qihuang Gong, “All-optical computing based on convolutional neural networks,” Opto-Electronic Advances 4, 200060–200060 (2021).
  * Bhusal  _et al._ (2021b) Narayan Bhusal, Mingyuan Hong, Nathaniel R. Miller, Mario A Quiroz-Juarez, Roberto de J. Leon-Montiel, Chenglong You, and Omar S. Magana-Loaiza, “Smart quantum statistical imaging beyond the abbe-rayleigh criterion,” arXiv:2110.05446 (2021b), arXiv:2110.05446 [quant-ph] .
  * Zhang and Zhao (2019) Chao Zhang and Yufei Zhao, “Orbital angular momentum nondegenerate index mapping for long distance transmission,” IEEE Transactions on Wireless Communications 18, 5027–5036 (2019).
  * Krenn  _et al._ (2016) Mario Krenn, Johannes Handsteiner, Matthias Fink, Robert Fickler, Rupert Ursin, Mehul Malik, and Anton Zeilinger, “Twisted light transmission over 143 km,” Proceedings of the National Academy of Sciences 113, 13648–13653 (2016).
  * Zhao  _et al._ (2011) S. M. Zhao, J. Leach, L. Y. Gong, J. Ding, and B. Y. Zheng, “Aberration corrections for free-space optical communications in atmosphere turbulence using orbital angular momentum states,” Optics Express 20, 452 (2011).
  * Zhou  _et al._ (2016) Zhi-Yuan Zhou, Zhi-Han Zhu, Shi-Long Liu, Shi-Kai Liu, Kai Wang, Shuai Shi, Wei Zhang, Dong-Sheng Ding, and Bao-Sen Shi, “Generation and reverse transformation of twisted light by spatial light modulator,” arXiv:1612.04482 (2016), arXiv:1612.04482 [physics.optics] .
  * Collins (1970) Stuart A. Collins, “Lens-system diffraction integral written in terms of matrix optics$\ast$,” J. Opt. Soc. Am. 60, 1168–1177 (1970).
  * Bos  _et al._ (2015) Jeremy P. Bos, Michael C. Roggemann, and V. S. Rao Gudimetla, “Anisotropic non-kolmogorov turbulence phase screens with variable orientation,” Appl. Opt. 54, 2039–2045 (2015).
  * Glindemann  _et al._ (1993) A. Glindemann, R.G. Lane, and J.C. Dainty, “Simulation of time-evolving speckle patterns using kolmogorov statistics,” Journal of Modern Optics 40, 2381–2388 (1993).
  * He  _et al._ (2016) Kaiming He, Xiangyu Zhang, Shaoqing Ren, and Jian Sun, “Deep residual learning for image recognition,” Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition (CVPR) (2016).
  * Kingma and Ba (2014) Diederik P. Kingma and Jimmy Ba, “Adam: A method for stochastic optimization,” arXiv:1412.6980 (2014), arXiv:1412.6980 [cs.LG] .
  * Fernando  _et al._ (2014) Basura Fernando, Amaury Habrard, Marc Sebban, and Tinne Tuytelaars, “Subspace alignment for domain adaptation,” arXiv:1409.5241 (2014), arXiv:1409.5241 [cs.CV] .
  * Patel  _et al._ (2015) Vishal M Patel, Raghuraman Gopalan, Ruonan Li, and Rama Chellappa, “Visual domain adaptation: A survey of recent advances,” IEEE Signal Processing Magazine 32, 53–69 (2015).
  * Krizhevsky  _et al._ (2017) Alex Krizhevsky, Ilya Sutskever, and Geoffrey E. Hinton, “ImageNet classification with deep convolutional neural networks,” Communications of the ACM 60, 84–90 (2017).
  * Liu  _et al._ (2019) Zhanwei Liu, Shuo Yan, Haigang Liu, and Xianfeng Chen, “Superhigh-resolution recognition of optical vortex modes assisted by a deep-learning method,” Physical Review Letters 123, 183902 (2019).
  * Zhang and Sabuncu (2018) Zhilu Zhang and Mert R. Sabuncu, “Generalized cross entropy loss for training deep neural networks with noisy labels,” in  _Proceedings of the 32nd International Conference on Neural Information Processing Systems_ , NIPS’18 (Curran Associates Inc., Red Hook, NY, USA, 2018) p. 8792 C8802.