now at: ]Faculty of Physics, Vienna University, Boltzmanngasse 5, 1090 Wien,
Austria

# Shortcut-to-Adiabatic Controlled-Phase Gate in Rydberg atoms

Luis S. Yagüe Bosch  Institute of Theoretical Physics, Heidelberg University,
Philosophenweg 12, 69120 Heidelberg, Germany  Tim Ehret  [  Institute of
Theoretical Physics, Heidelberg University, Philosophenweg 12, 69120
Heidelberg, Germany  Francesco Petiziol  Technische Universität Berlin,
Institut für Theoretische Physik, Hardenbergstrasse 36, Berlin 10623, Germany
Ennio Arimondo  Dipartimento di Fisica, Università di Pisa, Largo Pontecorvo
3, 56127 Pisa, Italy  Istituto Nazionale di Ottica, Consiglio Nazionale delle
Ricerche, Universitá di Pisa, Largo Pontecorvo 3, 56127 Pisa, Italy  Sandro
Wimberger<EMAIL_ADDRESS>Dipartimento di Scienze
Matematiche, Fisiche e Informatiche, Università di Parma, Parco Area delle
Scienze 7/A, 43124 Parma, Italy  INFN, Sezione di Milano Bicocca, Gruppo
Collegato di Parma, Parco Area delle Scienze 7/A, 43124 Parma, Italy

###### Abstract

A shortcut-to-adiabatic protocol for the realization of a fast and high-
fidelity controlled-phase gate in Rydberg atoms is developed. The adiabatic
state transfer, driven in the high-blockade limit, is sped up by compensating
nonadiabatic transitions via oscillating fields that mimic a counterdiabatic
Hamiltonian. High fidelities are obtained in wide parameter regions. The
implementation of the bare effective counterdiabatic field, without original
adiabatic pulses, enables to bypass gate errors produced by the accumulation
of blockade-dependent dynamical phases, making the protocol efficient also at
low blockade values. As an application toward quantum algorithms, how the
fidelity of the gate impacts the efficiency of a minimal quantum-error
correction circuit is analyzed.

##  I Introduction

Rydberg atoms [1] have been established as a platform for many applications in
quantum physics over the last 30 years, going from old beam experiments [2, 3,
4, 5] to much better controlled cold-atom setups in which single Rydberg atoms
are hold in position with optical tweezers [6, 7, 8, 9]. We are interested in
the use of cold Rydberg atoms for quantum information protocols [10, 11, 12,
6], in view of their ability to implement single well-controlled quantum gates
with high fidelity [13, 14, 7, 6, 15, 16, 17, 18, 19].

In particular, adiabatic protocols seem to be promising for large-scale
computations because of their intrinsic robustness [20, 21, 22]. In contrast
to sudden or quasi-resonant protocols no precise resonance tuning and pulse
timing is required [10, 11, 12, 13, 14, 7, 6, 15, 16, 17, 18, 19, 23]. A
specific adiabatic protocol for realizing a controlled-phase (CZ) gate has
recently been proposed by Saffman et al. [24]. An advantage of this adiabatic
gate is that it prevents large transient occupation of the Rydberg states,
such that the protocol is less exposed to spontaneous emission from such
states—a primary source of fidelity loss in this platform. This interacting
Rydberg gate works differently from a recent realization of a CZ gate with
superconducting qubits based on a different form of interaction [25].

In this work, we explore the possibility to overcome the typical slowness of
adiabatic quantum computation protocols by accelerating the CZ gate of [24]
via shortcuts to adiabaticity [26, 27] (STA). The idea of a shortcut is to
find a modified Hamiltonian that realizes the same population transfer with
high fidelity but in a nonadiabatic, faster way, and it has been successfully
applied to control in few-level systems in multiple experimental platforms,
see, e.g., [28, 29, 30, 31, 32]. In particular, we derive and analyse a
control protocol in which nonadiabatic transitions are counteracted via
additional oscillating components in the control fields that mimic a
counterdiabatic (CD) Hamiltonian [33, 31, 34]. While STA techniques typically
address individual state-to-state transfers, we here rather engineer the full
propagator by targeting the realization of a complete gate. This more
challenging task, which has been traditionally carried out using optimal
control [35, 36, 37, 38, 39], was recently also addressed with STA methods
[40, 41, 42, 43, 44, 45, 46, 47].

The focus of the CD protocols is usually on non-interacting systems, because
as a common feature the required CD corrections are hard to be implemented
experimentally, with a few exceptions of integrable many-body systems [26, 48,
27, 49, 50]. This is true, in principle, also for protocols with Rydberg atoms
interacting in their excited states, but we demonstrate that an approximate
Hamiltonian not requiring many-body controls nevertheless leads to good final
fidelities for a broad range of protocol times and system parameters.

##  II Adiabatic CZ gate

While the CD approach may be applied to any adiabatic protocol, we focus our
attention on the configuration depicted in Fig. 1 representing an adiabatic
Rydberg driving scheme slightly modified with respect to [24]. The two atoms,
having hyperfine levels $\ket{0}$ and $\ket{1}$ of the electronic ground
state, are symmetrically driven with lasers of Rabi frequency $\Omega(t)$
targeting the transition from $\ket{1}$ to a Rydberg state $\ket{r}$, with
detuning $\Delta(t)$, and they experience a Rydberg-Rydberg interaction of
strength $V$. The time-dependent Hamiltonian, in the frame rotating at the
laser frequency, reads

$H(t)=H_{d}(t)\otimes\mathds{1}+\mathds{1}\otimes
H_{d}(t)+V\ket{rr}\\!\bra{rr},$ (1)

where $H_{d}(t)$ describes single atom control,

$H_{d}(t)=\frac{\Omega(t)}{2}(\ket{r}\\!\bra{1}+\ket{1}\bra{r})+\Delta(t)\ket{r}\\!\bra{r}.$
(2)

Both pulses $\Omega(t)$ and $\Delta(t)$ are represented in Fig. 1(a). The Rabi
frequency has the following shape typical for adiabatic rapid passage [24, 46,
51, 21]

$\Omega(t)=\frac{\Omega_{\mathrm{max}}}{\mathcal{N}}\left[e^{-(t-t_{0})^{4}/\tau^{4}}-a-bt(t-2t_{0})\right],$
(3)

which is centered at $t_{0}$. The real constants $a$ and $b$ are chosen such
that the function and its derivative vanish at the start and end points.
$\mathcal{N}$ is a normalization constant such that $\Omega_{\mathrm{max}}$ is
the maximum of the function $\Omega(t)$. The detuning is sine shaped including
a continuous sign change in phase II to avoid a sudden jump of the function.
These pulses produce subsequent excitation (phase I) and de-excitation (phase
III) of the atoms to state $\ket{r}$, where the Rydberg blockade mechanism
sets in, as sketched in Fig. 1(b). During phase II, the detuning is smoothly
inverted, while the state remains unchanged. The instantaneous energy levels
of $H(t)$ as a function of time are shown in Fig. 1(c), where two avoided
crossings are visible (indicated by arrows), signaling a mixing of states
$\ket{01}$ and $\ket{0r}$ (red colour), $\ket{10}$ and $\ket{r0}$ (orange
colour) and $\ket{11}$ and $\ket{d_{+}}=(\ket{1r}+\ket{r1})/\sqrt{2}$ (blue
colour), respectively. The pulses are designed for initial states $\ket{01}$,
$\ket{10}$ and $\ket{11}$ to accumulate a geometric phase $\pi$ during the de-
excitation in the infinite blockade and adiabatic limit. Indeed, in this case
the basis states evolve according to

$\displaystyle\ket{1}\to\ket{r}\to-\ket{1},$ (4)
$\displaystyle\ket{11}\to\ket{d_{+}}\to-\ket{11}.$ (5)

Together with the fact that state $\ket{00}$ has no dynamics induced by the
control fields, the phase accumulation produces the CZ gate. An advantage of
the adiabatic gate is that it drives the atoms almost continuously, which has
been observed to be particularly robust against spontaneous emission from the
Rydberg state compared to having the atoms sit in state $\ket{r}$ [24].

Figure 1: (a) Time dependence of the adiabatic pulses. In phase II the
detuning is changed continuously. The explicit form of the curves is given in
the text. The peak values of the pulses are $\Omega_{\mathrm{max}}/2\pi=17$
MHz and $\Delta_{\mathrm{max}}/2\pi=23$ MHz at a total time of
$T_{\mathrm{tot}}=0.594\;\mu$s. The width is given by $\tau=0.0945\;\mu$s and
the detuning in phase II takes $\approx T_{\mathrm{tot}}/10$. (b) Energy
levels scheme and sketch of the driving sequence; the atoms are first excited
to their Rydberg state (phase I), and are then de-excited (phase III); the
adiabatic pulses interpolate the three phases in a continuous manner. (c)
Instantaneous energy levels of the driven Hamiltonian of Eq. (1) as a function
of time in phase I; Arrows indicate two avoided crossings between states
$\ket{01}$ and $\ket{0r}$ or $\ket{10}$ and $\ket{r0}$, and states $\ket{11}$
and $\ket{1r}+\ket{r1}$, respectively.

The goal of the present work is to develop a new driving protocol that
produces the same adiabatic evolution as the original protocol, but in a
shorter time, or with a higher fidelity in a given total duration. In order to
construct such a protocol, we apply the idea of effective counterdiabatic
driving (eCD) [33] based on Floquet-engineering [52, 53, 54]: the pulses are
modified with the introduction of additional Fourier components, whose
amplitude is chosen such that they compensate for nonadiabatic transitions, by
mimicking the effect of a counterdiabatic field [55, 56]. Differently from
other STA techniques, which aim at matching the adiabatic dynamics only at the
initial and final times [57, 58, 59, 60, 61], the eCD method strives for
tracking the adiabatic dynamics closely at all times. This is an advantage for
the present protocol since it allows to maintain the low occupation of the
Rydberg states guaranteed by the adiabatic sweep. A similar situation is
encountered for accelerated Stimulated Raman Adiabatic Passage (STIRAP) [46,
32], where adiabatic following of a dark state protects against emission from
the intermediate state. In the following, the states $\ket{00}$ and
$(\ket{r1}-\bra{1r})/\sqrt{2}$ will be neglected, since they are not involved
in the driven dynamics. Moreover, we consider the regime of large blockade,
$V\gg\Omega(t),\Delta(t)$, such that the $\ket{rr}$ state can be adiabatically
eliminated (please see Appendix A). In our numerical simulations, all states
are included (without adiabatic elimination) to gauge the effect of finite
values of $V$.

##  III Accelerated CZ gate

In order to construct the accelerating protocol, the first step is to
determine the exact counterdiabatic field [55, 56, 57, 27] for the Hamiltonian
of Eq. (1). This is defined as the correction $H_{\mathrm{CD}}(t)$ to the
Hamiltonian, such that $H(t)+H_{\mathrm{CD}}(t)$ generates the adiabatic
dynamics exactly without restrictions on the total evolution time. Indicating
with $E_{n}(t)$ and $\ket{n(t)}$ the instantaneous eigenvalues and
eigenvectors of the adiabatically driven Hamiltonian $H(t)$, this means that
an initial eigenstate $\ket{n(0)}$ is mapped by the evolution to the
corresponding instantaneous eigenstate $\ket{n(t)}$ for all times, with
additional dynamical and geometrical phase factors,

$\exp\left(-i\int_{0}^{t}dt^{\prime}E_{n}(t^{\prime})+\int_{0}^{t}dt^{\prime}\braket{n(t^{\prime})}{\partial_{t^{\prime}}n(t^{\prime})}\right)\ket{n(t)}.$
(6)

When instantaneous energies and eigenvectors are known, and in the absence of
degeneracies, the counterdiabatic field can be compactly expressed as [55]

$H_{\mathrm{CD}}(t)=i\sum_{n,m\neq
n}\frac{\ket{n(t)}\bra{n(t)}\partial_{t}H(t)\ket{m(t)}\bra{m(t)}}{E_{m}(t)-E_{n}(t)}.$
(7)

We remark that, if $H_{\mathrm{CD}}(t)$ is implemented, it ensures both
tracking of the instantaneous eigenstate $\ket{n(t)}$ and fixes the geometric
phase factor in Eq. (6), while the role of $H(t)$ is to imprint the dynamical
phase only. Since the gate is based on the accumulation of only geometric
phases, $H_{\mathrm{CD}}(t)$ would be sufficient to realise the CZ gate, even
in the absence of the original adiabatic driving $H(t)$. This is an
interesting property, since in the present protocol, dynamical phases are
actually responsible for unwanted phase shifts: the state $\ket{11}$
dynamically accumulates an excess phase [62, 63]

$\varphi=\int_{0}^{T_{\mathrm{tot}}}dt\;\Omega^{2}(t)/4V.$ (8)

This phase is non-negligible for large Rabi frequencies (which would actually
be the ideal regime for adiabaticity) and finite $V$. Therefore, by
implementing only $H_{\mathrm{CD}}(t)$, rather than $H(t)$ or
$H(t)+H_{\mathrm{CD}}(t)$, one can obtain both a perfect shortcut and avoids
such phase errors.

To derive $H_{\mathrm{CD}}(t)$ for the Hamiltonian (1), we first write $H(t)$
as the sum $H(t)=[H_{0}(t)\otimes\ket{0}\\!\bra{0}+\ket{0}\\!\bra{0}\otimes
H_{0}(t)]+H_{1}(t)$ of two mutually commuting terms corresponding to the
subspaces in which at least one of the two atoms is in its $\ket{0}$ ground
state (the term in square brackets) or neither is [$H_{1}(t)$], respectively.
Adopting the bases $\\{\ket{10},\ket{r0},\ket{01},\ket{0r}\\}$ and
$\\{\ket{11},\ket{d_{+}}\\}$ for the two separate subspaces, see Appendix A,
the Hamiltonians $H_{k}(t)$ read

$H_{k}(t)=\frac{1}{2}\begin{bmatrix}0&\Omega_{k}(t)\\\
\Omega_{k}(t)&2\Delta(t)\end{bmatrix},$ (9)

with $\Omega_{0}(t)=\Omega(t)$ and $\Omega_{1}(t)=\sqrt{2}\Omega(t)$. The
Hamiltonian $H(t)$ hence describes two uncoupled two-level problems, one
defined in each of the {$\ket{1},\ket{r}$} single-atom subspaces, and the
other associated with the two-qubit subspace $\\{\ket{11},\ket{d_{+}}\\}$.
These two two-level problems share the same detuning $\Delta(t)$ but the
latter has a Rabi frequency which is enhanced by a factor of $\sqrt{2}$, an
effect due to the Rydberg blockade [24, 23, 7]. Following this observation,
the counterdiabatic field $H_{\mathrm{CD}}(t)$ can then readily be expressed
using known results for two-level systems [55, 56], and, for real-valued
$\Omega_{k}(t)$ and $\Delta(t)$, it involves purely imaginary off-diagonal
couplings. It reads
$H_{\mathrm{CD}}(t)=H_{\mathrm{CD},0}(t)+H_{\mathrm{CD},1}(t)$, with

$\displaystyle H_{\mathrm{CD},0}(t)=$
$\displaystyle-\frac{f_{0}(t)}{2i}\big{[}\ket{1}\\!\bra{r}\otimes\ket{0}\\!\bra{0}$
$\displaystyle\quad+\ket{0}\\!\bra{0}\otimes\ket{1}\\!\bra{r}\big{]}+\mathrm{H.c.}\
,$ $\displaystyle H_{\mathrm{CD},1}(t)=$
$\displaystyle-\frac{f_{1}(t)}{2i}\ket{11}\\!\bra{d_{+}}+\mathrm{H.c.}\ ,$
(10)

where the time-dependent amplitudes $f_{k}(t)$, $k\in\\{0,1\\}$, are given by

$f_{k}(t)=\frac{\Delta(t)\partial_{t}\Omega_{k}(t)-\Omega_{k}(t)\partial_{t}{\Delta(t)}}{\Delta^{2}(t)+\Omega_{k}^{2}(t)}.$
(11)

The time profile of the these functions is shown in Fig. 2(a,b) for the
adiabatic pulse shapes considered here. A direct implementation of the CD
field would be experimentally challenging for two reasons. First, it requires
the introduction of additional pulses with imaginary Rabi frequency. Second,
it requires time-dependent control of two-qubit interactions, since different
pulses should act conditioned on whether one of the atoms is de-excited or
whether they are both excited. This is necessary in order to counteract
unwanted nonadiabatic excitations that can drive the system into entangled
states in the $\\{\ket{11},\ket{d_{+}}\\}$ subspace. We here solve the first
issue by means of the eCD driving technique [33]: the imaginary part of the
pulses is approximated by introducing additional real Fourier components in
the original Hamiltonian, oscillating quickly at a frequency $\omega\gg
f_{k}(t)$. The amplitude and phase of such oscillations is found by enforcing
the evolution operator at the end of the oscillation period to match, up to
first order in $1/\omega$, the exact counterdiabatic propagator (see Appendix
B). In other words, the eCD oscillations approximate the Hamiltonian
$H_{\mathrm{CD}}(t)$ with an error scaling with $O(\omega^{-3/2})$ [33].

A suitable eCD Hamiltonian is found to have the form

$\displaystyle H_{\mathrm{eCD}}(t)=$
$\displaystyle\big{[}g_{1}(t)(\ket{r}\\!\bra{1}\otimes
P_{0}+P_{0}\otimes\ket{r}\\!\bra{1})$
$\displaystyle+g_{2}(t)(\ket{r}\\!\bra{1}\otimes
P_{1}+P_{1}\otimes\ket{r}\\!\bra{1})+\mathrm{H.c.}\big{]}$
$\displaystyle+g_{3}(t)(P_{r}\otimes\mathds{1}+\mathds{1}\otimes P_{r}),$ (12)

with $P_{x}=\ket{x}\\!\bra{x}$, for $x=0,1,r$, and with oscillating functions

$\displaystyle g_{1}(t)=$ $\displaystyle\sqrt{\omega f_{0}(t)}\sin(\omega
t),\quad g_{2}(t)=\sqrt{\omega f_{1}(t)}\cos(\omega t),$ $\displaystyle
g_{3}(t)=$ $\displaystyle-\sqrt{\omega f_{1}(t)}\sin(\omega t)+\sqrt{\omega
f_{0}(t)}\cos(\omega t),$ (13)

To investigate the performance of this shortcut-to-adiabatic protocol, we
first analyse the final infidelity
$\mathcal{I}=1-|\\!\braket{\psi_{\mathrm{fin}}}{\psi(t_{f})}\\!|^{2}$ produced
starting from an initial state $\ket{\psi}=(\ket{00}+\ket{11})/\sqrt{2}$,
which should ideally be converted by the CZ gate into the final state
$\ket{\psi_{\mathrm{fin}}}=(\ket{00}-\ket{11})/\sqrt{2}$. This is compared
with the original adiabatic protocol in Fig. 2, where the infidelity is
depicted as a function of the total protocol time and of the Rydberg blockade
for the adiabatic $H(t)$ [Fig. 2(c)] and shortcut-to-adiabatic
$H_{\mathrm{eCD}}(t)$ [Fig. 2(d)] protocols. Fidelity regions of interest,
above $1-\mathcal{I}=0.999$, require the adiabatic protocol to work in the
parameter region of both large blockade and protocol time. When reducing the
protocol time for the adiabatic gate, the pulse amplitudes need to be
increased in order to maintain the same area. However, this is accompanied by
an increase of the phase error mentioned in Sec. II, causing significant
fidelity loss when accelerating the gate with finite blockade. On the
contrary, the eCD drive bypasses the generation of the unwanted phase, since
it implements only $H_{\mathrm{CD}}(t)$. It is thus able to attain large
fidelities for most of the parameter space explored, and in any case
outperforms the purely adiabatic scheme. This confirms the efficiency of the
modified protocol, on the one hand, and shows also that it is very stable with
respect to the choice of timing and blockade parameters. Considering as an
example the value of the Rydberg interaction of $V=100$ MHz, a fidelity above
0.998 can be reached in $T_{\mathrm{tot}}=0.26\;\mu$s, for pulse parameters
$\Omega_{\mathrm{max}}/2\pi=17$ MHz, $\Delta_{\mathrm{max}}/2\pi=23$ MHz, and
a driving frequency $\omega=300$ MHz.

Similarly to the adiabatic protocol, also the modified one necessitates a
minimal total protocol time to access a given high-fidelity region. This
phenomenon, reminiscent of quantum speed limits [64], descends indirectly from
the increasing violation of adiabaticity for shorter durations. Indeed, as the
protocol time is shortened, the exact CD field compensating for nonadiabatic
transitions varies on a faster timescale. Specifically, the time derivative of
$\partial_{t}H_{\mathrm{CD}}(t)$ is ruled by $\partial_{t}f_{0}(t)$ and
$\partial_{t}f_{1}(t)$, which scale like $1/T^{2}_{\mathrm{tot}}$ according to
Eq. (11). As a consequence, a larger driving frequency $\omega$ would be
needed to approximate ${H}_{\mathrm{CD}}$ with sufficient precision to
maintain a constant fidelity while the protocol time is reduced. Similar
effects, in the context of the Landau-Zener-Majorana-Stückelberg problem, have
been studied, e.g., in Ref. [33].

To fully characterize the performance of the gate, an experiment could perform
a full quantum process tomography (QPT) [65, 66]. According to Kraus’ theorem,
the qubit gate dynamics can be described as a quantum map
$\mathcal{E}(\cdot)=\sum_{j}K_{j}(\cdot)K_{j}^{\dagger},$ where the argument
is the system’s (pure or mixed) state, with Kraus operators $K_{j}$ satisfying
$\sum_{j}K_{j}^{\dagger}K_{j}=\mathds{1}$. Expanding the operators $K_{j}$ on
tensor products of single-qubit Pauli matrices,
$\\{\sigma_{i}\\}_{i=0,\ldots,3}=\\{\mathds{1},\sigma_{x},\sigma_{y},\sigma_{z}\\}$,
one obtains

$\mathcal{E}(\cdot)=\sum_{i}K_{i}(\cdot)K_{i}^{\dagger}=\sum_{ij}\chi_{ij}A_{i}(\cdot)A_{j}^{\dagger},$
(14)

which defines the $\chi$-matrix fully characterizing the quantum process [66],
where
$\\{A_{i}\\}=\\{\sigma_{j}\otimes\sigma_{k}\\}_{j=0,\dots,3;k=0,\dots,3}$.
Figure 2(e) shows the deviation of the ideal $\chi$-matrix from the here
proposed eCD-based CZ gate. While simulations include also the Rydberg states,
the process tomography shown in Fig. 2(e) is computed in the computational
subspace of the two qubits. The matrix elements of the eCD $\chi$-matrix
differ from the ideal case by less than $10^{-3}$, confirming the
effectiveness of the eCD gate for all possible input states. The QPT was
computed using the QuTip Python library [67, 68].

Figure 2: Time profile of the amplitudes $f_{0}(t)$ [(a)] and $f_{1}(t)$ [(b)]
of the counterdiabatic field $H_{\text{CD}}$. In phase II they vanish as the
Rabi pulses are exactly zero. Infidelity for the initial state
$(\ket{00}+\ket{11})/\sqrt{2}$ of the adiabatic protocol (c) and of the CD
protocol (d) as a function of the total protocol time $T_{\mathrm{tot}}$ and
of the strength of the Rydberg blockade. The change of the detuning in phase
II takes $\approx T_{\mathrm{tot}}/10$. The pulses have peak values
$\Omega_{\mathrm{max}}/2\pi=17$ MHz and $\Delta_{\mathrm{max}}/2\pi=23$ MHz at
$T_{\mathrm{tot}}=0.594\;\mu$s as in [24] and then scale as
$\Omega_{\mathrm{max}},\;\Delta_{\mathrm{max}}\propto 1/T_{\mathrm{tot}}$ in
order to maintain adiabatic evolution. The eCD frequency is chosen to be
$\omega=300$ MHz. (e) QPT of the proposed eCD method: difference between the
eCD $\chi$-matrix, $\chi_{\mathrm{eCD}}$, and the ideal one. The parameters
are the same as in (c) and (d) with fixed $V/2\pi=500$ MHz and
$T_{\text{tot}}=0.594\;\mu$s. The colorbar displays the phase of the matrix
elements. The labels $i,x,y,z$ refer to the indices of the operator basis,
while the odd (even) axes ticks refer to the upper (lower) row of
$\chi$-matrix labels.

We see from the structure of the effective Hamiltonian in the first two lines
of Eq. (12) that the required pulse is not separable with respect to the two
atoms. This can be overcome by approximating $f_{0}\approx f_{1}$, what would
lead to disentangled counterdiabatic pulses that are much easier to implement
in practice. In this case the eCD Hamiltonian can be realized as

$\displaystyle H_{\mathrm{eCD}}^{\prime}(t)=$
$\displaystyle\sqrt{\omega\bar{f}}\big{[}\sin(\omega
t)(\ket{r}\\!\bra{1}\otimes\mathds{1}+\mathds{1}\otimes\ket{r}\\!\bra{1})+\mathrm{H.c.}$
$\displaystyle+\cos(\omega t)(P_{r}\otimes\mathds{1}+\mathds{1}\otimes
P_{r})\big{]},$ (15)

where we take $\bar{f}$ to be the arithmetic mean between the two original
pulses $f_{0}$ and $f_{1}$. In contrast to the eCD protocol discussed above,
this now approximated protocol has to be combined with the original $H(t)$
from Eq. (1) in order to guarantee the following of the desired state
evolution.

Results corresponding to Fig. 2(c,d) are shown in Fig. 3 for initial states
involving either $\ket{01}$ or the Rydberg interaction sensitive state
$\ket{11}$. The structure along the $x$-axis in the fidelity in (a) and (c)
arises from population fluctuations produced by the micromotion of the
oscillating Hamiltonian. Similar variations on a smaller scale can already be
seen in Fig. 2(d). The results of the superposition state
$(\ket{00}+\ket{11})/\sqrt{2}$ (b) are obviously less good than in Fig. 2.
There is a parameter range, see the maximal fidelity area in Fig. 3(b), in
which the achieved fidelities are nevertheless high. This specific regime was
analysed further and we can understand the effect by separating the error
introduced by the infidelity of state $\ket{11}$, see Fig. 3(c), and by its
unwanted phase error
$\phi_{e}=\pi-|\mathrm{arg}\braket{11}{\psi_{\mathrm{fin}}}|$, see Fig. 3(d).
The fidelity of the superposition state $(\ket{00}+\ket{11})/\sqrt{2}$ in (b)
is determined mainly by the phase. This explains the overall lower values and
in particular the bright area of higher fidelity. The approximated
$H_{\mathrm{eCD}}$ generates an additional phase which has opposite sign
compared to the original phase error (8) and decreases for small Rydberg
blockades, therefore almost cancelling the total phase error along the bright
area. This interesting interference phenomenon allows us to find such regions
of high-fidelity.

Figure 3: Infidelity for the separable accelerated protocol for initial states
$(\ket{00}+\ket{01})/\sqrt{2}$ (a), $(\ket{00}+\ket{11})/\sqrt{2}$ (b) and
$\ket{11}$ (c) and relative phase error
$|\phi_{e}|/\pi=|1-|\mathrm{arg}\braket{11}{\psi_{\mathrm{fin}}}|/\pi|$
accumulated for initial state $\ket{11}$ (d). The arrows indicate, that the
infidelities of the superposition state can be decomposed into a population
error, see (c) and a relative phase error of the sensitive state $\ket{11}$,
see (d). We use $H(t)$ from (1) with $\Omega_{\mathrm{max}}/2\pi=17$ MHz and
$\Delta_{\mathrm{max}}/2\pi=23$ MHz, in combination with the approximated eCD
method with $\omega=1$ GHz.

Having characterized the very promising performance of the eCD gate at the
level of a single population transfer, we next analyze the impact of these
fidelities in an application to a practical proof-of-concept quantum
algorithm.

##  IV Application to a QEC code

Figure 4: Realization of a minimal quantum-error correction algorithm based on
the eCD CZ gate with pulses $\Omega_{\mathrm{max}}/2\pi=17$ MHz,
$\Delta_{\mathrm{max}}/2\pi=23$ MHz, protocol time
$T_{\mathrm{tot}}=0.594\;\mu$s, blockade strength $V/2\pi=500$ MHz and eCD
frequency $\omega=350$ MHz. The quantum circuit is represented in (h), while
in (a)-(c) and (d)-(e) the difference of the matrix elements of the recovered
density matrix for the corrected qubit is shown. The first and second row
corresponds to the case in which no flip or one flip have occurred,
respectively, for three different initial states, $\ket{0}$ [(a) and (d)],
$\ket{1}$ [(b) and (e)] and $(\ket{0}+\ket{1})/\sqrt{2}$ [(c) and (f)]. The
colorscale indicates the phase of the residual matrix elements. The infidelity
of the different configurations is shown in (g) for the circuit based on the
eCD gate (blue) compared to the same results using the adiabatic CZ gate in
[24] (red). We excluded the initial state $\ket{11}$ without a bitflip, as
this configuration is dark to the CZ gate. The dashed lines show the average
over the five states. In all cases the minimum fidelities achieved by the eCD
gate are larger than $0.999$.

As an exemplary application of the modified controlled-phase gate, we study a
minimal implementation of a quantum-error correcting scheme. We consider the
simplified three-qubit code circuit depicted in Fig. 3(h) [69], which can
correct a potential flip occurring to the first physical qubit only. It
involves a sequence of CNOT gates with two ancillary qubits which encodes an
input state $\ket{\psi_{\mathrm{in}}}=\alpha\ket{0}+\beta\ket{1}$, prepared in
the first qubit, into a superposition
$\ket{\psi_{L}}=\alpha\ket{0_{L}}+\beta\ket{1_{L}}$ of the logical code words
$\ket{0_{L}}=\ket{000}$ and $\ket{1_{L}}=\ket{111}$. Error correction is
implemented via syndrome measurements on the ancillary qubit, which is
preceded by additional CNOT gates, and a conditional recovery [69, 65]. In the
Rydberg realization considered here, the CNOT gates are implemented via
single-qubit Hadamard gates and the modified CZ gate, according to [65]

$\text{CNOT}=R_{c}^{x}(\pi)R_{t}^{y}(-\pi/2)\text{CZ}R_{t}^{y}(\pi/2)R_{c}^{x}(\pi),$
(16)

where the labels $c$ and $t$ indicate the control and target qubit,
respectively, and $R_{c/t}^{i}(\theta)=\exp(i\theta\sigma_{c/t}^{i}/2)$ is a
single-qubit rotation. To focus on the impact of CZ gate errors, other
operations such as measurement and recovery are assumed to be perfect. We
study the performance of the code by inspecting the difference between the
output density matrix post error correction and the input density matrix. This
is shown in Fig. 4, where the recovery in the case of no error occurrences
(first row) is compared with the case of a bit-flip occurrence (lower row),
for three different initial states $\ket{0}$ [(a) and (d)], $\ket{1}$ [(b) and
(e)] and $(\ket{0}+\ket{1})/\sqrt{2}$ [(c) and (f)]. For all input states, all
matrix elements of the reconstructed density matrix match the input $\rho$ up
to errors below $10^{-2}$. In Fig. 4(g), we further compare the fidelity
results obtained by the eCD gate with those obtained by the adiabatic gate
proposed in [24]. For almost all error configurations, the eCD method
outperforms the adiabatic approach, and ensures very high fidelity of
recovery, with similar values for all configurations.

These results thus confirm that the modified gate enables the implementation
of high-fidelity quantum operations, not creating errors that can become
unexpectedly amplified during the algorithm. For the case of no error, best
recovery is obtained for an initial state $\ket{1}$, which is because the
state $\ket{1}$ is dark to the CZ gate after application of the single-qubit
Hadamard gate. For the error case, the best recovery corresponds to the input
$\ket{0}$ state.

##  V Conclusion

We have proposed a shortcut-to-adiabatic controlled-phase gate for Rydberg
atoms based on effective counterdiabatic fields. The starting protocol from
[24] has the advantage that it avoids populating additional intermediate
states similarly to STIRAP-like protocols [70, 46]. Our proposed accelerated
version is less prone to decoherence occurring on longer time scales, while it
remains intrinsically stable to its adiabatic backbone. The enhanced
performance of the gate is verified in numerical simulations both in single-
state transfer fidelities and through full quantum process tomography,
confirming the effectiveness of the gate for any choice of input state.
Rydberg-interaction induced phase errors are irrelevant if only the
counterdiabatic field $H_{\mathrm{CD}}$ is implemented, making the protocol
interesting for a wide parameter range, without the need of additional phase
compensation [63]. We finally exemplify the application of the eCD gate to a
practical quantum computation, by discussing its use in a proof-of-concept
quantum error correction algorithm and demonstrating typically a large
fidelity gain as compared to the purely adiabatic approach. Future work will
be dedicated to disentangle the effective counterdiabatic pulses with respect
to the two atoms not only approximately as shown above, possibly following
ideas proposed in [23, 14].

###### Acknowledgements.

S.W. acknowledges funding by the National Recovery and Resilience Plan (NRRP),
Mission 4 Component 2 Investment 1.3 – Call for tender No. 341 of 15/03/2022
of Italian Ministry of University and Research funded by the European Union –
NextGenerationEU, Project number PE0000023, Concession Decree No. 1564 of
11/10/2022 adopted by the Italian Ministry of University and Research, CUP
D93C22000940001, Project title “National Quantum Science and Technology
Institute” (NQSTI).

##  Appendix A Adiabatic elimination

In the limit of large Rydberg interaction $V$, the state $\ket{rr}$ is barely
populated during the evolution and can thus be adiabatically eliminated as
detailed in the following. Then only the states $\ket{rr}$, $\ket{d_{+}}$ and
$\ket{11}$ are coupled to each other, and the time-dependent Schrödinger
equation thus gives the following system of equations in the manifold spanned
by such states,

$\displaystyle i\frac{\partial_{t}c_{rr}(t)}{dt}$
$\displaystyle=[\Delta(t)+V]c_{rr}(t)+\frac{\Omega(t)}{\sqrt{2}}c_{+}(t),$
(17a) $\displaystyle i\frac{\partial_{t}c_{+}(t)}{dt}$
$\displaystyle=\frac{\Omega(t)}{\sqrt{2}}[c_{rr}(t)+c_{11}(t)]+\frac{\Delta}{2}c_{+}(t),$
(17b) $\displaystyle i\frac{\partial_{t}c_{11}(t)}{dt}$
$\displaystyle=\frac{\Omega(t)}{\sqrt{2}}c_{+}(t),$ (17c)

for the corresponding coefficients $c_{rr}(t)$, $c_{+}(t)$ and $c_{11}(t)$ of
the wavefunction. Assuming that $V$ is sufficiently strong for the population
of level $\ket{rr}$ not to vary during the dynamics, we impose
$\partial_{t}c_{rr}(t)=0$. This gives

$c_{rr}(t)=-\frac{\Omega(t)}{\sqrt{2}[\Delta(t)+V]}c_{+}(t).$ (18)

Inserting this result in the remaining Eqs. (17b) and (17c), and keeping
leading-order terms in $V^{-1}$, one obtains the two-state Hamiltonian for
$\ket{11}$ and $\ket{d_{+}}$

$H=\frac{1}{2}\begin{bmatrix}0&\sqrt{2}\Omega(t)\\\
\sqrt{2}\Omega(t)&2\Delta(t),\end{bmatrix}$ (19)

as reported in Eq. (9).

##  Appendix B Derivation of the eCD field

The idea of the eCD scheme is to construct a Hamiltonian $H_{\mathrm{eCD}}(t)$
which contains fast oscillations in controllable parameters, whose ‘averaged’
effect is the realization of $H_{\mathrm{CD}}(t)$. Namely, one forces the
propagator produced by $H_{\mathrm{eCD}}(t)$ to approximately match the
propagator generated by $H_{\mathrm{CD}}(t)$. This procedure can be understood
as Floquet engineering [54, 52, 53] of a time-dependent Hamiltonian, and is
typically done by computing first terms of a high-frequency Magnus expansion
[33]. We show now that the proposed eCD field of Eq. (12) generates the
counterdiabatic dynamics given by the Hamiltoninan of Eqs. (10) up to order
$1/\omega$.

Starting from Eq. (12), the first terms $H_{\mathrm{eff}}^{(n)}$ of the
efffective Hamiltonian computed via the Magnus expansion [71] generated by
$H_{\mathrm{eCD}}(t)$ at the end of the oscillation period $T=2\pi/\omega$
read

$\displaystyle H_{\mathrm{eff}}^{(0)}=$
$\displaystyle\int_{0}^{T}H_{\mathrm{eCD}}(t)dt=0,$ (20) $\displaystyle
H_{\mathrm{eff}}^{(1)}=-$
$\displaystyle\frac{1}{2}\int_{0}^{T}dt_{1}\int_{0}^{t_{1}}dt_{2}[H_{\mathrm{eCD}}(t_{1}),H_{\mathrm{eCD}}(t_{2})]$
(21) $\displaystyle=$
$\displaystyle\frac{f_{0}(t)}{2}\big{[}(-i\ket{r}\\!\bra{1}+i\ket{1}\\!\bra{r})\otimes
P_{0}$ (22)
$\displaystyle\qquad+P_{0}\otimes(-i\ket{r}\\!\bra{1}+i\ket{1}\\!\bra{r})\big{]}$
(23)
$\displaystyle+\frac{f_{1}(t)}{2}\big{[}-i\ket{d_{+}}\\!\bra{11}+i\ket{11}\\!\bra{d_{+}}\big{]}.$
(24)

The oscillatory $H_{\mathrm{eCD}}(t)$ thus produces an effective end-of-period
Hamiltonian with the same desired structure of $H_{\mathrm{CD}}(t)$.

## References

  * Gallagher [2005] T. Gallagher,  _Rydberg Atoms_, Cambridge Monographs on Atomic, Molecular and Chemical Physics (Cambridge University Press, 2005).
  * Rempe  _et al._ [1987] G. Rempe, H. Walther, and N. Klein, Observation of quantum collapse and revival in a one-atom maser, Phys. Rev. Lett. 58, 353 (1987).
  * Varcoe  _et al._ [2000] B. T. H. Varcoe, S. Brattke, M. Weidinger, and H. Walther, Preparing pure photon number states of the radiation field, Nature 403, 743 (2000).
  * Raimond  _et al._ [2001] J. M. Raimond, M. Brune, and S. Haroche, Manipulating quantum entanglement with atoms and photons in a cavity, Rev. Mod. Phys. 73, 565 (2001).
  * Haroche [2013] S. Haroche, Nobel lecture: Controlling photons in a box and exploring the quantum to classical boundary, Rev. Mod. Phys. 85, 1083 (2013).
  * Graham  _et al._ [2019] T. M. Graham, M. Kwon, B. Grinkemeyer, Z. Marra, X. Jiang, M. T. Lichtman, Y. Sun, M. Ebert, and M. Saffman, Rydberg-mediated entanglement in a two-dimensional neutral atom qubit array, Phys. Rev. Lett. 123, 230501 (2019).
  * Levine  _et al._ [2019] H. Levine, A. Keesling, G. Semeghini, A. Omran, T. T. Wang, S. Ebadi, H. Bernien, M. Greiner, V. Vuletić, H. Pichler, and M. D. Lukin, Parallel implementation of high-fidelity multiqubit gates with neutral atoms, Phys. Rev. Lett. 123, 170503 (2019).
  * Scholl  _et al._ [2021] P. Scholl, M. Schuler, H. J. Williams, A. A. Eberharter, D. Barredo, K.-N. Schymik, V. Lienhard, L.-P. Henry, T. C. Lang, T. Lahaye, A. M. Läuchli, and A. Browaeys, Quantum simulation of 2D antiferromagnets with hundreds of Rydberg atoms, Nature 595, 233 (2021).
  * Ebadi  _et al._ [2021] S. Ebadi, T. T. Wang, H. Levine, A. Keesling, G. Semeghini, A. Omran, D. Bluvstein, R. Samajdar, H. Pichler, W. W. Ho, S. Choi, S. Sachdev, M. Greiner, V. Vuletić, and M. D. Lukin, Quantum phases of matter on a 256-atom programmable quantum simulator, Nature 595, 227 (2021).
  * Morgado and Whitlock [2021] M. Morgado and S. Whitlock, Quantum simulation and computing with Rydberg-interacting qubits, AVS Quantum Science 3, 10.1116/5.0036562 (2021), 023501.
  * Adams  _et al._ [2019] C. S. Adams, J. D. Pritchard, and J. P. Shaffer, Rydberg atom quantum technologies, Journal of Physics B: Atomic, Molecular and Optical Physics 53, 012002 (2019).
  * Saffman  _et al._ [2010] M. Saffman, T. G. Walker, and K. Mølmer, Quantum information with Rydberg atoms, Rev. Mod. Phys. 82, 2313 (2010).
  * Saffman [2016] M. Saffman, Quantum computing with atomic qubits and rydberg interactions: progress and challenges, Journal of Physics B: Atomic, Molecular and Optical Physics 49, 202001 (2016).
  * Su  _et al._ [2016] S.-L. Su, E. Liang, S. Zhang, J.-J. Wen, L.-L. Sun, Z. Jin, and A.-D. Zhu, One-step implementation of the Rydberg-Rydberg-interaction gate, Phys. Rev. A 93, 012306 (2016).
  * Pelegrì  _et al._ [2022] G. Pelegrì, A. J. Daley, and J. D. Pritchard, High-fidelity multiqubit Rydberg gates via two-photon adiabatic rapid passage, Quantum Science and Technology 7, 045020 (2022).
  * Pagano  _et al._ [2022] A. Pagano, S. Weber, D. Jaschke, T. Pfau, F. Meinert, S. Montangero, and H. P. Büchler, Error budgeting for a controlled-phase gate with strontium-88 Rydberg atoms, Phys. Rev. Res. 4, 033019 (2022).
  * Jandura and Pupillo [2022] S. Jandura and G. Pupillo, Time-Optimal Two- and Three-Qubit Gates for Rydberg Atoms, Quantum 6, 712 (2022).
  * Robicheaux  _et al._ [2021] F. Robicheaux, T. M. Graham, and M. Saffman, Photon-recoil and laser-focusing limits to Rydberg gate fidelity, Phys. Rev. A 103, 022424 (2021).
  * Chang  _et al._ [2023] T. Chang, T. N. Wang, H. H. Jen, and Y.-C. Chen, High-fidelity Rydberg control-$Z$ gates with time-optimal pulses, arXiv:2303.16395 (2023).
  * Brierley  _et al._ [2012] R. T. Brierley, C. Creatore, P. B. Littlewood, and P. R. Eastham, Adiabatic state preparation of interacting two-level systems, Phys. Rev. Lett. 109, 043002 (2012).
  * Vitanov  _et al._ [2017] N. V. Vitanov, A. A. Rangelov, B. W. Shore, and K. Bergmann, Stimulated Raman adiabatic passage in physics, chemistry, and beyond, Rev. Mod. Phys. 89, 015006 (2017).
  * Barends  _et al._ [2016] R. Barends, A. Shabani, L. Lamata, J. Kelly, A. Mezzacapo, U. L. Heras, R. Babbush, A. G. Fowler, B. Campbell, Y. Chen, Z. Chen, B. Chiaro, A. Dunsworth, E. Jeffrey, E. Lucero, A. Megrant, J. Y. Mutus, M. Neeley, C. Neill, P. J. J. O’Malley, C. Quintana, P. Roushan, D. Sank, A. Vainsencher, J. Wenner, T. C. White, E. Solano, H. Neven, and J. M. Martinis, Digitized adiabatic quantum computing with a superconducting circuit, Nature 534, 222 (2016).
  * Han  _et al._ [2016] R. Han, H. K. Ng, and B.-G. Englert, Implementing a neutral-atom controlled-phase gate with a single Rydberg pulse, Europhysics Letters 113, 40001 (2016).
  * Saffman  _et al._ [2020] M. Saffman, I. I. Beterov, A. Dalal, E. J. Páez, and B. C. Sanders, Symmetric Rydberg controlled-$Z$ gates with adiabatic pulses, Phys. Rev. A 101, 062309 (2020).
  * Wang  _et al._ [2019] T. Wang, Z. Zhang, L. Xiang, Z. Jia, P. Duan, Z. Zong, Z. Sun, Z. Dong, J. Wu, Y. Yin, and G. Guo, Experimental Realization of a Fast Controlled-Z Gate via a Shortcut to Adiabaticity, Phys. Rev. Appl. 11, 034030 (2019).
  * Torrontegui  _et al._ [2013] E. Torrontegui, S. Ibáñez, S. Martínez-Garaot, M. Modugno, A. del Campo, D. Guéry-Odelin, A. Ruschhaupt, X. Chen, and J. G. Muga, Shortcuts to adiabaticity, in  _Advances in Atomic, Molecular, and Optical Physics_, Advances In Atomic, Molecular, and Optical Physics, Vol. 62, edited by E. Arimondo, P. R. Berman, and C. C. Lin (Academic Press, 2013) pp. 117–169.
  * Guéry-Odelin  _et al._ [2019] D. Guéry-Odelin, A. Ruschhaupt, A. Kiely, E. Torrontegui, S. Martínez-Garaot, and J. G. Muga, Shortcuts to adiabaticity: Concepts, methods, and applications, Rev. Mod. Phys. 91, 045001 (2019).
  * Bason  _et al._ [2012] M. G. Bason, M. Viteau, N. Malossi, P. Huillery, E. Arimondo, D. Ciampini, R. Fazio, V. Giovannetti, R. Mannella, and O. Morsch, High-fidelity quantum driving, Nat. Phys. 8, 147?152 (2012).
  * An  _et al._ [2016] S. An, D. Lv, A. del Campo, and K. Kim, Shortcuts to adiabaticity by counterdiabatic driving for trapped-ion displacement in phase space, Nature Communications 7, 12999 (2016).
  * Zhou  _et al._ [2017] B. Zhou, A. Baksic, H. Ribeiro, C. Yale, F. Heremans, P. Jerger, A. Auer, G. Burkard, A. Clerk, and D. Awschalom, Accelerated quantum control using superadiabatic dynamics in a solid-state lambda system, Nature Physics 13, 330 (2017).
  * Boyers  _et al._ [2019] E. Boyers, M. Pandey, D. K. Campbell, A. Polkovnikov, D. Sels, and A. O. Sushkov, Floquet-engineered quantum state manipulation in a noisy qubit, Phys. Rev. A 100, 012341 (2019).
  * Vepsäläinen  _et al._ [2019] A. Vepsäläinen, S. Danilin, and G. S. Paraoanu, Superadiabatic population transfer in a three-level superconducting circuit, Science Advances 5, eaau5999 (2019).
  * Petiziol  _et al._ [2018] F. Petiziol, B. Dive, F. Mintert, and S. Wimberger, Fast adiabatic evolution by oscillating initial hamiltonians, Phys. Rev. A 98, 043436 (2018).
  * Claeys  _et al._ [2019] P. W. Claeys, M. Pandey, D. Sels, and A. Polkovnikov, Floquet-engineering counterdiabatic protocols in quantum many-body systems, Phys. Rev. Lett. 123, 090602 (2019).
  * Dorner  _et al._ [2005] U. Dorner, T. Calarco, P. Zoller, A. Browaeys, and P. Grangier, Quantum logic via optimal control in holographic dipole traps, Journal of Optics B: Quantum and Semiclassical Optics 7, S341 (2005).
  * Treutlein  _et al._ [2006] P. Treutlein, T. W. Hänsch, J. Reichel, A. Negretti, M. A. Cirone, and T. Calarco, Microwave potentials and optimal control for robust quantum gates on an atom chip, Phys. Rev. A 74, 022312 (2006).
  * Cozzini  _et al._ [2006] M. Cozzini, T. Calarco, A. Recati, and P. Zoller, Fast rydberg gates without dipole blockade via quantum control, Optics Communications 264, 375 (2006), quantum Control of Light and Matter.
  * Goerz  _et al._ [2011] M. H. Goerz, T. Calarco, and C. P. Koch, The quantum speed limit of optimal controlled phasegates for trapped neutral atoms, Journal of Physics B: Atomic, Molecular and Optical Physics 44, 154011 (2011).
  * 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, D. Sugny, and F. K. Wilhelm, Training Schrödinger’s cat: quantum optimal control, Eur. Phys. J. D 69, 279 (2015).
  * Martinis and Geller [2014] J. M. Martinis and M. R. Geller, Fast adiabatic qubit gates using only ${\sigma}_{z}$ control, Phys. Rev. A 90, 022307 (2014).
  * Chen  _et al._ [2015] Y.-H. Chen, Y. Xia, Q.-Q. Chen, and J. Song, Fast and noise-resistant implementation of quantum phase gates and creation of quantum entangled states, Phys. Rev. A 91, 012325 (2015).
  * Santos and Sarandy [2015] A. C. Santos and M. S. Sarandy, Superadiabatic Controlled Evolutions and Universal Quantum Computation, Scientific Reports 5, 15775 (2015).
  * Petiziol  _et al._ [2019] F. Petiziol, B. Dive, S. Carretta, R. Mannella, F. Mintert, and S. Wimberger, Accelerating adiabatic protocols for entangling two qubits in circuit qed, Phys. Rev. A 99, 042315 (2019).
  * Shen  _et al._ [2019] C.-P. Shen, J.-L. Wu, S.-L. Su, and E. Liang, Construction of robust Rydberg controlled-phase gates, Opt. Lett. 44, 2036 (2019).
  * Yan  _et al._ [2019] T. Yan, B.-J. Liu, K. Xu, C. Song, S. Liu, Z. Zhang, H. Deng, Z. Yan, H. Rong, K. Huang, M.-H. Yung, Y. Chen, and D. Yu, Experimental realization of nonadiabatic shortcut to non-abelian geometric gates, Phys. Rev. Lett. 122, 080501 (2019).
  * Petiziol  _et al._ [2020] F. Petiziol, E. Arimondo, L. Giannelli, F. Mintert, and S. Wimberger, Optimized three-level quantum transfers based on frequency-modulated optical excitations, Sci. Rep. 10, 2185 (2020).
  * Wang  _et al._ [2018] T. Wang, Z. Zhang, L. Xiang, Z. Jia, P. Duan, W. Cai, Z. Gong, Z. Zong, M. Wu, J. Wu, L. Sun, Y. Yin, and G. Guo, The experimental realization of high-fidelity ?shortcut-to-adiabaticity? quantum gates in a superconducting xmon qubit, New J. Phys. 20, 065003 (2018).
  * Kolodrubetz  _et al._ [2017] M. Kolodrubetz, D. Sels, P. Mehta, and A. Polkovnikov, Geometry and non-adiabatic response in quantum and classical systems, Physics Reports 697, 1 (2017).
  * Hegade  _et al._ [2021] N. N. Hegade, K. Paul, Y. Ding, M. Sanz, F. Albarrán-Arriagada, E. Solano, and X. Chen, Shortcuts to adiabaticity in digitized adiabatic quantum computing, Phys. Rev. Appl. 15, 024038 (2021).
  * Wurtz and Love [2022] J. Wurtz and P. J. Love, Counterdiabaticity and the quantum approximate optimization algorithm, Quantum 6, 635 (2022).
  * Morris [1964] R. J. Morris, Theory of adiabatic rapid passage for three equally spaced levels, Phys. Rev. 133, A740 (1964).
  * Eckardt [2017] A. Eckardt, Colloquium: Atomic quantum gases in periodically driven optical lattices, Rev. Mod. Phys. 89, 011004 (2017).
  * Goldman and Dalibard [2014] N. Goldman and J. Dalibard, Periodically driven quantum systems: Effective hamiltonians and engineered gauge fields, Phys. Rev. X 4, 031027 (2014).
  * Bukov  _et al._ [2015] M. Bukov, L. D’Alessio, and A. Polkovnikov, Universal high-frequency behavior of periodically driven systems: from dynamical stabilization to Floquet engineering, Adv. Phys. 64, 139 (2015).
  * Berry [2009] M. V. Berry, Transitionless quantum driving, Journal of Physics A: Mathematical and Theoretical 42, 365303 (2009).
  * Demirplak and Rice [2003] M. Demirplak and S. A. Rice, Adiabatic population transfer with control fields, The Journal of Physical Chemistry A 107, 9937 (2003).
  * Demirplak and Rice [2008] M. Demirplak and S. A. Rice, On the consistency, extremal, and global properties of counterdiabatic fields, The Journal of Chemical Physics 129, 154111 (2008).
  * Ibáñez  _et al._ [2012] S. Ibáñez, X. Chen, E. Torrontegui, J. G. Muga, and A. Ruschhaupt, Multiple Schrödinger Pictures and Dynamics in Shortcuts to Adiabaticity, Phys. Rev. Lett. 109, 100403 (2012).
  * Ibáñez  _et al._ [2013] S. Ibáñez, X. Chen, and J. G. Muga, Improving shortcuts to adiabaticity by iterative interaction pictures, Phys. Rev. A 87, 043402 (2013).
  * Martínez-Garaot  _et al._ [2014] S. Martínez-Garaot, E. Torrontegui, X. Chen, and J. G. Muga, Shortcuts to adiabaticity in three-level systems using lie transforms, Phys. Rev. A 89, 053408 (2014).
  * Deffner  _et al._ [2014] S. Deffner, C. Jarzynski, and A. del Campo, Classical and quantum shortcuts to adiabaticity for scale-invariant driving, Phys. Rev. X 4, 021013 (2014).
  * Ehret [2022] T. Ehret, Optimizing elementary quantum gates, Bachelor Thesis, Heidelberg University (2022).
  * Maller  _et al._ [2015] K. M. Maller, M. T. Lichtman, T. Xia, Y. Sun, M. J. Piotrowicz, A. W. Carr, L. Isenhower, and M. Saffman, Rydberg-blockade controlled-not gate and entanglement in a two-dimensional array of neutral-atom qubits, Phys. Rev. A 92, 022336 (2015).
  * Deffner and Campbell [2017] S. Deffner and S. Campbell, Quantum speed limits: from heisenberg?s uncertainty principle to optimal quantum control, J. Phys. A: Math. and Theor. 50, 453001 (2017).
  * Nielsen and Chuang [2010] M. A. Nielsen and I. L. Chuang,  _Quantum Computation and Quantum Information: 10th Anniversary Edition_ (Cambridge University Press, Cambridge, 2010).
  * Mohseni  _et al._ [2008] M. Mohseni, A. T. Rezakhani, and D. A. Lidar, Quantum-process tomography: Resource analysis of different strategies, Phys. Rev. A 77, 032322 (2008).
  * Johansson  _et al._ [2012] J. Johansson, P. Nation, and F. Nori, Qutip: An open-source python framework for the dynamics of open quantum systems, Computer Physics Communications 183, 1760 (2012).
  * Johansson  _et al._ [2013] J. Johansson, P. Nation, and F. Nori, Qutip 2: A python framework for the dynamics of open quantum systems, Computer Physics Communications 184, 1234 (2013).
  * Foot [2004] C. J. Foot,  _Atomic Physics_ (Oxford University Press, Oxford, 2004).
  * Møller  _et al._ [2008] D. Møller, L. B. Madsen, and K. Mølmer, Quantum Gates and Multiparticle Entanglement by Rydberg Excitation Blockade and Adiabatic Passage, Phys. Rev. Lett. 100, 170504 (2008).
  * Blanes  _et al._ [2009] S. Blanes, F. Casas, J. Oteo, and J. Ros, The Magnus expansion and some of its applications, Physics Reports 470, 151 (2009).


# On dynamics of gasless combustion in slowly varying periodic media: periodic
fronts, their stability and propagation-extinction-diffusion-reignition
pattern.

Amanda Matson Department of Mathematical Sciences, Kent State University,
Kent, OH 44242, USA<EMAIL_ADDRESS>Leonid Kagan School of
Mathematical Sciences, Tel Aviv University, Tel Aviv 69978, Israel. E-mail:
<EMAIL_ADDRESS>Claude-Michel Brauner Institut de Mathématiques de
Bordeaux UMR CNRS 5251, Université de Bordeaux, 33405 Talence, France. E-mail:
<EMAIL_ADDRESS>Gregory Sivashinsky School of Mathematical
Sciences, Tel Aviv University, Tel Aviv 69978, Israel.
<EMAIL_ADDRESS>Peter V. Gordon Department of Mathematical
Sciences, Kent State University, Kent, OH 44242, USA. E-mail:
<EMAIL_ADDRESS>
###### Abstract

In this paper we consider a classical model of gasless combustion in a one
dimensional formulation under the assumption of ignition temperature kinetics.
We study the propagation of flame fronts in this model when the initial
distribution of the solid fuel is a spatially periodic function that varies on
a large scale. It is shown that in certain parametric regimes the model
supports periodic traveling fronts. An accurate asymptotic formula for the
velocity of the flame front is derived and studied. The stability of periodic
fronts is also explored, and a critical condition in terms of parameters of
the problem is derived. It is also shown that the instability of periodic
fronts, in a certain parametric regimes, results in a propagation-extinction-
diffusion-reignition pattern which is studied numerically.

Keywords: Gasless combustion, periodic traveling fronts, instability,
propagation-extinction-diffusion-reignition pattern.

##  1 Introduction

The present study is concerned with combustion waves in gasless systems in
which the final products, as well as initial reactants, are in the solid
phase. The premixed reactants are typically powdered mixtures of metals with a
nonmetal oxide. When ignited, these systems support a self-propagating
deflagration. The combustion process results in a solid product commonly
referred to as self-propagating high temperature synthesis (SHS) due to its
application to the high-temperature synthesis of new refracting materials
[Gorosh2022, Mer2003, Mer1997]. A substantial part of theoretical and
experimental studies of gasless combustion are targeted to understand the
dynamics of flame fronts propagating toward unburned spatially homogeneous
solid fuel. It is well established by now that such fronts commonly suffer
from various instabilities. Apparently the first experimental study reporting
unstable oscillatory flames in solid samples is [bel50] and dates back to
1950’s. Later experimental studies revealed more complicated regimes of
propagation such as spinning fronts [mer73], other regimes were discovered in
[ex1, ex2]. The theoretical exploration of pulsating fronts in a one
dimensional configuration was performed in [matsiv78]. The spinning regime in
cylindrical samples was studied in [siv81], spiral waves in planer geometry
were discussed in [volmat]. Later theoretical studies (see e.g [BM90]) showed
existence of non-trivial chaotic patterns of propagations. Reviews of some of
these results can be found in [Marg1991, V3]. Systematic description of
experimentally observed modes of propagations is given in [Ivl3]. Earlier
studies of instabilities were based on a standard single temperature constant
density approximation adiabatic model with a first order reaction rate and
Arrhenius temperature kinetics. Some additional effects such as the impact of
temperature kinetics and the reaction order [vol2, gol], the presence of a
heat loss [kur] and other factors were later explored. Let us note that the
majority of these results were obtained by formal linear stability analysis of
flame fronts, via numerical simulations or combination of thereof. Rigorous
analysis of nonlinear stability of gasless flame fronts poses a very non-
trivial mathematical challenge as the absence of material diffusivity leads to
the analysis of highly degenerate systems. Apparently the most advanced
results in this direction, so far, were obtained in [ghaz10].

While the understanding of propagation regimes of flame fronts moving towards
homogeneous unburned solid reactant is of major importance, it is also of
interest to consider a situation in which the initial profile of an unreacted
solid varies, particularly in a periodic manner. Indeed, in real systems it is
quite common to have a certain level of inhomogeneity in solid reactants
which, in extreme case of solids consisting of separated solid grains, leads
to emergence of discrete reaction waves [Muk2008]. Unlike homogenous case, the
theoretical explorations of flame fronts in the solids of varying
concentration are more limited. It is apparent, however, that the presence of
variations in solid fuel concentration results in new effects such as, for
example, frequency locking [mdk, mdk2] and hump effect [CMB2].

In this paper we are interested in the analysis of flame fronts propagating
toward unburned solid fuel whose concentration is varying on a large scale
(scales much larger than the characteristic width of the flame front). We show
that the periodicity in fuel concentration leads to the emergence of periodic
flame fronts. We give a comprehensive asymptotic description of such fronts
and study their stability. We show that the instability of periodic fronts, in
certain parametric regimes, leads to a very interesting propagation-
extinction-diffusion-reignition pattern.

The paper is organized as follows. In section 2 we give a mathematical
formulation of the problem. In section 3 we derive an asymptotic expression
for the velocity of periodic traveling fronts and compare results with
numerical simulations. In section 4 we discuss stability of periodic traveling
fronts. In section 5 we discuss a propagation-extinction-diffusion-reignition
regime which appears as a result of flame front instability.

##  2 Formulation

The current study is based on the classical model of gasless combustion. This
model, in one dimensional formulation, reads:

$\displaystyle\left\\{\begin{array}[]{lll}\Theta_{t}=\Theta_{xx}+W(\Theta,\Phi),\\\
\Phi_{t}=-W(\Theta,\Phi),\end{array}\right.$ (2.3)

where $\Theta$ and $\Phi$ are appropriately scaled temperature and
concentration of deficient solid reactant, $x\in\mathbb{R},$ $t>0$ are
spatiotemporal coordinates and $W(\Theta,\Phi)$ is the reaction rate. In this
work we will assume that the reaction rate has a first order reaction and
ignition type temperature kinetics. Namely,

$\displaystyle
W(\Theta,\Phi):=\left\\{\begin{array}[]{ll}A\Phi&\Theta\geq\Theta_{i},\\\
0&\mbox{otherwise},\end{array}\right.$ (2.6)

where $\Theta_{i}\in(0,1)$ is an ignition temperature and

$\displaystyle A=\frac{\Theta_{i}}{1-\Theta_{i}}$ (2.7)

is a scaling factor. The scaling factor is chosen in such a way that the
velocity of the flame front in a spatially homogenous concentration field is
unity. A comprehensive study of such flame fronts and their stability in the
presence of material diffusion is given in [BGKS, CMB1]. A fully non-linear
reduced model for diffusive instabilities of such fronts was derived and
analyzed in [dima].

The goal of this paper is to analyze flame propagation regimes for model
(2.3)-(2.6) with front like initial conditions for the temperature and
spatially periodic slowly varying conditions for the concentration of the
deficient solid reactant. Specifically, we study problem (2.3)-(2.6)
complemented with boundary like conditions for the temperature $\Theta$ in far
fields

$\displaystyle\Theta(x,t)\to 0,\quad\mbox{as}\quad x\to\infty,$
$\displaystyle\Theta_{x}(x,t)\to 0,\quad\mbox{as}\quad x\to-\infty,$ (2.8)

and front like initial conditions:

$\displaystyle\Theta(x,t=0)=\Theta_{0}(x),$ (2.9)

where $\Theta_{0}$ is a smooth monotone decreasing function that approaches
zero and one as $x\to\pm\infty$ sufficiently fast. In addition we assume that
$\Theta_{0}(0)=\Theta_{i}.$ The initial concentration of solid deficient
reactant is prescribed as follows:

$\displaystyle\Phi(x,t=0)=\Phi_{0}(x)$ (2.10)

We set the following assumption on the initial concentration field $\Phi_{0}$.
The function $\Phi_{0}$ is an increasing function on $x\in(-\infty,0)$ that
approaches zero as $x\to-\infty$ sufficiently fast and $\Phi_{0}(0)=1$. For
$x\in(0,\infty)$ we set $\Phi_{0}(x)=\Psi_{\varepsilon}(x)=\Psi(\varepsilon
x)$. Here and below $0<\Psi<1$ is a prescribed continuous $2\pi$ periodic
function on $[0,\infty)$ satisfying $\Psi(0)=1$ and $0<\varepsilon\ll 1$. That
is $\Psi_{\varepsilon}$ is a slowly varying periodic function with
characteristic scale $1/\varepsilon$. In numerical simulations and formal
asymptotic constructions presented in this paper we set:

$\displaystyle\Psi(x)=1-a(1-\cos(x)),$ (2.11)

with $0<a<1/2$. Representative initial conditions for temperature and
concentration of the deficient reactant are depicted in Figure 1.

Figure 1: Profile of initial conditions for the temperature $\Theta_{0}$ (red)
and concentration of the deficient solid reactant $\Phi_{0}$ (blue) with
$\Theta_{i}=1/2$, $a=0.25$ and $\varepsilon=0.1$

To characterize the propagation phenomena for problem (2.3)-(2.10) it is
convenient to introduce the notion of an ignition interface which is a
leftmost position in the flame front such that the maximal temperature to the
right of this position is below the ignition temperature. Namely, we define
the position of the ignition interface $q(t)$ as follows:

$\displaystyle q(t)=\min\\{x~{}\mbox{such
that}\quad\Theta(y,t)<\Theta_{i}\quad\mbox{for}\quad y>x\\}.$ (2.12)

The velocity of the ignition interface is defined as the derivative of this
position, that is $dq(t)/dt$. It is easy to verify that when crossing the
ignition interface the temperature and its spatial first derivative remain
continuous:

$\displaystyle\Theta(q(t),t)=\Theta_{i},\quad[\Theta(q(t),t)]=0,\quad[\Theta_{x}(q(t),t)]=0,$
(2.13)

where $[\cdot]$ stands for the jump of the quantity. We note that $q(t)$ is a
well defined quantity even in the absence of the reaction.

As a first step of our study we show (on the formal level) that problem
(2.3)-(2.10) admits periodic traveling fronts solutions. We construct these
solutions and derive a formula for the speed of propagation in the next
section.

##  3  Periodic traveling fronts.

In this section we construct an approximation of periodic fronts that emerge
as a long term limit of solutions for problem (2.3)-(2.10) in certain
parametric regimes.

It is well known that the dynamics of flame fronts for model (2.3)-(2.10) are
predominately defined by their behavior in the vicinity of the ignition
temperature rather than far fields. In view of this observation and the fact
that the fuel concentration field varies only on a large scale
$1/\varepsilon,$ one may expect that in the first approximation the dynamics
of flame fronts can be captured by considering the problem (2.3)-(2.10) on the
intermediate scale $1/\sqrt{\varepsilon}$. On spacial scale
$1/\sqrt{\varepsilon}$ the initial fuel concentration field can be viewed as a
locally prescribed constant. Hence, introducing a traveling front ansatz
$(\Theta,\Phi)(x,t)=(\theta,\phi)(\xi)$ with
$\xi=x-{\varepsilon^{-1}}{q(\tau)},$ where $\tau=\varepsilon t$ and $q(\tau)$
are slow time and the position of the ignition interface respectively, and
substituting this ansatz into (2.3) for sufficiently small $\varepsilon$ one
ends up with the following auxiliary problem:

$\displaystyle\left\\{\begin{array}[]{lll}\theta^{\prime\prime}+V\theta^{\prime}=0,&V\phi^{\prime}=0,&\xi>0\\\
\theta^{\prime\prime}+V\theta^{\prime}+A\phi=0,&V\phi^{\prime}-A\phi=0,&\xi<0,\end{array}\right.$
(3.3)

where $V=dq(\tau)/d\tau$ is the velocity of the ignition interface and prime
stands for the derivative with respect to $\xi$. This problem is complemented
with the continuity of the solution when crossing the ignition interface which
is set to be located at $\xi=0$

$\displaystyle\theta=\Theta_{i},\quad[\theta]=[\theta^{\prime}]=[\phi]=0,\quad\mbox{at}\quad\xi=0,$
(3.4)

and far field conditions

$\displaystyle\left\\{\begin{array}[]{lll}\theta\to
0,&\phi\to\Sigma&\mbox{as}\quad\xi\to\infty,\\\ \theta^{\prime}\to 0,&\phi\to
0,&\mbox{as}\quad\xi\to-\infty,\end{array}\right.$ (3.7)

where $\Sigma$ is the local concentration of solid fuel.

One can easily verify that for arbitrary $0<\Theta_{i}<\Sigma\leq 1$ this
problem admits a unique solution that reads:

$\displaystyle\theta(\xi)=\left\\{\begin{array}[]{ll}\Theta_{i}\exp\left(-\sqrt{\frac{\Sigma-\Theta_{i}}{1-\Theta_{i}}}\xi\right),&\xi>0,\\\
\Sigma-\left(\Sigma-\Theta_{i}\right)\exp\left(\frac{\Theta_{i}}{\sqrt{(1-\Theta_{i})(\Sigma-\Theta_{i})}}\xi\right),&\xi\leq
0,\end{array}\right.$ (3.10)
$\displaystyle\phi(\xi)=\left\\{\begin{array}[]{ll}\Sigma,&\xi>0,\\\
\Sigma\exp\left(\frac{\Theta_{i}}{\sqrt{(1-\Theta_{i})(\Sigma-\Theta_{i})}}\xi\right),&\xi\leq
0,\end{array}\right.$ (3.13) $\displaystyle
V=V(\Sigma,\Theta_{i})=\sqrt{\frac{\Sigma-\Theta_{i}}{1-\Theta_{i}}}.$ (3.14)

A typical profile of this traveling front is depicted in Figure 2.

Figure 2: Traveling front profile given by (3.10), (3.13) with
$\Theta_{i}=1/2$ and $\Sigma=3/4$.

In the vicinity of characteristic size $1/\sqrt{\varepsilon}$ of the ignition
interface located at the point $x=\varepsilon^{-1}q(\tau)$ the shape of the
flame front is expected to be close to the one given by $\eqref{eq:fr4}$ and
the position of this interface is defined by the following first order initial
value problem

$\displaystyle\frac{d}{d\tau}q(\tau)=\sqrt{\frac{\Psi(q(\tau))-\Theta_{i}}{1-\Theta_{i}}},\quad
t>0,\quad q(0)=0.$ (3.15)

For any given function $\Psi$ such that $\Sigma_{min}:=\min_{s\geq
0}\Psi(s)\geq\Theta_{i}$ the solution of problem (3.14) can be obtained (at
least in the implicit form) by separation of variables. In particular, when
$\Psi$ is given by (2.11) we have the following, closed form, expressions for
the position of the ignition interface

$\displaystyle q(\tau)=2{\rm am}\left(\frac{\tau}{2},\lambda\right),$ (3.16)

and

$\displaystyle\frac{d}{d\tau}q(\tau)={\rm
dn}\left(\frac{\tau}{2},\lambda\right),$ (3.17)

for its velocity. Where,

$\displaystyle\lambda:=\frac{1-\Sigma_{min}}{1-\Theta_{i}}=\frac{2a}{1-\Theta_{i}},$
(3.18)

and ${\rm am}\left(\frac{\tau}{2},\lambda\right),{\rm
dn}\left(\frac{\tau}{2},\lambda\right),$ are Jacobi amplitude and Jacobi
elliptic functions respectively [AS].

Figure 3: The position $q(\tau)$ and velocity $dq(\tau)/d\tau$ of the ignition
interface as a function of slow time $\tau$ and parameter $\lambda$.

Hence the dynamics of the ignition front is fully determined by the specific
value of the parameter $\lambda\in[0,1)$. For fixed ignition temperature
$\Theta_{i}$ this parameter characterizes the amplitude of oscillations in the
initial data for the deficient reactant and increases from zero to one as the
amplitude of oscillations increases. We note that $\lambda\to 0$ corresponds
to the case when $\Sigma_{min}=1$ (that is the absence of the oscillations in
the initial data for the deficient reactant) in which case the ignition front
propagates with the constant velocity one. The case $\lambda\to 1$ corresponds
to a situation when oscillations in the initial value of the deficient
reactant are maximal. In the latter case the front is stagnant, and its
velocity becomes zero when the front reaches the position with minimal
concentration of solid fuel. The ignition front does not exist when
$\lambda>1$. The profiles for the position and velocity of the ignition
interface as a function of slow time $\tau$ and parameter $\lambda$ are
depicted in Figure 3.

Let us also note that the function $dq(\tau)/d\tau$ is a periodic function
with period

$\displaystyle T\left(\lambda\right)=4F\left(\frac{\pi}{2},\lambda\right),$
(3.19)

where $F\left(\frac{\pi}{2},\lambda\right)$ is the elliptic integral of the
first kind [AS]. It is easy to check that the period of oscillation for the
periodic front is an increasing function of the parameter $\lambda$. Moreover,
one has $T\left(\lambda\right)\to 2\pi$ when $\lambda\to 0$ and
$T\left(\lambda\right)\to\infty$ when $\lambda\to 1$. The dependency of the
period on $\lambda$ is depicted in Figure 4.

Figure 4: Period of oscillation of the velocity profile as a function of the
parameter $\lambda$.

The comparison of direct numerical simulation of problem (2.3)-(2.10) with the
asymptotic solution constructed in this section shows excellent agreement in
terms of position, velocity and local structure of the ignition front provided
$\Sigma_{min}$ is sufficiently larger than the ignition temperature. Figure 5
shows several snapshots of the propagating traveling front over one period of
oscillation as a solution of problem (2.3)–(2.10) with $\Theta_{i}=0.5,$
$\Sigma_{min}=0.6~{}(a=0.2)$ and $\varepsilon=0.05.$

Figure 5: Snapshots of periodic traveling front profiles for problem
(2.3)–(2.10) with $\Theta_{i}=1/2,$ $\Sigma_{min}=0.6~{}(a=0.2)$ and
$\varepsilon=0.05$. The red line represents the temperature profile and the
blue line represents the profile of the deficient reactant. Purple dashed line
indicates value of the ignition temperature.

The velocity of the ignition front obtained in this simulation is very close
to the one obtained from the asymptotic expression (the maximal relative error
is about $5\%$), see Figure 6. Moreover, as expected the local structure of
the traveling front is quite close to the one given by (3.10), (3.13) in a
$1/\sqrt{\varepsilon}$ vicinity of the ignition interface, see Figure 7. As
clearly seen from Figure 5, outside of the $1/\sqrt{\varepsilon}$ vicinity of
the ignition interface, the periodic front exhibits an oscillatory pattern in
burned fuel that decays due to thermal diffusion. This oscillatory tail does
not significantly influence the propagation of the ignition interface.

Figure 6: Velocity of the ignition interface for problem (2.3)–(2.10) with
$\Theta_{i}=1/2,$ $\Sigma_{min}=0.6~{}(a=0.2)$ and $\varepsilon=0.05$ (blue
solid line) and its asymptotic value obtained from expression (3.17) (black
dashed line). Figure 7: A snapshot of the traveling front profile for problem
(2.3)–(2.10) with $\Theta_{i}=1/2,$ $\Sigma_{min}=0.6~{}(a=0.2)$ and
$\varepsilon=0.05$ in $1/\sqrt{\varepsilon}$ vicinity of the ignition
interface at the instance when the front is located at the local minimum of
the initial concentration profile. Red and blue solid lines indicate the
profile of the temperature and concentration of the deficient reactant. Black
dashed line represent these quantities obtained by the asymptotic formulas
(3.10) and (3.13).

For smaller values of the parameter $a$ (that is for larger $\Sigma_{min}$)
the accuracy of the asymptotic predictions is even better. Figure 8 shows the
velocity of propagation of the flame front obtained by means of direct
numerical simulation of problem (2.3)–(2.10) and from asymptotic expression
(3.17) with $\Theta_{i}=0.5,$ $\Sigma_{min}=0.7~{}(a=0.15)$ and
$\varepsilon=0.05$. In this case, the relative error of asymptotic prediction
is below $1.5\%.$ Similar behavior is also seen for other values of the
ignition temperature.

Figure 8: Velocity of the ignition interface for problem (2.3)–(2.10) with
$\Theta_{i}=1/2,$ $\Sigma_{min}=0.7~{}(a=0.15)$ and $\varepsilon=0.05$ (blue
solid line) and its asymptotic value obtained from expression (3.17) (black
dashed line).

However, when the value $\Sigma_{min}$ gets closer to $\Theta_{i}$ the picture
changes substantially and the behavior of flames fronts become very different
from predictions of asymptotic formulas. This discrepancy, to no surprise, is
due to the instability of the periodic front that occurs when $\Sigma_{min}$
gets closer to $\Theta_{i}$. As will be shown in the next section, the
stability of the periodic flame front is controlled by a parameter

$\displaystyle\Lambda=\frac{\Theta_{i}}{\Sigma_{min}}.$ (3.20)

When the value of this parameter is below some critical value
$\Lambda^{*}\approx 6/7\approx 0.857$ the periodic front appears to be stable
but exhibits instabilities when $\Lambda$ exceeds $\Lambda^{*}$. Note that the
velocity of the periodic flame fronts depicted in Figure 8 correspond to value
of $\Lambda=0.714$ and hence represents a stable flame front well within the
stability region. Velocity profile depicted in Figure 6 corresponds to
$\Lambda=0.833$ that is slightly below the stability threshold $\Lambda^{*}.$

The behavior of flame fronts above the stability threshold depends quite
sensitively on specific values of controlling parameters $\Theta_{i}$ and
$\Sigma_{min}$. We will discuss this behavior in Section 5. In the next
section we will discuss the linear stability of the periodic flame fronts.

##  4  Linear stability of fronts

As evident from the discussion in the preceding sections, the local behavior
of periodic flame fronts is determined (in the first approximation) by
traveling front solutions of auxiliary problem (3.3)–(3.7). It is then
expected that the stability of the periodic fronts can be associated with the
stability of solutions of this problem given by (3.10)–(3.14) with
$\Sigma=\Sigma_{min}$. In what follows we briefly discuss the stability of
these solutions as solutions of problem (2.3).

Setting

$\displaystyle\Theta(x,t)=\theta(\eta)+\exp(\omega t)\vartheta(\eta),$ (4.1)
$\displaystyle\Phi(x,t)=\phi(\eta)+\exp(\omega t)\varphi(\eta),$ (4.2)
$\displaystyle\eta=x-Vt-\exp(\omega t)v,$ (4.3)

substituting this ansatz into (2.3) and assuming that $\vartheta,\varphi$ and
$v$ are infinitesimally small after some tedious, but straightforward
computations similar to these in [BGKS] one ends up with the following
dispersion repletion,

$\displaystyle{\cal H}(\Theta_{i},\Sigma,\omega)=0,\quad\omega\in\mathbb{C}.$
(4.4)

where

$\displaystyle{\cal
H}(\Theta_{i},\Sigma,\omega)=\left[\Theta_{i}-2V(1-\Theta_{i})p\right]\times$
$\displaystyle\left[A+\left(\frac{A+\omega}{V}\right)^{2}\right]+\frac{A\Sigma}{V}\left[p-\frac{A+\omega}{V}\right],$
(4.5)

and

$\displaystyle p=\frac{1}{2}\left[-V+\sqrt{V^{2}+4\omega}\right].$ (4.6)

The neutral stability condition is then obtained as a solution of (4.4)
restricted to $\omega$ with zero real part that is:

$\displaystyle{\cal H}(\Theta_{i},\Sigma,i\Omega)=0,\quad\Omega\in\mathbb{R}.$
(4.7)

The neutral stability curve $\Theta^{*}(\Sigma)$ obtained as a solution of
equation (4.7) turned out to be a straight line with the accuracy above
$0.5\%$ given as

$\displaystyle\Theta^{*}(\Sigma)\approx\Theta^{*}_{i}(\Sigma)\approx\Theta_{i}^{*}(1)\Sigma\approx
0.857\Sigma,\quad\Sigma\in(0,1].$ (4.8)

We also note that $\Theta^{*}(1)=6/7$ as follows from [CMB1, Section 4]. It is
also easy to check that for $\Theta_{i}^{*}(\Sigma)<\Sigma\leq\Theta_{i}$ the
traveling fronts are unstable. For $0<\Theta_{i}\leq\Theta_{i}^{*}(\Sigma)$
the traveling fronts are linearly stable and for $\Theta_{i}>\Sigma,$ thanks
to (3.15), the fronts do not exist, see Figure 9.

Figure 9:  Neutral stability curve and regions of stability, instability and
non existence of fronts.

The velocity of the front along the neutral stability curve and frequency
along the neutral stability curve, obtained as a solution of equation (4.7),
are shown on Figures 10 and 11.

Figure 10:  Velocity (blue line) along the neutral stability curve (red line).
Figure 11:  Frequency (blue line) along neutral stability curve (red line).

Thus, in full agreement with the numerical simulation, we conclude that the
instability of periodic fronts takes place when $\Sigma_{min}$ is below the
stability threshold. One may expect that, similar to the case of homogeneous
solid fuel concentration, there is a rich family of instability patterns which
get more and more complex as $\Sigma_{min}$ decreases. In the following
section we will discuss one of (apparently many possible) regimes of
propagation: propagation - extinction- diffusion - reignition regime. This
regime occurs exclusively due to the presence of inhomogeneity in solid
reactant and ,to the best of our knowledge, was not reported in the
literature.

##  5 Propagation - extinction - diffusion - reignition regime

In this section we will discuss the behavior of the combustion fronts above
the stability threshold. This behavior is rather sensitive to the specific
values of the controlling parameters $\Theta_{i}$ and $\Sigma_{min}$. When
$\Theta_{i}$ is relatively large $\Theta_{i}>\Theta_{i}^{*}\approx 0.37$ the
crossing of the stability threshold leads to extinction.

Figure 12:  The velocity of the ignition interface obtained form solution of
problem (2.3)–(2.10) with $\Theta_{i}=0.3,$ $\Sigma_{min}=0.32~{}(a=0.34)$ and
$\varepsilon=0.05$.

However, when $\Theta_{i}$ is smaller than $\Theta^{*}$ and $\Lambda$ is
slightly above the stability threshold $\Lambda^{*},$ an interesting
propagation - extinction - diffusion - reignition regime emerges.

Figure 13: Representative snapshots of flame front profile (rad) and
concentration profile (blue) for problem (2.3)–(2.10) with $\Theta_{i}=0.3$
and $\Sigma_{min}=0.32~{}(a=0.34)$ and $\varepsilon=0.05$.

In this regime comparable parts of the initial concentration curve are above
and below the stability threshold. The dynamics of this regime can be
described as follows. Initially the ignition front propagates as a stable
periodic front in the region when concentration of the deficient reactant is
high. As the front propagates, the concentration decreases and the front
enters in the fuel lean instability region which eventually leads to complete
extinction. Theses processes proceed at a relatively fast reactive time scale.
After the extinction takes place, the only mechanism that impacts the flame
front is the diffusion of temperature. This slow process puts dynamics on a
characteristic diffusive time scale. The flame front starts to diffuse and
widens. This results in a slow temperature increase in the region where
concentration of the solid reactant is higher. This dynamic finally results in
re-ignition of the reactant that triggers a local explosion which reignites
the flame. The flame then starts to propagate up to the point when the low
concentration region is reached. This process repeats itself in an
approximately periodic manner.

Figure 12 shows the velocity of the ignition interface when $\Theta_{i}=0.3$
and $\Sigma_{min}=0.32~{}(a=0.34)$. Figure 13 depicts snapshots of dynamic of
the flame front in this regime. The further increase of oscillations in the
initial concentration leads to extinction. For $\Theta_{i}=0.3,$ for example,
when $a$ exceeds value $0.413$ ($\Sigma_{min}=0.174$) the extinction takes
place.

Let us also note that, as evident from Figure 12, the flame front in
propagation - extinction - diffusion - reignition regime exhibits propagation
with negative velocity. This type of behavior was reported for a model of
nearly equi-diffusive combustion [add1], but known to be absent in gasless
combustion of homogeneous samples [add2].

Acknowledgments. The work of AM, LK, GS and PVG was supported, in part, by the
US-Israel Bination Science Foundation (Grant 2020-005).

Declaration of Competing Interest. Authors declare no competing interests.

## References