# Fast, high precision dynamics in quantum optimal control theory

Mogens Dalgaard  Department of Physics and Astronomy, Aarhus University, Ny
Munkegade 120, 8000 Arhus C, Denmark  Felix Motzoi  Forschungszentrum Jülich,
Institute of Quantum Control (PGI-8), D-52425 Jülich, Germany

(July 2019; July 2019)

###### Abstract

Quantum optimal control theory is becoming increasingly crucial as quantum
devices become more precise, but the need to quickly optimize these systems
classically remains a significant bottleneck in their operation. Here we
present a new theoretical quantum control framework for much faster
optimization than the state of the art by replacing standard time propagation
with a product of short-time propagators, each calculated using the Magnus
expansion. The derived formulas for exact series terms and their gradients,
based on earlier approximate integrals in a simulation setting, allow us to
subsume the high cost of calculating commutators and integrals as an initial
overhead. This provides an order of magnitude speedup for quantum control
optimization.

##  I Introduction

Quantum technologies may solve societally relevant problems within areas such
as optimization Farhi  _et al._ (2014), equation solving Harrow  _et al._
(2009), drug design Peruzzo  _et al._ (2014); McClean  _et al._ (2016), and
machine learning Schuld  _et al._ (2015); Biamonte  _et al._ (2017) in
combinatorially difficult regimes where classical computers are expected to
fail Dalzell  _et al._ (2020). In recent years, we have witnessed the small-
scale implementation of some of these ideas Boixo  _et al._ (2016); Kandala
_et al._ (2017); Kokail  _et al._ (2019); Omran  _et al._ (2019); Harrigan
_et al._ (2021). However, large-scale implementations still require dramatic
improvement in our ability to accurately simulate and control subparts of
these systems.

Quantum control optimization is formally treated within quantum control theory
d’Alessandro (2007); Werschnik and Gross (2007); Glaser  _et al._ (2015);
Goerz (2015), where recent progress includes: frequency domain optimization
Dong  _et al._ (2020); Sørensen  _et al._ (2020); Müller  _et al._ (2021),
Hessian-based optimization Goodwin and Kuprov (2016); Dalgaard  _et al._
(2020a); Jensen  _et al._ (2021), optimization of many-body matrix product
states Jensen  _et al._ (2020); Omran  _et al._ (2019), noise-resilient
control Müller  _et al._ (2018); Liu  _et al._ (2019); Khani  _et al._ (2012);
Gupta  _et al._ (2020); Dalgaard  _et al._ (2021a), reinforcement learning-
based optimization Palittapongarnpim  _et al._ (2017); Bukov  _et al._ (2018);
Xu  _et al._ (2019); An and Zhou (2019); Dalgaard  _et al._ (2020b); Yao  _et
al._ (2020); Baum  _et al._ (2021), circuit optimization Mitarai  _et al._
(2018); Motzoi  _et al._ (2017); Arrazola  _et al._ (2019), feedback control
Chen  _et al._ (2020); Magrini  _et al._ (2021); Motzoi  _et al._ (2016);
Basilewitsch  _et al._ (2019), and global cost functional landscape
optimization Kokail  _et al._ (2019); Koczor and Benjamin (2020); Dalgaard
_et al._ (2021b).

In the context of quantum simulations, a widely used approach is the rich
class of numerical differential equation solvers such as Runge-Kutta Butcher
(1996). However, state-of-the-art quantum control optimization relies on
gradient-based minimization, which often outperforms gradient-free
alternatives Sørensen  _et al._ (2018). Although calculating the gradient
with, e.g., Runge-Kutta is possible Machnes  _et al._ (2018), it is usually a
numerically expensive and slow approach. For this reason, we explore in this
paper gradient-based optimization that relies on the Magnus expansion Magnus
(1954); Blanes  _et al._ (2009); Auer  _et al._ (2018); Ma  _et al._ (2020),
which benefits from significantly better convergence properties than, e.g.,
Taylor or Dyson series. However, like these other series, going beyond a few
orders is computationally (or analytically) taxing, while lower-order results
only provide approximate solutions. Therefore a standard approach is to
resolve the total simulation duration into shorter, high-accuracy time
intervals, which naturally leads to a trade-off between accuracy per time step
and the total number of time steps.

Within quantum control optimization, a current state-of-the-art gradient-based
algorithm, GRAPE Khaneja  _et al._ (2005); De Fouquieres  _et al._ (2011), is
a low-accuracy (per time step) method, viewable as an approximation of the
first term in the Magnus expansion. In this work, we demonstrate that one may
easily improve GRAPE by including higher-order terms in the Magnus expansion.
Our contribution is thus two-fold: (i) we derive the control gradient for
three different integrations schemes based on the Magnus expansion. (ii) We
compare the performance of these four (i.e., including standard GRAPE) control
optimization methods in equal computational resource comparisons on both slow
and fast dynamics, relative to the control fields’ modulation rate.

A crucial insight of our work is the ability to pre-calculate commutators and
integrals, allowing for the fast inclusion of higher-order terms in the Magnus
expansion. Our method also leads to the option of whether to approximate or
exact calculate the arising integrals. The latter approach increases accuracy
per time step but at the cost of overhead in either human labor or
computational resources. This overhead is seldom justified for standard
“single-shot” simulations, making Magnus-based integration rarely more
efficient than other frequently used integration techniques. However, quantum
control optimization generally relies on many time propagations, making this
overhead overall insignificant.

As a testbed for our methods, we optimize the global control of an Ising-type
spin chain, recently studied also in Refs. Dalgaard  _et al._ (2021b, a).
Locally interacting spin systems are compelling to study because of their
importance within solid-state physics Plischke and Bergersen (1994). In
addition, many quantum computational architectures depend on adjacently
coupled qubits Cirac and Zoller (1995); Bloch (2008); Blais  _et al._ (2020);
Wu  _et al._ (2020); Kinos  _et al._ (2021), while also having potential
applications in surface error correction codes Fowler  _et al._ (2012);
Barends  _et al._ (2014) and transport of quantum information Christandl  _et
al._ (2004); Yung (2006). Our results demonstrate up to a factor of ten in
differences in wall time consumption between the best- and worst-performing
methods in equal-accuracy comparisons on both fast and slow dynamics.

The paper is organized as follows: In Section II we demonstrate how to use the
Magnus expansion to derive different integration schemes and how to utilize
these for quantum control optimization. In Section III, we benchmark these
against each other with respect to both speed and accuracy. In Section IV, we
present concluding remarks.

##  II Theory

Figure 1: Optimization results with the rotating wave approximation. (a) The
wall time consumption for calculating the infidelity versus the numerical
error of the infidelity averaged over 100 optimized pulses for the different
methods presented in the text. We obtained the data by scanning over different
values of $\Delta t$. From this data, we find the $\Delta t$ that leads to an
average error of $10^{-6}$. In (b), we compare the optimization from the found
discretization of the two approximation methods. We repeat in (c) the
comparison for the two exact methods. Finally, in (d), we depict the smallest
infidelity found at each wall time duration for each of the different methods.

Quantum control optimization allows for shaping externally applied control
fields in order to realize dynamical operations such as state-to-state
transfers, unitary gate synthesis, or preparation of a density matrix in both
closed and open quantum systems Machnes  _et al._ (2011). One of the most
successful approaches in this regard, is the gradient pulse engineering
algorithm (GRAPE) Khaneja  _et al._ (2005); De Fouquieres  _et al._ (2011) and
its many extensions Motzoi  _et al._ (2011); Schulte-Herbrüggen  _et al._
(2011); Goodwin and Kuprov (2016); Lucarelli (2018); Sørensen  _et al._
(2018); Dalgaard  _et al._ (2020a). The general idea is to reformulate the
control problem as a minimization task of a cost-function
$\mathcal{F}[u^{(k)}]$ with respect to a set of control fields
$\\{u^{(k)}\\}_{k}$ that steers the system via a control Hamiltonian

$\displaystyle H(t)=H_{0}+\sum_{k}u^{(k)}(t)H_{k}.$ (1)

Here $H_{0}$ and $H_{k}$ denote drift and control Hamiltonians, respectively,
which we assume are time-independent. In this work, we specifically seek to
optimize pulses expressed in a set of basis functions
$u^{k}(t)=\sum_{n}b_{n}^{(k)}\phi_{n}^{(k)}(t)$. This type of parametrization
has been popularized in literature through the chopped random basis approach
Caneva  _et al._ (2011); Müller  _et al._ (2021); van Frank  _et al._ (2016);
Motzoi  _et al._ (2011); Sørensen  _et al._ (2018, 2020), which is suitable to
obtain pulses that are easy to implement experimentally with vanishing effects
Sørensen  _et al._ (2020) of filtering. With this parametrization, we now seek
to minimize the cost-function $\mathcal{F}(\\{b_{n}^{(k)}\\})$ over the pulse
parameters $\\{b_{n}^{(k)}\\}$. For gradient-based optimization we also need
the derivative $\partial_{b_{n}^{(k)}}\mathcal{F}$, which for a state-to-state
transfer depends on the end state and $\ket{\psi(T)}$ its derivative
$\partial_{b_{n}^{(k)}}\ket{\psi(T)}$ Khaneja  _et al._ (2005), where $T$
denotes the total control duration. By using the Schrödinger equation and its
derivative we obtain the coupled equations $(\hbar=1)$ Goerz (2015); Machnes
_et al._ (2018)

$\displaystyle\partial_{t}\ket{\psi(t)}$ $\displaystyle=-iH(t)\ket{\psi(t)}$
(2) $\displaystyle\partial_{t}\partial_{b_{n}^{(k)}}\ket{\psi(t)}$
$\displaystyle=-i\phi_{n}^{(k)}(t)H_{k}\ket{\psi(t)}-iH(t)\partial_{b_{n}^{(k)}}\ket{\psi(t)},$
(3)

which we may solve with, e.g., Runge-Kutta. Unfortunately, this requires
propagating an additional state $\partial_{b_{n}^{(k)}}\ket{\psi(t)}$ per
pulse parameter $b_{n}^{(k)}$, which is much more costly than the usual
alternative Khaneja  _et al._ (2005), thereby making Runge-Kutta or similar
techniques less tractable for quantum control optimization. In addition,
Runge-Kutta does not preserve normalization, which may cause a nonphysical
loss of state norm, which is only negligible at a sufficiently high number of
time steps.

Instead we may numerically solve, e.g., the Schrödinger equation, Eq. (2), for
a duration $T$ in a series of $N$ equidistant steps, $\Delta t=T/N$, via the
time evolution operator $U(\Delta t,t)\ket{\psi(t)}=\ket{\psi(t+\Delta t)}$.
The time evolution operator for the $j$th time step $t_{j}=j\Delta t$ is
formally given by $U(\Delta t,t_{j-1})=U_{j}=e^{-i\Omega_{j}}$, where
$\Omega_{j}=\sum_{n=1}^{\infty}\Omega_{j}^{[n]}$ denotes the Magnus expansion
Magnus (1954); Blanes  _et al._ (2009), where the first term reads

$\displaystyle\Omega_{j}^{[1]}=\int_{t_{j}}^{t_{j}+\Delta
t}dt_{1}H(t_{1})=\Delta tH_{0}+\sum_{k}c^{(1)}_{k,j}H_{k}.$ (4)

Note in the above equation we have inserted Eq. (1) and implicitly carried out
the integral

$\displaystyle c^{(1)}_{k,j}=\int_{t_{j}}^{t_{j}+\Delta
t}dt_{1}u^{(k)}(t_{1})$ (5)

The second term in the Magnus expansion reads

$\displaystyle\Omega_{j}^{[2]}$
$\displaystyle=-\frac{i}{2}\int_{t_{j}}^{t_{j}+\Delta
t}dt_{1}\int_{t_{j}}^{t_{1}}dt_{2}[H(t_{2}),H(t_{1})]$
$\displaystyle=-i\sum_{k}c^{(2)}_{k,j}[H_{0},H_{k}]-i\sum_{k<k^{\prime}}c^{(3)}_{k,k^{\prime},j}[H_{k},H_{k}^{\prime}],$
(6)

where we similarly to before have collected the integrals

$\displaystyle c^{(2)}_{k,j}=$
$\displaystyle\frac{1}{2}\int_{t_{j}}^{t_{j}+\Delta
t}dt_{1}\int_{t_{j}}^{t_{1}}dt_{2}\big{(}u^{(k)}(t_{1})-u^{(k)}(t_{2})\big{)}$
(7) $\displaystyle c^{(3)}_{k,k^{\prime},j}=$
$\displaystyle\frac{1}{2}\int_{t_{j}}^{t_{j}+\Delta
t}dt_{1}\int_{t_{j}}^{t_{1}}dt_{2}\bigg{(}u^{(k)}(t_{2})u^{(k^{\prime})}(t_{1})$
$\displaystyle-u^{(k^{\prime})}(t_{2})u^{(k)}(t_{1})\bigg{)}.$ (8)

Note, that using the parametrization
$u^{k}(t)=\sum_{n}b_{n}^{(k)}\phi_{n}^{(k)}(t)$, we may precalculate the
various integrals over the basis functions $\phi_{n}^{(k)}$, which allows for
significantly faster computation. We elaborate further on this in Appendix A.
Each increasing term in the Magnus expansion grows in complexity, hence, we
must for all practical purposes truncate the Magnus expansion at some finite
term $m$. Here the per step error
$||\Omega_{j}-\sum_{n=1}^{m}\Omega_{j}^{[n]}||$ scales as $\mathcal{O}(\Delta
t^{2m+1})$ and the accumulated error in the final state $\ket{\psi(T)}$ scales
as $\mathcal{O}(\Delta t^{2m})$ Blanes  _et al._ (2009). Note that written in
the form of Eq. (4) and (6), we may also pre-calculate the commutators, which
again allows for faster computation of the evolution.

We may obtain the derivatives by using the multi-variate chain-rule similarly
to Ref. Motzoi  _et al._ (2011); Sørensen  _et al._ (2018),

$\displaystyle\frac{\partial\mathcal{F}}{\partial
b_{n}^{(k)}}=\sum_{j}\Bigg{(}\frac{\partial\mathcal{F}}{\partial
c_{k,j}^{(1)}}\frac{\partial c_{k,j}^{(1)}}{\partial
b_{n}^{(k)}}+\frac{\partial\mathcal{F}}{\partial c_{k,j}^{(2)}}\frac{\partial
c_{k,j}^{(2)}}{\partial b_{n}^{(k)}}$
$\displaystyle+\sum_{k^{\prime}<k}\frac{\partial\mathcal{F}}{\partial
c_{k^{\prime},k,j}^{(3)}}\frac{\partial c_{k^{\prime},k,j}^{(3)}}{\partial
b_{n}^{(k)}}+\sum_{k^{\prime}>k}\frac{\partial\mathcal{F}}{\partial
c_{k,k^{\prime},j}^{(3)}}\frac{\partial c_{k,k^{\prime},j}^{(3)}}{\partial
b_{n}^{(k)}}+\ldots\Bigg{)}.$ (9)

Here the derivatives $\partial\mathcal{F}/\partial c_{k,j}^{(1)}$,
$\partial\mathcal{F}/\partial c_{k,j}^{(2)},\ldots$ are calculable using
standard GRAPE propagation Khaneja  _et al._ (2005); De Fouquieres  _et al._
(2011), which requires two propagations: one forward and one backwards. This
is significantly smaller than the alternative approach, Eq. (2) and (3), which
requires one additional propagation of $\partial_{b_{n}^{(k)}}\ket{\psi(t)}$
per pulse parameter $b_{n}^{(k)}$ thereby making the alternative much more
expensive.

In this work, we benchmark four different integration schemes for quantum
control: two-second order methods with one explicitly calculating
$\Omega_{j}^{[1]}$ (M2exact) and one approximating it via midpoint
interpolation (M2approx). As well, we study two fourth-order methods with one
explicitly calculating $\Omega_{j}^{[1]}+\Omega_{j}^{[2]}$ (M4exact) and one
approximating it via Gauss–Legendre quadrature (M4approx). Note, that M2approx
and M4approx have previously been used as numerical integration schemes Blanes
_et al._ (2009), while M2approx may readily be seen as the current standard
(GRAPE) in the quantum control literature Khaneja  _et al._ (2005); De
Fouquieres  _et al._ (2011); Machnes  _et al._ (2011); Sørensen  _et al._
(2018); Dalgaard  _et al._ (2020a). Besides M2approx, none of these methods
have, to the best of our knowledge, been used to perform gradient-based
quantum control optimization before. As mentioned in the introduction, this
paper aims to derive the control gradient needed for state-of-the-art
optimization for these three alternatives to standard GRAPE and benchmark all
of these methods against each other for control optimization. We elaborate on
the more technical details of the former in Appendix A and provide the
comparison in the following section.

Figure 2: Optimization results without the rotating wave approximation. (a)
The wall time consumption for calculating the infidelity for a state transfer
on a spin chain without the rotating wave approximation plotted against the
numerical error of the infidelity averaged over 100 optimized pulses for the
different methods presented in the text. We obtained the data by scanning over
different values of $\Delta t$. From this data, we find the $\Delta t$ that
leads to an average error of $10^{-6}$. In (b), we compare the optimization
from the found discretization of the two approximation methods. We repeat in
(c) the comparison for the two exact methods. In (d), we depict the smallest
infidelity at each wall time duration for each of the different methods.

##  III Results

###  III.1 Controlling a spin chain within the RWA

To compare the different methods, we consider control of a one-dimensional
spin chain where we model an Ising-type spin-spin interaction for nearest and
next-nearest neighboring spins, which was also studied recently in the context
of optimal control theory in Refs. Dalgaard  _et al._ (2021b, a). We assume
access to two independent global control fields $u^{(x)}$ and $u^{(y)}$
orthogonal to the chain and driven at qubit-resonance. With this, the drift
Hamiltonian becomes

$\displaystyle
H_{0}=\frac{1}{2}\sum_{j}\omega_{j}\sigma_{j}^{z}-J\sum_{j}\sigma_{j}^{z}\sigma_{j+1}^{z}-g\sum_{j}\sigma_{j}^{z}\sigma_{j+2}^{z},$
(10)

where $J$ is the nearest spin-spin interaction, $g=J/10$ is the next-nearest
spin-spin interaction, and $\sigma$ denotes the Pauli spin operator. Here we
model periodic boundary conditions and assume an isotropic spin chain
($\omega=\omega_{j}$ for all $j$). The control part of the Hamiltonian is

$\displaystyle H_{c}(t)$ $\displaystyle=2\sum_{k=x,y}u^{(k)}(t)\cos(\omega
t)\sum_{j}\sigma_{j}^{k}$ (11)

such that $H(t)=H_{0}+H_{c}(t)$. We may derive an effective Hamiltonian by
using the rotating wave approximation (RWA)

$\displaystyle H_{\text{RWA}}(t)$
$\displaystyle=-J\sum_{j}\sigma_{j}^{z}\sigma_{j+1}^{z}-g\sum_{j}\sigma_{j}^{z}\sigma_{j+2}^{z}$
$\displaystyle+\sum_{k=x,y}u^{(k)}(t)\sum_{j}\sigma_{j}^{k}.$ (12)

We express the control pulses in a modulated Fourier basis
$\phi^{(k)}_{n}(t)$=$s(t)\cos(\pi nt/T)$ and $\phi^{(k)}_{n}(t)=s(t)\sin(\pi
nt/T)$ for even and odd values of $n=1,2\ldots 8$, respectively, and $k=x,y$.
Here the modulation function $s(t)$ denotes a shape function that enforces the
pulse to smoothly start and end at zero

$\displaystyle
s(t)=\begin{cases}&\frac{1}{2}[\cos\left(\pi(\frac{t}{\tau}-1)\right)+1],\text{
for }0\leq t<\tau,\\\ &1\text{ for }\tau\leq t<T-\tau,\\\
&\frac{1}{2}[\cos\left(\pi(\frac{t-T}{\tau}+1)\right)+1]\text{ for }T-\tau\leq
t\leq T.\\\ \end{cases}$ (13)

For optimization we use a sequential programming (least squares) algorithm
Kraft (1988) implemented in the Python library Scipy Virtanen  _et al._
(2020). This optimization algorithm allows us to handle non-linear constraints
of the form

$\displaystyle u_{\min}^{(k)}\leq u^{(k)}(t)\leq u_{\max}^{(k)}\text{ for all
}k,t,$ (14)

with $u_{\min}^{(k)}$ and $u_{\max}^{(k)}$ denoting the minimal and maximal
permissible control fields, respectively, which in this work we choose as
$u_{\max}^{(k)}=-u_{\min}^{(k)}=J$ for $k=x,y$.

We seek to make a state transfer between two degenerate ground states from an
initial state $\ket{\psi_{0}}=\ket{0,0,\ldots,0}$ and a target state
$\ket{\psi_{t}}=\ket{1,1,\ldots,1}$, which requires temporarily populating a
set of excited states by the control fields. This task is solvable by
minimizing the infidelity of the transfer

$\displaystyle\mathcal{F}=1-|\braket{\psi_{t}}{\psi(T)}|^{2},$ (15)

where $\ket{\psi(T)}$ denotes the solution to the Schrödinger equation from
the initial state $\ket{\psi_{0}}=\ket{\psi(t=0)}$.

In the following, we compare the various integration methods for gradient-
based control optimization. In order to facilitate such an analysis, we
initiate a three step procedure: i) We draw $100$ seeds (initial random
guesses) and optimize these using a large number of time steps, $N$, in order
to estimate the “true” infidelity $\mathcal{F}_{\text{true}}$ at $TJ=2.9$,
which for the best solutions lead to an infidelity around $10^{-4}$. ii) Using
the optimized pulses, we now scan over different values of $\Delta t$ in order
to calculate the infidelity $\mathcal{F}_{\text{calculated}}(\Delta t)$. For
each $\Delta t$ we further calculate the simulation error
$|\mathcal{F}_{\text{calculated}}(\Delta t)-\mathcal{F}_{\text{true}}|$ and
monitor the wall time (i.e., computational time). Here the wall time is
without the initialization time, which we discuss further in Sec. III.3. We
depict in Fig. 1(a) the average simulation error as a function of the average
wall time for the various methods. From the figure, we observe an excellent
improvement in using the two fourth-order methods, which provides around a
factor of $10$ in speed up for an average accuracy of $10^{-6}$.
Interestingly, M2exact performs worse than M2approx, while M4exact performs
better than M4approx. iii) At last, we also compare the various methods for
control optimization. From the obtained data, we estimate for each integration
scheme the $\Delta t$ that on average leads to a simulation error of
$10^{-6}$. This error threshold is a few orders lower than current
experimental reachable infidelities and, therefore, should suffice for a fair
comparison. We note that the higher-order methods generally perform better in
the low-error regime, i.e., closer to machine precision.

Using the found $\Delta t$, we then optimize the 100 seeds used in step i)
again, where we now save the infidelity and wall time for each incremental
update made by GRAPE. We depict the optimization trajectories (reached
infidelity versus wall time) in Fig. 1(b) and Fig. 1(c) for the two
approximation and the two exact methods, respectively. For an easy comparison,
we further depict in Fig. 1(d) the smallest infidelity at each wall time
instance for each of the different methods. As expected, the different methods
converge to the same infidelity-minimum but at different speeds. We observe
around a factor of 10 speedups with the fourth-order methods compared to the
second-order methods.

###  III.2 Controlling a spin chain without the RWA

Figure 3: The error of time propagation at various values of $\Delta t$ for
(a) with and (b) without the RWA. Here we observe that each method is
consistent with its analytically predicted error scaling. We also depict the
initialization times for (c) with and (d) without the RWA. For the two exact
methods, we pre-calculate the integrals (4)-(8) numerically and note that,
where possible, calculating these analytically further reduces this
initialization overhead.

We now make an additional comparison, where we seek to simulate and optimize
the spin chain without the rotating wave approximation This is easily
achievable by including the cosine drives from Eq. (11) into our basis
functions $\phi_{n}^{(k)}(t)\rightarrow\phi_{n}^{(k)}(t)\cos(\omega t)$. Not
having the rotating wave approximation requires a larger number of time steps
$N$ to resolve the fast oscillations from the external fields, making
simulations typically slow. In this section, we model a resonance frequency of
$\omega=20J$, which is in a regime Motzoi  _et al._ (2011) where the rotating
wave approximation typically leads to errors compatible with the previously
optimized infidelities. In other words, this is a regime where we need the
more exact model in realistic simulations.

We now repeat the procedure from the previous section in this new scenario. In
Fig. 2(a), we depict the simulation error versus the wall time for each of the
four methods averaged over 100 optimized pulses. Overall we see the same
tendencies as before: we obtain around a factor of 10 speed up with the
fourth-order methods relative to the second-order ones. We still observe the
counter-intuitive tendency that M2exact performs worse than M2approx. In Fig.
2(b) and Fig. 2(c) we depict the optimization trajectories for the approximate
and the exact methods, respectively. Similarly to before, we also depict in
Fig. 2(d) the smallest infidelity at each wall time instance. We again observe
a factor of 10 in speed improvement with the two fourth-order methods for
optimization.

###  III.3 Further comparison and discussion

Before concluding, we briefly elaborate on the differences between the various
methods introduced and benchmarked in the preceding sections. First, we depict
the aforementioned $\Delta t$-scans in Fig. 3(a) and (b) with and without the
RWA, respectively. Here we observe that each method is consistent with the
second and fourth-order scaling also marked in the figures.

We also depict the initialization time for each method in Fig. 3(c) and (d)
with and without the RWA, respectively. The initialization time was not
included in the wall time consumption depicted in the previous two plots, but
included here for a more fair comparison. Note that we pre-evaluate the
sampling points for the two approximation methods, while for the two exact
methods, we pre-calculate the integrals. In addition, we also pre-calculate
the commutators for the two fourth-order methods. The most expensive method is
by far M4exact. However, we believe the initialization time is heavily
reducible for a large variety of control settings by pre-calculating the
integrals analytically but at the price of overhead in human labor.
Nonetheless, since M4exact is the consistently fastest method, there are many
scenarios where this additional overhead is far smaller than the optimization
time, including where many noise trajectories must be sampled, where search
complexity requires many different initial seeds, or where Hilbert space
complexity makes matrix exponentiation far more expensive than the integrals
themselves. For instance, control optimization with random sampling (as used
in this work to sufficiently explore the control space) may easily take
between 24 and 48 hours on a standard computer. Here, a subsequent ten-fold
speed up may easily justify an additional initialization cost of one or two
hours. Otherwise, M4approx avoids this high initialization cost at the price
of a relatively small slow-down.

In this work, we have only considered errors stemming from the numerical
integration, which is justified by relatively low numerical errors associated
with standard matrix exponentiation techniques. However, there also exists a
rich class of approximations methods to the matrix exponential De Fouquieres
_et al._ (2011), such as a second-order $\mathcal{O}(\Delta t^{2})$ Suzuki-
Trotter expansion Jensen  _et al._ (2021), which may also reduce computation
time with pre-calculations (e.g., eigen-decomposition of Hamiltonian terms).
Thus, where Suzuki-Trotter methods are most applicable (large exponentiation
cost), it may be viable to instead choose much smaller $\Delta t$ to enable
their use, while having to make do with $\mathcal{O}(\Delta t^{2})$ error
scaling and associated slowdown.

##  IV Conclusion

In this work, we have derived a new quantum optimal control framework based on
improved time propagation using a discretized Magnus expansion. We obtained
several high-performing quantum control propagators within this framework,
which we tested in equal-accuracy comparisons to state-of-the-art approaches.
In doing so, we demonstrated around an order of magnitude in speedup for both
fast and slow system dynamics relative to the timescale of the control
modulation. In particular, we derived an exact fourth-order Magnus propagator
(M4exact) that is consistently faster than the previous methods for individual
runs, though at an upfront overhead cost. At the same time, the Gauss-Legendre
quadrature approximation (M4approx) offers a slight reduction in speedup while
dramatically reducing this overhead. The overhead in precomputing commutators
and integrals typically represents a small fraction of the entire computation
time for standard optimal control problems. Still, the appropriate method can
always be chosen accordingly. We foresee this new approach to have widespread
applicability in optimal control-theoretic settings, and wherever many
dynamical trajectories must be calculated, such as for Monte Carlo sampling.

One current limitation of our framework is the difficulty in generalizing the
Magnus expansion beyond the first couple of terms. One possible solution would
be to combine our framework with the work of Arnal et al. Arnal  _et al._
(2018), which expresses the Magnus expansion in right-nested integrals and
commutators. Whether or not such generalization is adaptable with our
framework remains an open question and could potentially be the focus of
future work.

###### Acknowledgements.

We thank Francesco Preti and Michael Schilling for helpful discussions. This
work was funded by the Carlsberg Foundation, the Deutsche
Forschungsgemeinschaft (DFG, German Research Foundation) under Germany’s
Excellence Strategy – Cluster of Excellence Matter and Light for Quantum
Computing (ML4Q) EXC 2004/1 – 390534769, and through the European Union’s
Horizon 2020 research and innovation programme under Grant Agreements No.
817482 (PASQuanS) and No. 820394 (ASTERIQS). The numerical results presented
in this work were obtained at the Centre for Scientific Computing, Aarhus,
phys.au.dk/forskning/cscaa.

##  Appendix A Different methods

###  A.1 M2exact

The first propagation scheme we consider is calculating the exact first term
in the Magnus expansion. Here (and in the following methods), we start by
shifting the time scale $[t_{j},t_{j}+\Delta t]$ to $[0,\Delta t]$. In this
case, we may evaluate Eq. (5) by using the parametrization
$u^{k}(t)=\sum_{n}b_{n}^{(k)}\phi_{n}^{(k)}(t)$

$\displaystyle c_{k,j}^{(1)}=\sum_{n}b_{n}^{(k)}\int_{0}^{\Delta
t}\phi_{n}^{(k)}(t_{j}+t_{1})dt_{1},$ (16)

with derivative

$\displaystyle\frac{\partial c_{k,j}^{(1)}}{b_{n}^{(k)}}=\int_{0}^{\Delta
t}\phi_{n}^{(k)}(t_{j}+t_{1})dt_{1}.$ (17)

This integrals may be calculated in advance either analytically or numerically
to arbitrary precision. In this work we calculated the integrals numerically.

###  A.2 M2approx (standard GRAPE)

The second scheme we consider is the current standard in the literature:
GRAPE, which is done by using midpoint interpolation of Eq. (16) Blanes  _et
al._ (2009), which gives the coefficient

$\displaystyle c_{k,j}^{(1)}=\Delta
t\sum_{n}b_{n}^{(k)}\phi_{n}^{(k)}(t_{j}+\Delta t/2),$ (18)

with derivative given as in Motzoi  _et al._ (2011) by

$\displaystyle\frac{\partial c_{k,j}^{(1)}}{b_{n}^{(k)}}=\Delta
t\phi_{n}^{(k)}(t_{j}+\Delta t/2).$ (19)

###  A.3 M4exact

The first term in the Magnus expansion has already been calculated so we just
need to evaluate Eq. (7) and (8). Then,

$\displaystyle c_{k,j}^{(2)}=\frac{1}{2}\sum_{n}b_{n}^{(k)}\int_{0}^{\Delta
t}\Big{(}\int_{0}^{t_{1}}\phi_{n}^{(k)}(t_{j}+t_{2})dt_{2}-\phi_{n}^{(k)}(t_{j}+t_{1})t_{1}\Big{)}dt_{1},$
(20)

with derivative

$\displaystyle\frac{\partial c_{k,j}^{(2)}}{\partial
b_{n}^{(k)}}=\frac{1}{2}\int_{0}^{\Delta
t}\Big{(}\int_{0}^{t_{1}}\phi_{n}^{(k)}(t_{j}+t_{2})dt_{2}-\phi_{n}^{(k)}(t_{j}+t_{1})t_{1}\Big{)}dt_{1}.$
(21)

Lastly we have

$\displaystyle c_{k,k^{\prime},j}^{(3)}=$
$\displaystyle\frac{1}{2}\sum_{n,m}\Big{[}b_{n}^{(k)}b_{m}^{(k^{\prime})}\int_{0}^{\Delta
t}\Big{(}\phi_{n}^{(k)}(t_{j}+t_{1})\int_{0}^{t_{1}}\phi_{m}^{(k^{\prime})}(t_{j}+t_{2})dt_{2}\Big{)}dt_{1}$
$\displaystyle-b_{n}^{(k^{\prime})}b_{m}^{(k)}\int_{0}^{\Delta
t}\Big{(}\phi_{m}^{(k^{\prime})}(t_{j}+t_{1})\int_{0}^{t_{1}}\phi_{n}^{(k)}(t_{j}+t_{2})dt_{2}\Big{)}dt_{1}\Big{]},$
(22)

with derivative

$\displaystyle\frac{\partial c_{k,k^{\prime},j}^{(3)}}{\partial
b_{l}^{(h)}}=\begin{cases}\frac{1}{2}\sum_{m}b_{m}^{(k^{\prime})}\int_{0}^{\Delta
t}\Big{(}\phi_{n}^{(k)}(t_{j}+t_{1})\int_{0}^{t_{1}}\phi_{m}^{(k^{\prime})}(t_{j}+t_{2})dt_{2}\Big{)}dt_{1}&\text{
if }l=n\text{ and }h=k\\\ \frac{1}{2}\sum_{n}b_{n}^{(k)}\int_{0}^{\Delta
t}\Big{(}\phi_{n}^{(k)}(t_{j}+t_{1})\int_{0}^{t_{1}}\phi_{m}^{(k^{\prime})}(t_{j}+t_{2})dt_{2}\Big{)}dt_{1}&\text{
if }l=m\text{ and }h=k^{\prime}\\\
-\frac{1}{2}\sum_{m}b_{m}^{(k)}\int_{0}^{\Delta
t}\Big{(}\phi_{m}^{(k^{\prime})}(t_{j}+t_{1})\int_{0}^{t_{1}}\phi_{n}^{(k)}(t_{j}+t_{2})dt_{2}\Big{)}dt_{1}&\text{
if }l=n\text{ and }h=k^{\prime}\\\
-\frac{1}{2}\sum_{n}b_{n}^{(k^{\prime})}\int_{0}^{\Delta
t}\Big{(}\phi_{m}^{(k^{\prime})}(t_{j}+t_{1})\int_{0}^{t_{1}}\phi_{n}^{(k)}(t_{j}+t_{2})dt_{2}\Big{)}dt_{1}&\text{
if }l=m\text{ and }h=k\end{cases}$ (23)

###  A.4 M4approx

The last method we consider is an approximation scheme of fourth order based
on approximating the integral with Gaussian-Legendre quadratures Blanes  _et
al._ (2009). We define $H_{1,2}=H(t_{j}+c_{1,2}\Delta t)$ and
$u_{1,2}^{(k)}=u^{(k)}(t_{j}+c_{1,2}\Delta t)$ where
$c_{1,2}=1/2\mp\sqrt{3}/6$. With this the Magnus expansion reads

$\displaystyle\Omega_{j}=$ $\displaystyle\frac{\Delta
t}{2}\Big{(}H_{1}+H_{2}\Big{)}-i\frac{\sqrt{3}}{12}\Delta t^{2}[H_{2},H_{1}]$
$\displaystyle=$ $\displaystyle\Delta tH_{0}+\Delta
t\sum_{k}\frac{u_{1}^{(k)}+u_{2}^{(k)}}{2}H_{k}$ (24)
$\displaystyle-i\frac{\sqrt{3}}{12}\Delta
t^{2}\Big{(}\sum_{k}\big{(}u_{1}^{(k)}-u_{2}^{(k)}\big{)}[H_{0},H_{k}]+\sum_{k<k^{\prime}}\big{(}u_{2}^{(k)}u_{1}^{(k^{\prime})}-u_{2}^{(k^{\prime})}u_{1}^{(k)}\big{)}[H_{k},H_{k^{\prime}}]\Big{)}.$
(25)

With this we have the coefficient

$\displaystyle c_{k,j}^{(1)}=\Delta t\frac{u_{1}^{(k)}+u_{2}^{(k)}}{2}=\Delta
t\sum_{n}b_{n}^{(k)}\frac{\phi_{n}^{(k)}(t_{1})+\phi_{n}^{(k)}(t_{2})}{2},$
(26)

where we have introduced the notation $t_{1,2}=t_{j}+c_{1,2}\Delta t$. The
derivative is

$\displaystyle\frac{\partial c_{k,j}^{(1)}}{b_{n}^{(k)}}=\Delta
t\frac{\phi_{n}^{(k)}(t_{1})+\phi_{n}^{(k)}(t_{2})}{2}.$ (27)

Then

$\displaystyle c_{k,j}^{(2)}=\frac{\sqrt{3}}{12}\Delta
t^{2}\Big{(}u_{1}^{(k)}-u_{2}^{(k)}\Big{)}=\frac{\sqrt{3}}{12}\Delta
t^{2}\sum_{b}b_{n}^{(k)}\Big{(}\phi_{n}^{(k)}(t_{1})-\phi_{n}^{(k)}(t_{2})\Big{)},$
(28)

with derivative

$\displaystyle\frac{\partial
c_{k,j}^{(2)}}{b_{n}^{(k)}}=\frac{\sqrt{3}}{12}\Delta
t^{2}\Big{(}\phi_{n}^{(k)}(t_{1})-\phi_{n}^{(k)}(t_{2})\Big{)}.$ (29)

And

$\displaystyle c_{k,k^{\prime},j}^{(3)}$
$\displaystyle=\frac{\sqrt{3}}{12}\Delta
t^{2}\big{(}u_{2}^{(k)}u_{1}^{(k^{\prime})}-u_{2}^{(k^{\prime})}u_{1}{(k)}\big{)}$
(30) $\displaystyle=\frac{\sqrt{3}}{12}\Delta
t^{2}\sum_{n,m}b_{n}^{(k)}b_{m}^{(k^{\prime})}\phi_{n}^{(k)}(t_{2})\phi_{m}^{(k^{\prime})}(t_{1})-b_{n}^{(k^{\prime})}b_{m}^{(k)}\phi_{n}^{(k^{\prime})}(t_{2})\phi_{m}^{(k)}(t_{1}),$
(31)

with derivative

$\displaystyle\frac{\partial c_{k,k^{\prime},j}^{(3)}}{\partial
b_{l}^{(h)}}=\begin{cases}\frac{\sqrt{3}}{12}\Delta
t^{2}\sum_{m}b_{m}^{(k^{\prime})}\phi_{n}^{(k)}(t_{2})\phi_{m}^{(k^{\prime})}(t_{1})&\text{
if }l=n\text{ and }h=k\\\ \frac{\sqrt{3}}{12}\Delta
t^{2}\sum_{n}b_{n}^{(k)}\phi_{n}^{(k)}(t_{2})\phi_{m}^{(k^{\prime})}(t_{1})&\text{
if }l=m\text{ and }h=k^{\prime}\\\ -\frac{\sqrt{3}}{12}\Delta
t^{2}\sum_{m}b_{m}^{(k)}\phi_{n}^{(k^{\prime})}(t_{2})\phi_{m}^{(k)}(t_{1})&\text{
if }l=n\text{ and }h=k^{\prime}\\\ -\frac{\sqrt{3}}{12}\Delta
t^{2}\sum_{n}b_{n}^{(k^{\prime})}\phi_{n}^{(k^{\prime})}(t_{2})\phi_{m}^{(k)}(t_{1})&\text{
if }l=m\text{ and }h=k\end{cases}$ (32)

## References

  * Farhi  _et al._ (2014) E. Farhi, J. Goldstone, and S. Gutmann, arXiv preprint arXiv:1411.4028 (2014).
  * Harrow  _et al._ (2009) A. W. Harrow, A. Hassidim, and S. Lloyd, Phys. Rev. Lett. 103, 150502 (2009).
  * Peruzzo  _et al._ (2014) A. Peruzzo, J. McClean, P. Shadbolt, M.-H. Yung, X.-Q. Zhou, P. J. Love, A. Aspuru-Guzik, and J. L. O’brien, Nature communications 5, 1 (2014).
  * McClean  _et al._ (2016) J. R. McClean, J. Romero, R. Babbush, and A. Aspuru-Guzik, New J. Phys. 18, 023023 (2016).
  * Schuld  _et al._ (2015) M. Schuld, I. Sinayskiy, and F. Petruccione, Contemporary Physics 56, 172 (2015).
  * Biamonte  _et al._ (2017) J. Biamonte, P. Wittek, N. Pancotti, P. Rebentrost, N. Wiebe, and S. Lloyd, Nature 549, 195 (2017).
  * Dalzell  _et al._ (2020) A. M. Dalzell, A. W. Harrow, D. E. Koh, and R. L. La Placa, Quantum 4, 264 (2020).
  * Boixo  _et al._ (2016) S. Boixo, V. N. Smelyanskiy, A. Shabani, S. V. Isakov, M. Dykman, V. S. Denchev, M. H. Amin, A. Y. Smirnov, M. Mohseni, and H. Neven, Nature communications 7, 1 (2016).
  * Kandala  _et al._ (2017) A. Kandala, A. Mezzacapo, K. Temme, M. Takita, M. Brink, J. M. Chow, and J. M. Gambetta, Nature 549, 242 (2017).
  * Kokail  _et al._ (2019) C. Kokail, C. Maier, R. van Bijnen, T. Brydges, M. K. Joshi, P. Jurcevic, C. A. Muschik, P. Silvi, R. Blatt, C. F. Roos, _et al._ , Nature 569, 355 (2019).
  * Omran  _et al._ (2019) A. Omran, H. Levine, A. Keesling, G. Semeghini, T. T. Wang, S. Ebadi, H. Bernien, A. S. Zibrov, H. Pichler, S. Choi, _et al._ , Science 365, 570 (2019).
  * Harrigan  _et al._ (2021) M. P. Harrigan, K. J. Sung, M. Neeley, K. J. Satzinger, F. Arute, K. Arya, J. Atalaya, J. C. Bardin, R. Barends, S. Boixo, _et al._ , Nature Physics 17, 332 (2021).
  * d’Alessandro (2007) D. d’Alessandro,  _Introduction to quantum control and dynamics_ (Chapman and Hall/CRC, 2007).
  * Werschnik and Gross (2007) J. Werschnik and E. Gross, Journal of Physics B: Atomic, Molecular and Optical Physics 40, R175 (2007).
  * Glaser  _et al._ (2015) S. J. Glaser, U. Boscain, T. Calarco, C. P. Koch, W. Köckenberger, R. Kosloff, I. Kuprov, B. Luy, S. Schirmer, T. Schulte-Herbrüggen, _et al._ , The European Physical Journal D 69, 1 (2015).
  * Goerz (2015) M. H. Goerz,  _Optimizing Robust Quantum Gates in Open Quantum Systems_ , Ph.D. thesis, Universität Kassel (2015).
  * Dong  _et al._ (2020) D. Dong, C.-C. Shu, J. Chen, X. Xing, H. Ma, Y. Guo, and H. Rabitz, IEEE Transactions on Control Systems Technology (2020).
  * Sørensen  _et al._ (2020) J. J. Sørensen, J. S. Nyemann, F. Motzoi, J. Sherson, and T. Vosegaard, The Journal of Chemical Physics 152, 054104 (2020).
  * Müller  _et al._ (2021) M. M. Müller, R. S. Said, F. Jelezko, T. Calarco, and S. Montangero, arXiv preprint arXiv:2104.07687 (2021).
  * Goodwin and Kuprov (2016) D. Goodwin and I. Kuprov, The Journal of chemical physics 144, 204107 (2016).
  * Dalgaard  _et al._ (2020a) M. Dalgaard, F. Motzoi, J. H. M. Jensen, and J. Sherson, Phys. Rev. A 102, 042612 (2020a).
  * Jensen  _et al._ (2021) J. H. M. Jensen, F. S. Møller, J. J. Sørensen, and J. F. Sherson, Phys. Rev. A 103, 062612 (2021).
  * Jensen  _et al._ (2020) J. H. M. Jensen, F. S. Møller, J. J. Sørensen, and J. F. Sherson, arXiv preprint arXiv:2008.06076 (2020).
  * Müller  _et al._ (2018) M. M. Müller, S. Gherardini, and F. Caruso, Scientific reports 8, 1 (2018).
  * Liu  _et al._ (2019) B.-J. Liu, X.-K. Song, Z.-Y. Xue, X. Wang, and M.-H. Yung, Phys. Rev. Lett. 123, 100501 (2019).
  * Khani  _et al._ (2012) B. Khani, S. T. Merkel, F. Motzoi, J. M. Gambetta, and F. K. Wilhelm, Phys. Rev. A 85, 022306 (2012).
  * Gupta  _et al._ (2020) R. S. Gupta, C. L. Edmunds, A. R. Milne, C. Hempel, and M. J. Biercuk, npj Quantum Information 6, 1 (2020).
  * Dalgaard  _et al._ (2021a) M. Dalgaard, C. A. Weidner, and F. Motzoi, arXiv preprint arXiv:2107.11388 (2021a).
  * Palittapongarnpim  _et al._ (2017) P. Palittapongarnpim, P. Wittek, E. Zahedinejad, S. Vedaie, and B. C. Sanders, Neurocomputing 268, 116 (2017).
  * Bukov  _et al._ (2018) M. Bukov, A. G. Day, D. Sels, P. Weinberg, A. Polkovnikov, and P. Mehta, Phys. Rev. X 8, 031086 (2018).
  * Xu  _et al._ (2019) H. Xu, J. Li, L. Liu, Y. Wang, H. Yuan, and X. Wang, npj Quantum Information 5, 1 (2019).
  * An and Zhou (2019) Z. An and D. Zhou, EPL (Europhysics Letters) 126, 60002 (2019).
  * Dalgaard  _et al._ (2020b) M. Dalgaard, F. Motzoi, J. J. Sørensen, and J. Sherson, npj Quantum Information 6, 1 (2020b).
  * Yao  _et al._ (2020) J. Yao, M. Bukov, and L. Lin, in  _Mathematical and Scientific Machine Learning_ (PMLR, 2020) pp. 605–634.
  * Baum  _et al._ (2021) Y. Baum, M. Amico, S. Howell, M. Hush, M. Liuzzi, P. Mundada, T. Merkh, A. R. Carvalho, and M. J. Biercuk, arXiv preprint arXiv:2105.01079 (2021).
  * Mitarai  _et al._ (2018) K. Mitarai, M. Negoro, M. Kitagawa, and K. Fujii, Phys. Rev. A 98, 032309 (2018).
  * Motzoi  _et al._ (2017) F. Motzoi, M. P. Kaicher, and F. K. Wilhelm, Phys. Rev. Lett. 119, 160503 (2017).
  * Arrazola  _et al._ (2019) J. M. Arrazola, T. R. Bromley, J. Izaac, C. R. Myers, K. Brádler, and N. Killoran, Quantum Science and Technology 4, 024004 (2019).
  * Chen  _et al._ (2020) H. Chen, H. Li, F. Motzoi, L. Martin, K. B. Whaley, and M. Sarovar, New J. Phys. 22, 113014 (2020).
  * Magrini  _et al._ (2021) L. Magrini, P. Rosenzweig, C. Bach, A. Deutschmann-Olek, S. G. Hofer, S. Hong, N. Kiesel, A. Kugi, and M. Aspelmeyer, Nature 595, 373 (2021).
  * Motzoi  _et al._ (2016) F. Motzoi, E. Halperin, X. Wang, K. B. Whaley, and S. Schirmer, Phys. Rev. A 94, 032313 (2016).
  * Basilewitsch  _et al._ (2019) D. Basilewitsch, F. Cosco, N. L. Gullo, M. Möttönen, T. Ala-Nissilä, C. P. Koch, and S. Maniscalco, New J. Phys. 21, 093054 (2019).
  * Koczor and Benjamin (2020) B. Koczor and S. C. Benjamin, arXiv preprint arXiv:2008.13774 (2020).
  * Dalgaard  _et al._ (2021b) M. Dalgaard, F. Motzoi, and J. Sherson, arXiv preprint arXiv:2107.00008 (2021b).
  * Butcher (1996) J. C. Butcher, Applied numerical mathematics 20, 247 (1996).
  * Sørensen  _et al._ (2018) J. Sørensen, M. Aranburu, T. Heinzel, and J. Sherson, Phys. Rev. A 98, 022119 (2018).
  * Machnes  _et al._ (2018) S. Machnes, E. Assémat, D. Tannor, and F. K. Wilhelm, Physical review letters 120, 150401 (2018).
  * Magnus (1954) W. Magnus, Communications on pure and applied mathematics 7, 649 (1954).
  * Blanes  _et al._ (2009) S. Blanes, F. Casas, J.-A. Oteo, and J. Ros, Physics reports 470, 151 (2009).
  * Auer  _et al._ (2018) N. Auer, L. Einkemmer, P. Kandolf, and A. Ostermann, Computer Physics Communications 228, 115 (2018).
  * Ma  _et al._ (2020) A. Ma, A. B. Magann, T.-S. Ho, and H. Rabitz, Phys. Rev. A 102, 013115 (2020).
  * Khaneja  _et al._ (2005) N. Khaneja, T. Reiss, C. Kehlet, T. Schulte-Herbrüggen, and S. J. Glaser, Journal of magnetic resonance 172, 296 (2005).
  * De Fouquieres  _et al._ (2011) P. De Fouquieres, S. Schirmer, S. Glaser, and I. Kuprov, Journal of Magnetic Resonance 212, 412 (2011).
  * Plischke and Bergersen (1994) M. Plischke and B. Bergersen,  _Equilibrium statistical physics_ (World Scientific, 1994).
  * Cirac and Zoller (1995) J. I. Cirac and P. Zoller, Phys. Rev. Lett. 74, 4091 (1995).
  * Bloch (2008) I. Bloch, Nature 453, 1016 (2008).
  * Blais  _et al._ (2020) A. Blais, S. M. Girvin, and W. D. Oliver, Nature Physics 16, 247 (2020).
  * Wu  _et al._ (2020) X. Wu, X. Liang, Y. Tian, F. Yang, C. Chen, Y.-C. Liu, M. K. Tey, and L. You, Chinese Physics B (2020).
  * Kinos  _et al._ (2021) A. Kinos, D. Hunger, R. Kolesov, K. Mølmer, H. de Riedmatten, P. Goldner, A. Tallaire, L. Morvan, P. Berger, S. Welinski, _et al._ , arXiv preprint arXiv:2103.15743 (2021).
  * Fowler  _et al._ (2012) A. G. Fowler, M. Mariantoni, J. M. Martinis, and A. N. Cleland, Phys. Rev. A 86, 032324 (2012).
  * Barends  _et al._ (2014) R. Barends, J. Kelly, A. Megrant, A. Veitia, D. Sank, E. Jeffrey, T. C. White, J. Mutus, A. G. Fowler, B. Campbell, _et al._ , Nature 508, 500 (2014).
  * Christandl  _et al._ (2004) M. Christandl, N. Datta, A. Ekert, and A. J. Landahl, Phys. Rev. Lett. 92, 187902 (2004).
  * Yung (2006) M.-H. Yung, Phys. Rev. A 74, 030303 (2006).
  * Machnes  _et al._ (2011) S. Machnes, U. Sander, S. J. Glaser, P. De Fouquières, A. Gruslys, S. Schirmer, and T. Schulte-Herbrüggen, Phys. Rev. A 84, 022305 (2011).
  * Motzoi  _et al._ (2011) F. Motzoi, J. M. Gambetta, S. Merkel, and F. Wilhelm, Phys. Rev. A 84, 022307 (2011).
  * Schulte-Herbrüggen  _et al._ (2011) T. Schulte-Herbrüggen, A. Spörl, N. Khaneja, and S. Glaser, Journal of Physics B: Atomic, Molecular and Optical Physics 44, 154013 (2011).
  * Lucarelli (2018) D. Lucarelli, Phys. Rev. A 97, 062346 (2018).
  * Caneva  _et al._ (2011) T. Caneva, T. Calarco, and S. Montangero, Phys. Rev. A 84, 022326 (2011).
  * van Frank  _et al._ (2016) S. van Frank, M. Bonneau, J. Schmiedmayer, S. Hild, C. Gross, M. Cheneau, I. Bloch, T. Pichler, A. Negretti, T. Calarco, _et al._ , Scientific reports 6, 34187 (2016).
  * Kraft (1988) D. Kraft, DFVLR Obersfaffeuhofen, Germany (1988).
  * Virtanen  _et al._ (2020) P. Virtanen, R. Gommers, T. E. Oliphant, M. Haberland, T. Reddy, D. Cournapeau, E. Burovski, P. Peterson, W. Weckesser, J. Bright, _et al._ , Nature methods 17, 261 (2020).
  * Arnal  _et al._ (2018) A. Arnal, F. Casas, and C. Chiralt, Journal of Physics Communications 2, 035024 (2018).


# Analysevariabilite

Diop Amadou  Diop Amadou1 and Barro Diakarya2
1Université de Ségou (République du Mali<EMAIL_ADDRESS>
2Université Thomas Sankara (Burkina Faso<EMAIL_ADDRESS>
# Analysis of the evolution of agro-climatic risks in a context of climate
variability in the region of Ségou in Mali.

Diop Amadou  Diop Amadou1 and Barro Diakarya2
1Université de Ségou (République du Mali<EMAIL_ADDRESS>
2Université Thomas Sankara (Burkina Faso<EMAIL_ADDRESS>
###### Abstract

In the Sahel region the population depends largely on rain-fed agriculture. In
West Africa in particular, climate models turn to be unable to capture some
basic features of present-day climate variability. This study proposes a
contribution to the analysis of the evolution of agro-climatic risks in a
context of climate variability. Some statistical tests are used on the main
variables of the rainy season to determine the trends and the variabilities
described by the data series. Thus, the paper provides a statistical modeling
of the different agro-climatic risks while the seasonal variability of agro-
climatic parameters were analized as well as their inter annual variability.
The study identifies the probability distributions of agroclimatic risks and
the characterization of the rainy season was clarified.

Keywords : climate variability, agro-climatic risks,seasonal evolution,
variability parameters, tests.

MSC 2020 : 62F05, 91G70, 86A08

##  1 Introduction

Climate change is at the center of current concerns. Indeed, these
consequences directly impact the socio-economic development of all countries.
In particular, in West Africa where ecosystems are fragile [15], climatic
changes, combined with anthropogenic pressure, lead to significant ecological
imbalances [4, 1]. For the IPCC, 2007, the continent’s vulnerability to
climate change would be much greater. Thus, an overall negative effect on the
production of millet, maize and sorghum could occur [16, 5], with critical
thresholds of food insecurity in many regions [12]. Frequent droughts, longer
dry spells, land degradation and significant loss of biodiversity are also
predicted among others.

These are the challenges facing West African agriculture that is mainly
rainfed and therefore very sensitive to climatic variability. In addition, as
the main lever of national economies, its peasant practice is still based on
production systems with performances subject to climatic hazards [2]. In
addition, this vulnerability was accentuated at the end of the 20th century by
a decrease in rainfall and an increase in the population[16]. Indeed, a change
of regime in the rainfall series is observed over the period 1968-1990 [11, 7,
3]. The droughts of the 1970s and 1980s also caused a sharp decrease in
rainfall of 20 to 40% between 1931-1960 and 1968-1990 in the Sahel against 15%
in the tropical and humid regions. A southward movement of isohyets of the
order of 150 to 200 km is also observed in the Sudano-Sahelian zone [6].
According to [14, 17], the decrease in rainfall amounts can explain a 35 to
45% drop in crop yields in the Sahel.

These trends greatly impact the variability of the characteristics of the
agricultural season. Indeed, one of the characteristics of the Sudano-Sahelian
climate is the existence of a close and linear relationship between the date
of onset of rains and the length of the cropping period [13]. The variation in
the onset and cessetion dates of the rainy season, causing a fluctuation in
the length of the rainy season resulting in delays and premature stops of the
rains.

One of the risks associated with these modifications is the difficulty of
adapting plants to the local ecological disturbances induced. Particularly for
certain species with relatively long reproduction cycles [10]. Thus, a good
characterization of the onset and cessetion dates of the season is essential
to optimize the results of a cropping season.

The main objective of this study is to analyze the temporal evolution of agro-
climatic risks in a context of climatic variability in the region of Ségou in
Mali. To this end, the paper is organized as follows. The section 2 describes
the study area. The section 3 is concerned with an overview of the main
statistical tests used in the study. Our results sing at section 4. We use
tests on variables characterizing the rainy season to identify the trends and
variabilities contained in the data series. Then, we propose a statistical
modeling of these different agro-climatic risks. The data analysis techniques
and tools retained for this study are presented below. The results on breaks
in the series, the evolution of seasonal and intra-seasonal agro-climatic
risks and their possible consequences are then discussed.

##  2 Study area

Located in central Mali, the Ségou region 13∘22′05′′nord, 5∘16′2′′ouest covers
an area of 64,947 km2 (approximately 5% of Mali). It is bounded to the south
by the Sikasso region, to the south-east by Burkina Faso, to the east by the
Mopti region to the north by Mauritania and the Timbuktu region and to the
west by the Koulikoro region. With the regions of Mopti and Koulikoro, it
forms what is commonly referred to as the Center of Mali. It is mainly located
in the Sahelian zone where it enjoys a semi-arid climate (average annual
rainfall: 513 mm). The presence of several rivers (it is crossed by the Niger
river (over 292 km) as well as the Bani river) allows irrigated crops. Indeed,
the area is home to the Office du Niger and the Office riz de Ségou, which
have large developed and irrigated areas. The Ségou region has 16 classified
forests covering an area of 78.860 ha.

Figure 1: Ségou region location in Mali (source : the author)

##  3 Methods

In this study, daily data for around thirty years (1990-2019) from the
synoptic station of Ségou Latitude: $13{{}^{\circ}}24^{\prime}N;$ Longitude:
$06{{}^{\circ}}09^{\prime}N$; Altitude of the station: $288.05m$ on the
variables of temperature, humidity, rains, insolation and wind speed were
used. The missing data were estimated and replaced by the daily averages
obtained over the same period.

Initially, a set of graphical representations was used to visualize the trends
of the various climatic parameters and the temperature and rainfall anomalies.
The Mann-Kendall test was implemented to measure the significance of the
trends observed on the different variables and the Pettitt test to identify
the different break dates in the data. A break is characterized by a change in
the probability law of the time series at a given instant [8]. Agronomic
criteria were used to determine the onset and cessetion dates of the season
(Sivakumar and PRESAO). Correlation analyzes were also carried out between the
different agro-climatic parameters.

  * •

Mann-Kendall test

The Mann-Kendall test [9] has been widely used to detect trends in
meteorological, hydrological and agro-meteorological time series. It is used
to determine with a nonparametric test whether a trend is identifiable in a
time series that possibly includes a seasonal component. The null hypothesis
$H_{0}$ of these tests is that there is no trend. The three alternative
hypotheses of negative, non-zero or positive trend can be chosen.

Mann-Kendall tests are based on the calculation of Kendall’s tau measuring the
association between two samples and itself based on the ranks within the
samples.

Consider a dataset consisting of $x$ values with sample size $N$. The Mann-
Kendall test (MK) calculation begins by estimating the $S$ statistic:

$S={\displaystyle\sum\limits_{i=1}^{N-1}}{\displaystyle\sum\limits_{j=i+1}^{N}}sgn\left(x_{j}-x_{i}\right)\text{
\ \ \ }for\text{ }j>I$ (1)

KENDALL and STUART (1967); MANN (1945) states that when $N\geq 8$, the
distribution of $S$ approaches the Gaussian form with mean $E(S)=0$ and
variance $V(S)$ given by:

$V\left(S\right)=\frac{N\left(N-1\right)\left(2N+5\right){\displaystyle\sum\limits_{m=1}^{SS}}t_{i}\left(m-1\right)\left(2m+5\right)m}{18}$
(2)

Where $SS$ is the number of tied groups and $t_{i}$ is the length of the
$SS^{th}$ group.

The statistic $S$ is then standardized (MK), and its significance is estimated
from the normal cumulative distribution function.

$MK=\left\\{\begin{array}[c]{c}\frac{S-1}{\sqrt{V\left(S\right)}}\rightarrow
S>0\\\ 0\rightarrow S=0\\\ \frac{S+1}{\sqrt{V\left(S\right)}}\rightarrow
S<0\end{array}\right.$ (3)

  * •

Pettitt test

The Pettitt test is non-parametric and is derived from the Mann-Withney test
one. The absence of a break in the series $\left(x_{i}\right)$ of size
$N$constitutes the null hypothesis.

Pettitt defines the variable $U_{t,N}$:

$U_{t,N}={\displaystyle\sum\limits_{i=1}^{t}}{\displaystyle\sum\limits_{j=t+1}^{N}}D_{i,j}$
(4)

Where: $D_{i,j}=sgn\left(x_{i}-x_{j}\right)$ with $sgn\left(Z\right)=1$ if
$Z>0,$ $0$ if $Z=0$ and $-1$ if $Z<0.$

He proposes to test the null hypothesis using the $K_{N}$ statistic defined by
the maximum in absolute value of $U_{t,N}$ for $t$ varying from $1$ to $N-1$
From rank theory, Pettitt shows that if $k$ denotes the value of $K_{N}$ taken
by the series studied, under the null hypothesis, the probability of exceeding
$k$ is given approximately by:

$\Pr ob\left(K_{N}>k\right)\approx
2\exp\left[\frac{-6k^{2}}{\left(N^{3}+N^{2}\right)}\right]$ (5)

For a given first-kind $\alpha$ risk, if the estimated probability is less
than $\alpha$, the null hypothesis is rejected. The series then includes a
localized break at the moment $\tau$ when $K_{N}$ is observed.

  * •

Sivakumar criterion

The agronomic criterion proposed by [13] was retained to determine the dates
of the start of the season. Indeed, Sivakumar considers the onset of the rainy
season in the Sahelian and Sudanese regions as the date from May 01 collecting
a water depth of at least 20 mm over 3 consecutive days, without any dry
sequences of more than 7 days in the following 30 days. The date of the end of
the season is the day when, after September 01, there is no more rain greater
than or equal to 5 mm for at least 20 successive days or two weeks.

  * •

PRESAO criterion for cessetion date of the season

A comparison of the cessetion dates of season was carried out using the
criterion retained in the PRESAO111PREvision Saisonnière en Afrique de l’Ouest
(Seasonal PREvision in West Africa) process set up in 1998 by the
ACMAD222African Centre of Meteorological Application for Development-
AGRHYMET333AGRométéorologie-HYdrologie-METéorologie-ABN444Autorité du Bassin
du Niger-ICRISAT555International Crops Research Institute for the Semi-Arid
Tropics consortium. The end of the season date for this criterion is when from
01 September, soil capable of holding 70 mm of available water is completely
taken up by a daily evapotranspiration loss of 5 mm.

For the fit of data to a normal distribution, we applied a Shapiro-Wilk test,
associated with a quantile-quantile plot.

  * •

Shapiro-Wilk test

Published in 1965 by Samuel Sanford Shapiro and Martin Wilk, the Shapiro-Wilk
test tests the null hypothesis that a sample $x_{1},x_{2},...x_{n}$ is from a
normally distributed population.

The $W$ test statistic is defined by:

$W=\frac{\left({\displaystyle\sum\limits_{i=1}^{n}}a_{i}x_{\left(i\right)}\right)^{2}}{{\displaystyle\sum\limits_{i=1}^{n}}\left(x_{i}-\overline{x}\right)^{2}}$
(6)

Where

  * •

$x_{\left(i\right)}$ denote the $i^{th}$ smallest number in the sample

  * •

$\overline{x}$ is the sample mean

$\left(a_{1},a_{2},...a_{n}\right)=\frac{m^{T}V^{-1}}{\left(m^{T}V^{-1}V^{-1}m\right)^{1/2}}$

Where $m=\left(m_{1},m_{2},...m_{n}\right)^{T}$ are the expectations of the
order statistics of iid variables according to a normal distribution and $V$
is the variance-covariance matrix of these order statistics.

If the p-value is greater than the chosen alpha level, then we accept the null
hypothesis.

  * •

The quantile-quantile plot (Q-Q plot) is a graphical tool that makes it
possible to assess the relevance of the fit of a given distribution to a
theoretical model. From the observed statistical series, we then calculate a
certain number of quantiles $x_{i}$. If the statistical series follows the
chosen theoretical distribution well, we should have the observed quantiles
$x_{i}$ equal to the quantiles $x_{i}^{\ast}$ associated with the theoretical
model.

##  4 Results

###  4.1 Seasonal variability of agro-climatic parameters

####  4.1.1 Seasonal temperature variability

Analysis of the seasonal evolution of minimum, maximum and average
temperatures (Figure 2) over the period 1990-2019 shows a ripple throughout
the year. Minimum temperatures present maximums in April-May with an average
of 27.56∘C for the month of May and minimums in the months of January and
December with an average of 17.26∘C for the month of January. These low
temperatures would certainly be due to the cool and dry Harmattan winds coming
from the Sahara and sweeping the entire West African sub-region. The months of
July, August and September show a very low dispersion of minimum temperatures
which correspond to the winter period of the area. Maximum temperatures also
show two peaks; the first observed during the months of March-April-May with
temperatures reaching an average of 40∘C and the second during the months of
October and November. The months of July, August and September are the cooler
months with temperatures in the range of 30-36∘C. The period of March-April-
May corresponds to the dry season, which explains the high temperatures
observed. The average temperatures vary between 25 and 28∘C for the minimums
and 30 and 34∘C for the maximums.

Figure 2: Temperature regime

####  4.1.2 Seasonal rainfall and evapotranspiration variability

Rainfall regime in the area is uni modal (Figure 3). The monthly average
rainfall shows a high variability over the year. The period from June to
September is the wettest with quantities on average over 100 mm. August is the
wettest month with average heights reaching 200 mm. However, the monthly total
shows a strong dispersion during this month. Evapotranspiration regime seems
to draw the opposite of the trend of that of rain From this analysis, we can
say that the rainy season in the area is from June to September with more than
80% of the year to date.

Figure 3: Rainfall and ETP regime

Daily precipitation levels show (Figure 4) very variable frequencies over the
study period (1990-2019). Daily heights of less than 5 mm are very frequently
observed in the area. However, precipitation greater than 40 mm although very
rare was also observed.

Figure 4 : Rainfall occurancy

####  4.1.3 Seasonal humidity variability

Change in relative humidity over the season (Figure 5) is shown as an
asymmetric bell curve to the right. Minimum humidity is below 25% during the
months of November to April and is above 40% between July and September.
Maximum humidity has the same evolution as the minimum with rates below 60% in
the dry season and above 80% between July and September. However, it should be
noted that the two types of humidity behave in an opposite way in terms of
dispersion. In fact, the minimum humidity has a strong dispersion in the wet
season and a weak dispersion between October and April, unlike the maximum
humidity which shows a weak dispersion between July and October and a strong
dispersion the rest of the year.

Figure 5 : Relative humidity regime

####  4.1.4 Seasonal wind speed variability

Wind speed is generally between 0.5 m/s and 2.5 m/s (Figure 6). Strong winds
are recorded during the months of January, February, May and June with a speed
between 1.5 m/s and 2.5 m/s. This period also records a strong dispersion of
the wind speed which could be explained by the installation of the rainy
season with the predominance of thunderstorms for the months of May and June
and the cold season which extends from December to February.

Figure 6 : Wind speed regime

The strong monthly wind speed occurrences are winds below 2 m/s. Strong winds
of 2 to 3 m/s are less frequent in the area (Figure 7).

Figure 7 : Wind speed occurancy

####  4.1.5 Seasonal insolation variability

Monthly average insolation and its monthly cumulative (Figure 8) experience
the same type of dispersion during the year and over the period 1990-2019. In
fact, there is little dispersion during the wet season and great dispersion
the rest of the year. Monthly averages vary between 6 and 9 hours and the
cumulatives are generally greater than 225 hours per month. The smallest
averages are observed during the months of July, August and September due to
the frequency of cloud cover during the rainy season which prevents solar
radiation from reaching the earth’s surface and the strongest in the dry
season.

Figure 8 : Insolation regime

###  4.2 Inter annual variability of agro-climatic parameters

####  4.2.1 Inter annual temperature variability

Inter annual variation of temperatures (Figure 9) shows an increasing trend
for both minimum and maximum. However, the increase is greater but not
significant (Mann-Kendal test used with a significance level at the 5%
threshold) for minimum temperatures (0.14∘C / 10 years) than maximum
temperatures (0.12∘C /10 years). The Pettitt test at the 5% threshold made it
possible to detect the years of ruptures. For the minimum temperatures, the
year 2001 was detected as the year of rupture with a p-value$<$0.05 on the
other hand there is no significant rupture for the maximum temperatures. The
difference between the rupture periods (before 2001 and after 2001) is of the
order of 0.4∘C. Temperature rise becomes continuous and more significant on
minimum temperatures than on maximum temperatures.

Figure 9: Temperature inter annual variabilities

Figure 10 : Temperature anomalies

Overall, the analysis of minimum, maximum and average temperature anomalies
(Figures 9;10) shows a succession of hot decades (1990-1999 and 2010-2019) and
cold decades (2000-2009) throughout the study period. If this trend continues,
we could see more warm years in the next 10 years.

Figure 11: Mean temperature anomalies

####  4.2.2 Inter annual rainfall variability

The evolution of annual cumulative and annual maximum rainfall (Figure 12) is
marked by a high variability in cumulatives (500- and over 1000 mm) with an
increasing trend which, however, is not significant in the Mann-Kendall test.
Also the Pettitt test was found to be insignificant with p-values $>$0.05.
Nevertheless the difference in average calculated between rupture periods is
worth approximately 75 mm. In other words, the last 10 years have seen a drop
in rainfall amounts of around 75 mm. In the distribution of annual maximum
rainfall, 50% of the rainfall is greater than 65 mm and it exceeds 73 mm every
8 out of 30 years.

Figure 12: Rainfall variability

The analysis of the variability of cumulative rainfall (Figure 13) shows a
marked trend of years in deficit over the period 2001-2008 and two periods of
excess years: a shorter one from 1992 to 1995 and a longer one from 2009 to
2019. The remarkable deficit years are 2002 and 2007.

Figure 13: Cumulative rainfall variability

####  4.2.3 Inter annual wind speed variability

Annual wind speed averages are between 0.5 and 2 m / s (Figure 14). They show
a downward trend. Maximum daily wind speed exceeds in more than 75% of the
observations the 04 m / s.

Figure 14: Inter annual wind speed evolution

####  4.2.4 Inter annual insolation and humidity variability

Insolation shows a decreasing trend and varies between 6 and 8.5 hours on
average per year (Figure 15). Relative humidity shows low values for 2002 and
2007, which correspond to deficit years .

Figure 15: Insolation and humidity evolution

####  4.2.5 Inter annual evapotranspiration variability

The overall trend of evapotranspiration calculated from the Penman-Monteith
formula is slightly downward over the period and despite a high variability
observed, the trend and break tests were not significant (Figure 16). However,
the difference observed between the two averages of the two periods is of the
order of 0.22 mm.

Figure 16: ETP inter annual evolution

###  4.3 Characterization of the rainy season

####  4.3.1 Analysis of the onset date of the rainy season

The distribution of the onset dates of the season is symmetrical and presents
a large dispersion and a very high variability (Figure 17). They take place
between April 12 at the earliest and June 18 at the latest. The average onset
date of the season is around 07 June. In 75% of cases, the onset dates are
less than June 15. The farmer runs less risk of re-sowing by sowing after this
date because it is observed 3 years out of 4. Also any farmer who sows before
May 15th runs a great risk of re-sowing because its occurrence is 1 year out
of 4 .

Figure 17: Onset dates evolution

The late start of the season dates are marked by high variability over the
entire study period and remain strong from 2015 onwards (Figure 18).

Figure 18: Onset dates anomalies

####  4.3.2 Analysis of the cessetion dates of the rainy season

The cessetion dates of season vary between September 10 and October 20 for the
Sivakumar criterion and between September 2 and October 17 for the PRESAO
criterion(Figure 19). In addition, the two criteria meet on average with dates
around October 30 for Sivakumar and October 28 for PRESAO.

Figure 19: Cessetion dates variabilities

Overall, the end of season dates calculated with the SIVAKUMAR criterion are
slightly higher than those obtained by the PRESAO criterion. Indeed, 3 years
out of 4, the dates for SIVAKUMAR are lower than October 11 and those of
PRESAO are lower than October 06 (Figure 20).

Figure 20: Cessetion dates distribution

####  4.3.3 Analysis of the length of the rainy season

In rainfed cultivation, knowing the length of the rainy season in an area is
very crucial in order to make the right decision regarding the variety choice
for the success of a given crop. Thus for the two criteria retained in this
study, the lengths of the season vary approximately between 60 and 165 days.
For an average season length of 125 days, the use of variety of 100 to 120
days is recommended for the study area (Figure 21).

Figure 21: Season length variabilities

For both criteria, the tendency for long or short season lengths is very
variable. Nevertheless, the tendency of long seasons seems more marked for the
PRESAO criterion (Figure 22).

Figure 22: Season length anomalies

####  4.3.4 Analysis of the number of rainy days in the season

The number of rainy days in the season varies between 20 and 40 days except in
extreme cases against 40 and 60 days in the year (Figure 23). Thus, the number
of rains on the fringes of the rainy season is less important in the area.

Figure 23: Rainy dates evolution

###  4.4 Correlation analysis

####  4.4.1 Agro-climatic parameters correlation

In general, variables of the same type show a strong correlation with each
other (Figure 24). Indeed, the groups of variables such as (Tmax, Tmin, Tmean)
and (Umax, Umin) naturally exhibit a strong positive correlation. On the other
hand, it is important to underline certain more or less strong connections but
in the opposite direction. This is the case, for example, with maximum
temperature and relative humidity. Thus, hotter it is, less is water vapor in
the atmospheric air. The inverse correlation between minimum temperature and
insolation can be explained by the fact that the minimum temperature is
recorded at 06h UT and therefore before the appearance of sunlight. The
negative correlation between maximum temperature and rain is explained by the
fact that the occurrence of precipitation reduces temperatures in general. The
paradox with the minimum temperatures and the rain at this station could be
explained by the fact that most of the rains occur during the day and
therefore have more impact on the maximum temperatures than the minimum. Wind
and relative humidity correlate in the opposite direction since the wind is a
factor of evaporation and therefore participates in the drying of air masses.

Figure 24: Ago-climatic parameters correlation

####  4.4.2 Characteristic parameters of the season correlation

In general, there is a significantly negative correlation between the onset
date and the other parameters characteristic of the rainy season (Figure 25).
A result which corroborates the results of previous studies as in [13].
Moreover, the PRESAO method for the end of season date seems to contradict
this assertion on our data by showing a significantly positive correlation
with the onset date. This last result justifies our choice for the Sivakumar
method in the rest of our study.

Figure 25: Caracteristic parameters of season correlation

###  4.5 Identification of the probability distributions of agro-climatic
risks

The normal law seems to be the law of probability which best adjusts the
various key parameters of the rainy season (onset and cessetion date of the
season, length of the season). Indeed, the Shapiro-Wilks tests and the
quantile-quantile plots lead to the same conclusions in general.

####  4.5.1 Onset dates probability distribution

Figure 26: Onset dates adequation

####  4.5.2 Length of the season probability distribution

##### Sivakumar criterion

Figure 27: Length of season adequation (Sivakumar)

##### PRESAO criterion

Figure 28: Length of season adequation (PRESAO)

####  4.5.3 Cessetion dates probability distribution

##### Sivakumar criterion

Figure 29: Cessetion dates adequation (Sivakumar)

##### PRESAO criterion

Figure 30: Cessetion dates adequation (PRESAO)

####  4.5.4 Dry sequences distribution

Over a period of 30 years, the most recurrent dry sequences in the area have a
length less than or equal to 7 days (Figure 31). They constitute 93% of the
dry sequences observed. However, dry sequences of greater duration are also
observed. This last result is probably due to the false start observed over
the years 2002, 2005, 2013 and 2019.

Figure 31: Dry sequences distribution

##  5 Conclusion

This study allowed us to contribute to the analysis of the evolution of agro-
climatic risks in a context of climate variability in the region of Ségou in
Mali. So, we have used statistical tests on variables characterizing the rainy
season and we identified the trends and variabilities contained in the data
series. We have proposed a statistical modeling of these different agro-
climatic risks. Tthe seasonal variability of agro-climatic parameters
(temperature, rainfall , humidity, wind speed and insolation) have been
analized as well as their inter annual variability. Furthermore the
probability distributions of agroclimatic risks have been identified and the
characterization of the rainy season was clarified.

## References

  * Assemian et al. [2013] E. Assemian, F. Kouame, É. Djagoua, K. Affian, J. P. R. Jourda, M. Adja, T. Lasm, and J. Biemi.  Étude de l’impact des variabilités climatiques sur les ressources hydriques d’un milieu tropical humide: cas du département de bongouanou (est de la côte d’ivoire).  _Revue des sciences de l’eau/Journal of Water Science_ , 26(3):247–261, 2013.  doi: 10.7202/1018789ar.
  * Bazzaz and Sombroek [1997] F. Bazzaz and W. Sombroek.  Changements du climat et production agricole.  _FAO et Polytechnica, pp1-10_ , 1997.
  * Bell and Lamb [2006] M. A. Bell and P. J. Lamb.  Integration of weather system variability to multidecadal regional climate change: The west african sudan-sahel zone, 1951-98.  _Journal of Climate_ , 19(20):5343–5365, 2006\.  ISSN 08948755, 15200442.  doi: https://doi.org/10.1175/JCLI4020.1.
  * Benoît [2008] É. Benoît.  Les changements climatiques: vulnérabilité, impacts et adaptation dans le monde de la médecine traditionnelle au burkina faso.  _VertigO-la revue électronique en sciences de l’environnement_ , 8(1), 2008.  doi: https://doi.org/10.4000/vertigo.1467.
  * Berg et al. [2013] A. Berg, N. de Noblet-Ducoudré, B. Sultan, M. Lengaigne, and M. Guimberteau.  Projections of climate change impacts on potential c4 crop productivity over tropical regions.  _Agricultural and Forest Meteorology_ , 170:89–102, 2013\.  doi: https://doi.org/10.1016/j.agrformet.2011.12.003.
  * DIOUF et al. [2001] M. DIOUF, A. NONGUIERMA, A. AMANI, and A. Royer.  Lutte contre la sécheresse au sahel: résultats, acquis et perspectives au centre régional agrhymet.  _Science et changements planétaires/Sécheresse_ , 11(4):257–66, 2001.  ISSN 1147-7806.
  * Le Barbé et al. [2002] L. Le Barbé, T. Lebel, and D. Tapsoba.  Rainfall variability in west africa during the years 1950–90.  _Journal of climate_ , 15(2):187–202, 2002.  doi: https://doi.org/10.1175/1520-0442(2002)015<0187:RVIWAD>2.0.CO;2.
  * Lubes et al. [1994] H. Lubes, J. Masson, E. Servat, J. Paturel, B. Kouame, and J. Boyer.  Caractérisation de fluctuations dans une série chronologique par application de tests statistiques.  _Etude bibliographique, Programme ICCARE, Rapport_ , 3:1–21, 1994.  URL https://www.documentation.ird.fr/hor/fdi:010020654.
  * Mann [1945] H. B. Mann.  Nonparametric tests against trend.  _Econometrica: Journal of the econometric society_ , pages 245–259, 1945.  doi: https://doi.org/10.2307/1907187.  URL https://www.jstor.org/stable/1907187.
  * Moorcroft et al. [2006] P. Moorcroft, S. W. Pacala, and M. Lewis.  Potential role of natural enemies during tree range expansions following climate change.  _Journal of Theoretical Biology_ , 241(3):601–616, 2006.  doi: https://doi.org/10.1016/j.jtbi.2005.12.019.
  * Nicholson [2001] S. Nicholson.  Climatic and environmental change in africa during the last two centuries.  _Climate Research - CLIMATE RES_ , 17:123–144, 08 2001\.  doi: 10.3354/cr017123.
  * Roudier et al. [2012] P. Roudier, B. Sultan, P. Quirion, C. Baron, A. Alhassane, S. B. Traoré, and B. Muller.  An ex-ante evaluation of the use of seasonal climate forecasts for millet growers in sw niger.  _International Journal of Climatology_ , 32(5):759–771, 2012.  doi: https://doi.org/10.1002/joc.2308.  URL https://rmets.onlinelibrary.wiley.com/doi/abs/10.1002/joc.2308.
  * Sivakumar [1988] M. Sivakumar.  Predicting rainy season potential from the onset of rains in southern sahelian and sudanian climatic zones of west africa.  _Agricultural and forest meteorology_ , 42(4):295–305, 1988.  doi: https://doi.org/10.1016/0168-1923(88)90039-1.
  * Sivakumar [1992] M. V. K. Sivakumar.  Empirical analysis of dry spells for agricultural applications in west africa.  _Journal of Climate_ , 5(5):532–539, 1992.  ISSN 08948755, 15200442.  doi: https://doi.org/10.1175/1520-0442(1992)005<0532:EAODSF>2.0.CO;2.  URL http://www.jstor.org/stable/26197057.
  * Solomon et al. [2007] S. Solomon, M. Manning, M. Marquis, D. Qin, et al.  _Climate change 2007-the physical science basis: Working group I contribution to the fourth assessment report of the IPCC_ , volume 4.  Cambridge university press, 2007.
  * Sultan [2012] B. Sultan.  Global warming threatens agricultural productivity in africa and south asia.  _Environmental Research Letters_ , 7(4):041001, 2012.  doi: 10.1088/1748-9326/7/4/041001.  URL https://doi.org/10.1088/1748-9326/7/4/041001.
  * Sultan et al. [2005] B. Sultan, C. Baron, M. Dingkuhn, B. Sarr, and S. Janicot.  Agricultural impacts of large-scale variability of the west african monsoon.  _Agricultural and forest meteorology_ , 128(1-2):93–110, 2005.  doi: https://doi.org/10.1016/j.agrformet.2004.08.005.