For computing the hash tables, we can apply the construction of Alon and Naor
[3]. There, they provide a deterministic hash table of $t$ elements into
$O(t)$ space with read access of $O(1)$ time. While they did not discuss the
parallel implementation of their construction, their algorithm can be
implemented in $\mathrm{poly}(\log n)$ depth in a straightforward way. Their
approach is derandomizing a randomized hash function using the method of
conditional expectation on epsilon-biased spaces. They define a potential
function (see section 3.1. of [3]) which is a simple aggregation and can be
parallelized. We do not discuss the full details as the implementation is
straightforward. ∎

## Acknowledgments

M.G., C.G., S.I., and V.R. were supported in part by the European Research
Council (ERC) under the European Unions Horizon 2020 research and innovation
program (grant agreement No. 853109) and the Swiss National Science Foundation
(project grant 200021_184735). B.H. was supported in part by NSF grants
CCF-1814603, CCF-1910588, NSF CAREER award CCF-1750808, a Sloan Research
Fellowship, funding from the European Research Council (ERC) under the
European Union’s Horizon 2020 research and innovation program (grant agreement
949272), and the Swiss National Science Foundation (project grant
200021_184735).

## References

  * ABCP [96] B. Awerbuch, B. Berger, L. Cowen, and D. Peleg.  Fast network decompositions and covers.  J. of Parallel and Distributed Computing, 39(2):105–114, 1996\.
  * AGLP [89] B. Awerbuch, A. V. Goldberg, M. Luby, and S. A. Plotkin.  Network decomposition and locality in distributed computation.  In Proc. 30th IEEE Symp. on Foundations of Computer Science (FOCS), pages 364–369, 1989.
  * AN [96] Noga Alon and Moni Naor.  Derandomization, witnesses for boolean matrix multiplication and construction of perfect hash functions.  Algorithmica, 16(4):434–449, 1996.
  * AP [90] B. Awerbuch and D. Peleg.  Sparse partitions.  In Proc. 31st IEEE Symp. on Foundations of Computer Science (FOCS), pages 503–513, 1990.
  * BE [13] L. Barenboim and M. Elkin.  Distributed Graph Coloring: Fundamentals and Recent Developments.  Morgan & Claypool Publishers, 2013.
  * BEG+ [22] Marcel Bezdrighin, Michael Elkin, Mohsen Ghaffari, Christoph Grunau, Bernhard Haeupler, Saeed Ilchi, and Václav Rozhoň.  Deterministic distributed sparse and ultra-sparse spanners and connectivity certificates.  arXiv preprint arXiv:2204.14086, 2022.
  * BRS [89] Bonnie Berger, John Rompel, and Peter W Shor.  Efficient nc algorithms for set cover with applications to learning and geometry.  In Foundations of Computer Science, 1989., 30th Annual Symposium on, pages 54–59. IEEE, 1989.
  * BS [07] Surender Baswana and Sandeep Sen.  A simple and linear time randomized algorithm for computing sparse spanners in weighted graphs.  Random Structures & Algorithms, 30(4):532–563, 2007.
  * CG [21] Yi-Jun Chang and Mohsen Ghaffari.  Strong-diameter network decomposition.  In Proceedings of the 2021 ACM Symposium on Principles of Distributed Computing, pages 273–281, 2021.
  * EHRG [22] Michael Elkin, Bernhard Haeupler, Václav Rozhoň, and Christoph Grunau.  Deterministic low-diameter decompositions for weighted graphs and distributed and parallel applications.  In Proceedings of the Symposium on Foundations of Computer Science (FOCS), pages to appear, arXiv:2204.08254, 2022.
  * EN [16] M. Elkin and O. Neiman.  Distributed strong diameter network decomposition.  In Proc. 35th ACM Symp. on Principles of Distributed Computing (PODC), pages 211–216, 2016.
  * FGG+ [22] Salwa Faour, Mohsen Ghaffari, Christoph Grunau, Fabian Kuhn, and Václav Rozhoň.  Local distributed rounding: Generalized to mis, matching, set cover, and beyond.  In Manuscript, page to appear, 2022.
  * GGR [21] M. Ghaffari, C. Grunau, and V. Rozhon.  Improved deterministic network decomposition.  In Proc. 33rd ACM-SIAM Symp. on Discrete Algorithms (SODA), pages 2904–2923, 2021.
  * GHK [18] M. Ghaffari, D. Harris, and F. Kuhn.  On derandomizing local distributed algorithms.  In Proceedings of the Symposium on Foundations of Computer Science (FOCS), pages 662–673, 2018.
  * GK [18] Mohsen Ghaffari and Fabian Kuhn.  Derandomizing distributed algorithms with small messages: Spanners and dominating set.  In 32nd International Symposium on Distributed Computing (DISC 2018). Schloss Dagstuhl-Leibniz-Zentrum fuer Informatik, 2018.
  * GK [21] M. Ghaffari and F. Kuhn.  Deterministic distributed vertex coloring: Simpler, faster, and without network decomposition.  In Proc. 62nd IEEE Symp. on Foundations of Computer Science (FOCS), pages 1009–1020, 2021.
  * GKM [17] M. Ghaffari, F. Kuhn, and Y. Maus.  On the complexity of local distributed graph problems.  In Proc. 49th ACM Symp. on Theory of Computing (STOC), pages 784–797, 2017.
  * GY [20] Parikshit Gopalan and Amir Yehudayoff.  Concentration for limited independence via inequalities for the elementary symmetric polynomials.  Theory of Computing, 16(1):1–29, 2020.
  * Lin [92] N. Linial.  Locality in distributed graph algorithms.  SIAM Journal on Computing, 21(1):193–201, 1992.
  * Lub [86] M. Luby.  A simple parallel algorithm for the maximal independent set problem.  SIAM Journal on Computing, 15:1036–1053, 1986.
  * Lub [93] M. Luby.  Removing randomness in parallel computation without a processor penalty.  J. of Computer and System Sciences, 47(2):250–286, 1993.
  * MPX [13] Gary L Miller, Richard Peng, and Shen Chen Xu.  Parallel graph decompositions using random shifts.  In Proceedings of the twenty-fifth annual ACM symposium on Parallelism in algorithms and architectures, pages 196–203, 2013.
  * NS [95] M. Naor and L. Stockmeyer.  What can be computed locally?  SIAM Journal on Computing, 24(6):1259–1277, 1995.
  * Pel [00] D. Peleg.  Distributed Computing: A Locality-Sensitive Approach.  SIAM, 2000.
  * Pet [10] Seth Pettie.  Distributed algorithms for ultrasparse spanners and linear size skeletons.  Distributed Computing, 22(3):147–166, 2010.
  * PS [92] A. Panconesi and A. Srinivasan.  Improved distributed algorithms for coloring and network decomposition problems.  In Proc. 24th ACM Symp. on Theory of Computing (STOC), pages 581–592, 1992.
  * PY [18] Merav Parter and Eylon Yogev.  Congested clique algorithms for graph spanners.  In 32nd International Symposium on Distributed Computing, page 3, 2018.
  * RG [20] V. Rozhoň and M. Ghaffari.  Polylogarithmic-time deterministic network decomposition and distributed derandomization.  In Proc. 52nd ACM Symp. on Theory of Computing (STOC), pages 350–363, 2020.
  * RGH+ [22] Václav Rozhoň, Christoph Grunau, Bernhard Haeupler, Goran Zuzic, and Jason Li.  Undirected $1+\varepsilon$-shortest paths via minor-aggregates: Near-optimal deterministic parallel & distributed algorithms.  arXiv preprint arXiv:2204.05874, 2022.
  * RTZ [05] Liam Roditty, Mikkel Thorup, and Uri Zwick.  Deterministic constructions of approximate distance oracles and spanners.  In International Colloquium on Automata, Languages, and Programming, pages 261–272. Springer, 2005.
  * TZ [05] Mikkel Thorup and Uri Zwick.  Approximate distance oracles.  Journal of the ACM (JACM), 52(1):1–24, 2005.


# A single-species population in a refuge subject to heterogeneous diffusion

G. G. Piva1, C. Anteneodo1,2 1Department of Physics, Pontifical Catholic
University of Rio de Janeiro, PUC-Rio, and 2National Institute of Science and
Technology for Complex Systems, INCT-CS, Rua Marquês de São Vicente 225,
22451-900, Rio de Janeiro, RJ, Brazil

###### Abstract

We consider a modified nonlocal FKPP equation to describe the spatial
organization of a single-species population in a refuge (high-quality region
immersed in a hostile environment). In addition to the spatial-dependence of
the growth rate, which is the main distinguishing feature of a refuge, it is
realistic to consider that also the mobility can be heterogeneous. We focus on
two cases where the diffusivity responds to the population density, either
decreasing or increasing with population density. This state dependency can
produce a feedback that is missing in homogeneous environments. For
comparison, we also address space-dependent diffusivity with different values
inside and outside the refuge. The objective is to characterize how the
diffusion variability affects the spatial organization of the population
inside and outside the refuge. We note that the conditions for pattern
formation within the refuge are not strongly affected by the heterogeneity in
the diffusivity profile, but this affects the shape and stability of the
patterns.

##  I Introduction

A remarkable property of biological systems is the formation of spatial
structures. Patterns can emerge by self-organization, as a consequence of the
interactions between individuals, without external drives von2001diversity ;
rietkerk2008regular ; rietkerk2002self . To understand how self-organization
arises, simple models for the dynamics of the population density $u(x,t)$,
taking into account only elementary processes (such as reproduction,
competition and dispersal) have been very useful. In particular, the FKPP
equation murray2002 , which simply includes logistic growth plus normal
diffusion, is a classical model for the dynamics of a single-species
population. In the minimalistic framework of the FKPP equation and its
variants, spatially extended (nonlocal) competition is the main component to
produce self-organized spatial oscillations in the distribution of a
population murray2002 , whose characteristic wavelength is ruled by the range
of the intraspecies competitive interactions fuentes2003nonlocal ;
fuentes2004analytical ; hernandez2004clustering . As an illustrative example,
plant competition for water, known for generating spatial patterns, can be
considered a nonlocal process due to root spatial structure or water diffusive
dynamics escaff2015localized ; fernandez2019front . The nonlocality of other
elementary processes, such as reproduction and random dispersal, has a less
central role but can interfere constructively or destructively with the
possibility of pattern formation piva2021interplay .

Although patterns can emerge solely from interactions among individuals, they
are naturally affected by environmental conditions, as these control the rates
of the biological processes smith2003animal ; corley1977influence . Effects
can be observed at ecological scales in nature, as well as experimentally in
the laboratory. An example is the formation of regular vegetation patterns
tarnita2017theoretical  in a landscape constituted by patches characterized by
different growth rates pimm2019connecting ; maciel2013individual ;
wintle2019global . This scenario can also be produced artificially in the lab,
as when bacterial colonies are subjected to adverse conditions such as ultra-
violet light, except for a protected area (refuge) perry2005experimental .
Furthermore, growth-rate heterogeneity can induce pattern formation
dornelas2021landscape , in conditions under which it will not occur in
homogeneous media.

In previous works, normal diffusion with constant coefficient is typically
assumed, however, in real environments, non only growth rates but mobility can
also be heterogeneous, e.g., in active suspensions where movement is
influenced by nutrient gradients diluzio2005escherichia , or due to structural
features of the environment, as the movement of bacteria in porous media
sosa2017motility ; muskat1937flow . Furthermore, heterogeneous diffusion can
also emerge as a reaction arlt2019dynamics  that the individuals manifest, for
instance, in response to overcrowding or sparsity of a population, favoring,
or not, the random motion among other individuals colombo2012nonlinear ;
Cates2010 ; murray2002 ; lopez2006macroscopic ; kenkre2008nonlinearity ;
lenzi2001escape ; anteneodo2005non ; anteneodo2007brownian ; newman1980some ;
gurtin1977diffusion ; kareiva1983local ; birzu2019genetic . The specific way
organisms respond to concentration depends on several conditions and vary from
species to species gurtin1977diffusion ; kareiva1983local ; newman1980some ;
murray2002 . For instance, in populations of grasshoppers, the diffusion
coefficient is enhanced at high densities, where encounters between
individuals are more frequent, but in other species, this occurs at low
densities murray2002 . In any case, the state-dependency can feedback,
mitigating or reinforcing pattern formation. This scenario motivates us to
investigate the possible impact that, mainly, density-dependent diffusivity
might have on the spatial organization of a population in a refuge, whose
dynamics is described by a generalized nonlocal FKPP equation. This may bring
insights, for instance, on observations made in experiments with bacteria
perry2005experimental , where puzzling results cannot be explained by
considering only the simple form of the FKPP equation.

Aiming to investigate the effects of heterogeneous diffusion on pattern
formation, we focus on two functional forms of the state-dependent
diffusivity, namely, decay or increase with the population density, reflecting
enhanced mobility in response to sparseness and overcrowding, respectively.
For this purpose, exponential dependencies are studied as paradigm. For
comparison, uncoupling the diffusivity from the population density, we also
analyze the effect of the space-dependency where the diffusion coefficients
differ in the interior and exterior of the refuge.

##  II Model

We consider a single-species population living in a focal patch of size $L$
immersed in a large hostile environment in one dimension. In this scenario the
population is subjected to a positive growth rate $r_{in}>0$ inside the
refuge, and a negative one $r_{out}<0$ outside. Under these conditions, the
FKPP dynamics with spatially-extended competition defined by an interaction
kernel $\gamma$ is given by

$\partial_{t}u=\partial_{x}(D(u,x)\partial_{x}u)+r(x)u-u(\gamma\star u),$ (1)

where $r(x)$ is dornelas2021landscape

$r(x)=r_{in}+(r_{out}-r_{in})\Theta(|x|-L/2)\,,$ (2)

being $\Theta$ the Heaviside step function. The symbol “$\star$” stands for
the convolution operation that provides nonlocality through a kernel $\gamma$,
which for simplicity we consider to be a normalized rectangular shape of width
$2w$. This kind of system was considered before, for constant diffusion
coefficient $D$, that is, for homogeneous diffusivity dornelas2021landscape .
The extension we propose below, inspired in previous literature Cates2010 ;
tailleur2008statistical ; von2001diversity ; rietkerk2008regular ;
rietkerk2002self , assumes that the diffusivity can depend on the density
and/or on the spatial coordinate directly, $D(u,x)$, to reflect a reaction of
the mobility in response to the distribution of other individuals or to a
hostile medium.

### Space-dependent diffusivity

As a reference case, we will analyze in Sec. IV the impact on pattern
formation of a diffusivity profile $D(x)$ taking the values $D_{in}$ and
$D_{out}$ inside and outside the refuge, respectively, namely

$D(x)=D_{out}+(D_{in}-D_{out})\Theta_{s}[L/2-|x|],$ (3)

where $\Theta_{s}$ is a smoothed Heaviside step function, and the jump has
width $s$, such that the usual Heaviside is recovered in the limit $s\to 0$.
Similar settings have been used to study the role of space-dependent diffusion
on the critical patch size cantrell1999diffusion ; maciel2013individual ;
dos2020critical .

### State-dependent diffusivity

In the previous setting, mobility is determined by the degree of viability or
hostility of the environment, independently of the population level. But we
are mainly interested in variations of the diffusivity that are self-
generated, as response to the population level.

First in Sec. V, we will investigate a decreasing function of the population
density, namely,

$D_{1}(u)=d\exp(-u/\sigma)\,,$ (4)

where $d$ and $\sigma$ are positive parameters, and particularly $\sigma$
controls the decay with density, recovering a homogeneous diffusivity profile
in the limit $\sigma\to\infty$. This choice was motivated by previous works
which assume that density has a negative impact on diffusion Cates2010 ;
tailleur2008statistical . This functional form of $D_{1}(u)$ reflects a
reaction to sparsity, with greater mobility the more rarefied the population
is. In the case of a refuge within a hostile environment, it means a lower
diffusivity inside, where the population is more dense since $r_{out}<r_{in}$.

For the opposite possibility of enhanced response to overcrowding, we will use
as counterpart of Eq. (4),

$D_{2}(u)=d[(1-\exp(-u/\sigma)],$ (5)

for which homogeneity is obtained in the opposite limit $\sigma\to 0$.

##  III Methods

Results for the different scenarios described above, focusing on pattern
formation and mode stability, will be shown in the next sections. In all
cases, Eq. (1) was numerically integrated using a forward-time centered-space
algorithm (typically, $\Delta x=0.02$ and $\Delta t=10^{-5}$), with periodic
boundary conditions in a grid much larger than the refuge width, starting from
the initial condition corresponding to the homogeneous solution plus small
random fluctuations, namely, $u(x,t=0)\simeq u_{0}+\xi(x)$, being $\xi$ an
uncorrelated uniformly distributed variable of amplitude much smaller than
$u_{0}=r_{in}$. In numerical simulations, the refuge always spans the interval
$(-L/2,L/2)$, with fixed size $L=10$. The width of the rectangular kernel
$\gamma$ is also fixed ($w=1$). The numerical results are complemented by
analytical considerations based on linear stability analysis.

### Survival

The critical size of a refuge, for population survival, is known to depend on
growth and diffusion coefficients, both inside and outside the refuge, when
discontinuous binary forms of $r(x)$ and $D(x)$ are considered
maciel2013individual ; colombo2016intermittent ; colombo2018refuge ; ludwig79
; skellam91 ; kenkre03 . Succinctly, the critical size is derived assuming
small $u$, hence solving Eq. (1) in each domain after discarding
nonlinearities, and coupling the solutions through continuity conditions. As
the diffusivity is a continuous function in all the studied cases cases, then
we impose continuity for $u$ and its derivative at the interfaces. Hence, the
critical threshold for survival is approximately

$L_{c}\simeq
2\sqrt{\frac{D_{0}}{r_{in}}}\arctan\left(\sqrt{\frac{-r_{out}}{r_{in}}}\right)\,,$
(6)

where $D_{0}$ is the diffusivity at the interface for low density.

Since we are interested in studying the patterns that can appear in these
systems, we will work with values of the parameters far beyond the critical
values for survival, i.e., $L\gg L_{c}$, to warrant that the population does
not go extinct.

### Linear stability

For small perturbations around the homogeneous state $u_{0}$, we can
substitute $u(x,t)\simeq u_{0}+\epsilon(x,t)=u_{0}+\varepsilon\exp(ikx+\lambda
t)$, with $\varepsilon\ll u_{0}$, into Eq. (1), obtaining the linearized form

$\displaystyle\partial_{t}\epsilon$ $\displaystyle=$
$\displaystyle\bar{D}\partial_{xx}\epsilon-r\,\gamma\star\epsilon+\mathcal{O}(\epsilon^{2})\simeq\lambda(k)\epsilon,$
(7)

where $\bar{D}$ is the diffusion constant in the region of interest, that in
the state-dependent cases becomes $\bar{D}=D(u_{0})$. For the rectangular
kernel $\gamma(x)=\Theta(w-|x|)/(2w)$, Eq. (7) gives

$\lambda(k)=-\bar{D}\,k^{2}-r\frac{\sin{wk}}{wk}\,,$ (8)

which is the growth rate of mode $k$. Since at first order the diffusivity is
nearly constant, let us review the picture known for uniform diffusivity with
periodic boundary conditions. If $\lambda(k)>0$, the perturbation grows, with
characteristic mode $k^{*}$ which maximizes $\lambda(k)$, providing a
characteristic wavelength $\Lambda^{*}=2\pi/k^{*}$. Considering the extreme
value of the sinc function, one obtains the critical condition for pattern
formation kenkre03 ; colombo2012nonlinear

$\bar{D}\lesssim(1.43\pi)^{3}rw^{2}\simeq 0.011rw^{2}\equiv D_{c1}\,,$ (9)

for $r>0$; otherwise, no oscillations arise. This condition is crucial to
produce patterns in a homogeneous environment.

However, in a landscape with heterogeneous growth rate, such as in the case of
the refuge, and constant diffusivity $\bar{D}$, the uniform density state
$u_{0}$ can be destabilized even for $\lambda(k)<0$ and damped modes are
obtained from the complex roots of $\lambda(k)=0$ dornelas2021landscape . Let
us summarize the known results. The dominant mode (with less damping) is given
by the root $k=k_{R}\pm ik_{I}$ with the smallest imaginary part (in absolute
value). Its imaginary part gives the damping rate while its real part gives
the associated wavelength $\Lambda=2\pi/|k_{R}|$. Depending on the real and
imaginary parts of the relevant root, we can observe:

(I) sustained oscillations (when the imaginary part is zero), which
corresponds to condition (9);

(II) damped oscillations (when real and imaginary parts are non-null);

(III) absence of oscillations (when the real part is null).

These roots can be obtained numerically, but for the selected rectangular
kernel $\gamma$ it is possible to estimate the critical values analytically. A
mathematical estimate can be obtained from the Taylor expansion of
$\lambda(k)$ in Eq. (8), by finding the roots of the truncated series at fifth
order, which allows to obtain an explicit relation between the parameters for
the condition at which the real part of the relevant root becomes non null,
yielding $\bar{D}/(rw^{2})\lesssim 1/6+\sqrt{1/30}\simeq 0.35$, and with
higher order terms of the Taylor expansion, we arrive at

$\bar{D}\lesssim 0.37rw^{2}\equiv D_{c2}.$ (10)

This is the condition for the appearance of oscillations induced by growth-
rate heterogeneity, for constant diffusion coefficient $\bar{D}$. Inside the
refuge, we have $\bar{D}=D_{in}$ for the binary profile and $\bar{D}\simeq
D(u_{0})$ for the state-dependent ones.

##  IV Results for space-dependent diffusivity

Figure 1: (a) Stationary profiles of $u(x)$, $D(x)$ and $J$, under space-
dependent diffusivity in Eq. (3), for a refuge of size $L=10$ signaled by
vertical dotted lines, with different values of the interface width $s$
indicated in the legend, for $D_{in}=0.01$ in (a)-(c) and for $D_{in}=0.02$ in
(d)-(f). Fixed parameters are $D_{out}=1$, $r_{in}=1$, $r_{out}=-0.02$.

Let us consider the case where a distinct diffusion coefficient can occur
inside or outside the refuge , independently of the population level. Instead
of imposing continuity conditions maciel2013individual , we opted for a
smoothed form $\Theta_{s}(x)$ of the Heaviside function, in Eq. (3),
connecting the two levels by a narrow but finite interface of width $s\ll w$,
by defining
$\Theta_{s}(x)=\frac{\tanh\left(\frac{x-L/2}{s}\right)+\tanh\left(\frac{x+L/2}{s}\right)}{2\tanh(\frac{L}{2s})}$,
which tends to the Heaviside step function when $s\to 0$. In analogous way, we
might smooth the growth rate profile, but since its derivatives do not enter
the evolution equation, this smoothing is not significant in this case.

Typical stationary profiles of $u(x)$, $D(x)$ and the flux
$J=-D(x)\partial_{x}u$ are shown in Fig. 1, changing the interface width $s$,
for two different values of $D_{in}$. The width $s$ has an effect on the
spatial distribution of the population. When $s$ decreases, the flux at $\pm
L$ decreases in absolute value, allowing the growth of the peak closer to the
interface, while outside the refuge, the density becomes lower near the
interface. Then in the limit $s\to 0$, a density jump is expected, as well as
a jump in the first derivative of the density, such that the current
$J=-D(x)\partial_{x}u$ remains continuous at the interfaces. The fluxes are
maximal at the interfaces where they point outwards, always in the direction
of lower density. Inside the refuge, the fluxes tend to depopulate the peaks
but the logistic growth balances the fluxes and a steady state is attained.
Independently of the width $s$, for $D_{in}=0.01$ (a)-(c), Eq. (9)
approximately holds and oscillations are nearly sustained, except for the
external peak, while they are strongly damped for $D_{in}=0.02$ (d)-(f). Since
small values of $s$ produce qualitatively similar results, we fix $s=0.2$.

In Fig. 2 we show the effects of varying $D_{in}$ and $D_{out}$, as well as
$r_{in}$ and $r_{out}$. In (a) we observe the expressive reduction of damping
and increase of the oscillation amplitude associated to the increase of the
inner growth rate $r_{in}$. A similar effect is produced by decreasing
$D_{in}$, in (c), as predicted by the linear analysis. Furthermore, notice
that the profiles can become fragmented inside the refuge, due to regions of
nearly null density.

Figure 2: Stationary profiles of the population density in a refuge, for
$s=0.2$, and varying the parameters indicated in each panel, fixing (a)
$r_{out}=-0.02$, $D_{in}=0.01$, $D_{out}=1$, (b) $r_{in}=1$, $D_{in}=0.01$,
$D_{out}=1$, (c) $r_{in}=1$, $r_{out}=-0.02$, $D_{out}=1$, (d) $r_{in}=1$,
$r_{out}=-0.02$, $D_{in}=0.01$.

Values of $r_{out}<0$, varying in a wide range, have little effect on the
internal population, but the external population decays faster from the
interface with increasing $|r_{out}|$, without fluctuations, vanishing in the
limit of fully lethal conditions $r_{out}\to-\infty$. Increasing $D_{out}$
reduces the population in the refuge close to the interface, changing the
wavenumber of the oscillations, and makes the external population decay more
slowly with distance. In fact, in the external region, the population density
decays nearly exponentially from the interface, as
$u(x)\sim\exp\left(-\sqrt{\frac{r_{out}}{D_{out}}}[x-L/2]\right)$, as results
from the small $u$ approximate solution of Eq. (1). When $D_{out}<D_{in}$, the
population vanishes already near the interface.

We summarize the effects produced on pattern formation by the two more
relevant parameters $D_{in}$ and $r_{in}$, in the phase diagrams of Fig. 3.
The predictions of the linear analysis in Sec. III (for $\bar{D}=D_{in}$ and
$r=r_{in}$) are in accord with the kind of solutions (symbols) obtained from
numerical integration of Eq. (1). The dashed and full lines correspond to the
conditions given by Eqs. (9) and (10), respectively, which delimit the
different phases. The color plots represent the real (a) and imaginary (b)
parts of the dominant eigenvalue (minimum positive imaginary part, see Sec.
III), obtained numerically. Recall that, in (a) the white region (zero real
part) means no oscillations, in (b) the white region (zero imaginary part)
predicts non-damped oscillations, otherwise damped oscillations with multiple
peaks emerge within the refuge. Filled and open symbols correspond to the
classification of steady profiles obtained by numerical integration of the
FKPP Eq. (1), for which oscillations are observed, or not, respectively, in
good accordance to the regions predicted by the theoretical critical curves.

For the ranges of the parameters in this figure, a steady state is always
attained, and sustained oscillations are observed near the expected frontier
(dashed line),

Figure 3:  Phase diagrams in the plane $D_{in}-r_{in}$ for the binary profile
$D(x)$. The regions I-III correspond to those described in Sec. III. The heat
plots represent the real $k_{R}$ (a) and imaginary $k_{I}$ (b) parts of the
dominant root of $\lambda(k)$. The symbols correspond to results of the
numerical integration of Eq.-(1): oscillations are observed (filled), or not
(hollow). For the plus-circles a steady state with sustained oscillations is
reached. The dashed and solid lines correspond to the prediction given by Eqs.
(9) and (10). Fixed parameters are $D_{out}=1$, $r_{out}=-0.02$.

##  V Results for state-dependent diffusivity

In this section, we focus on diffusivity coefficients that are ruled by the
density. Since the diffusivity depends on the population density, it evolves
in time concomitantly with the density, in a self-consistent way.

### Enhanced reaction to sparsity: $D_{1}(u)=d\exp(-u/\sigma)$

Fig. 4 provides an illustrative example of the time evolution of the system
under adverse conditions except in the refuge, starting from a noisy uniform
density around $u_{0}$, until reaching a long-time profile (which occurs for
$t\gtrsim 100$ in the case of the figure).

Figure 4: Time evolution of population density (a), the density-dependent
diffusion coefficient $D_{1}(u)$ (b), and flux $J=-D_{1}(u)\partial_{x}u$ (c),
for times $t$ indicated in the legend, fixing $\sigma=0.3$, $d=1$, $r_{in}=1$,
for $r_{out}=-0.02$.

Initially (at $t=0$), the density, hence the diffusivity, are nearly constant
(with $D_{c1}<d\exp(u_{0}/\sigma)\simeq 0.036<D_{c2}$). In a first regime
($t\lesssim 1$), the two-level profile of the growth rate $r(x)$ induces
essentially a two-level population density, which in turn molds the
diffusivity profile. As soon as Eq. (10) holds, steady oscillations emerge,
with the characteristic wavelength predicted by the linear analysis, close to
$w/0.715$ colombo2012nonlinear .

In Fig. 4, we used $r_{out}=-0.02$, implying a weakly hostile environment. In
the more lethal case $r_{out}\to-\infty$, the external populations will go
extinct. Inside the refuge, the higher peaks of the density profile shift
towards the center when $|r_{out}|$ grows (not shown), but the qualitative
features remain essentially the same. Then we will choose $r_{out}=-0.02$ so
that the variations outside the refuge are amplified.

With regards to the parameters that define $D_{1}(u)$, while the pre-factor
$d$ regulates the maximum of diffusivity, $\sigma$ controls how sensitive to
density is the response of the diffusion coefficient. Let us analyze their
influence on the stationary profiles, looking at Fig. 5.

Figure 5: Stationary profiles $u$, and corresponding
$D_{1}(u)=d\exp(-u/\sigma)$, (a)-(b) with different values of $\sigma$
indicated in legend, keeping $d=1$; (c)-(d), with different values of $d$,
setting $\sigma=0.3$. Fixed parameters are $r_{in}=1$, $r_{out}=-0.02$.

Approximately, in the refuge $D_{in}\simeq d\exp(-u_{0}/\sigma)$ and outside
far from the interface $D_{out}\simeq d$. Therefore, $\sigma$ affects the
internal diffusivity, producing a level that can impact pattern formation, and
also the population near outside the refuge, affecting the persistence of the
population in this region. For increasing $\sigma$, when the internal
diffusivity level exceeds the critical value $D_{c2}\simeq 0.37$, oscillations
are spoiled, and the dependence on $u$ is flattened, producing homogeneous
diffusivity In the opposite case of small $\sigma$, the diffusivity is more
sensitive to the variations of density and results a density profile with
sharp peaks near the interfaces.

Moreover, from the viewpoint of the maintenance of the population level, a
larger $\sigma$, which tends to homogenize the diffusivity, reduces the
internal population, mainly the high crowding near the interface, but enhances
the outer population, which despite the negative growth rate is fed by the
outward fluxes.

The pre-factor $d$ affects both $D_{in}$ and $D_{out}$, as can be seen in Fig.
5(c)-(d), influencing the oscillations inside the refuge and the decay
outside. In this case, a variation of $d$ tending to homogenize the
diffusivity, which occurs at a low level, benefits the increase of the
internal population to the detriment of the external one.

The effects of both parameters of $D_{1}(u)$ on the oscillations inside the
refuge are summarized in Fig. 6, where we depict phase diagrams in the plane
$\sigma-d$, identifying the different structures that can emerge.

Figure 6: Phase diagrams in the plane $\sigma-d$ of $D_{1}(u)$, as in Fig. 3.
For the crossed circles a steady state is not observed. Fixed parameters are
$r_{in}=1$, $r_{out}=-0.02$.

The dashed and full lines correspond to the conditions given by Eqs. (9) and
(10), respectively, using $\bar{D}=D(u_{0})$, which delimit the different
phases I-III (see Sec. III). Despite the fluctuations of the density, the
approximation of a constant diffusivity $D(u_{0})$ in the internal region
gives a prediction in reasonable agreement with the results from simulations.
We verified that pattern wavelength agrees with the prediction
$\Lambda=2\pi/|k_{R}|$. Nonetheless, in the region where the imaginary part
vanishes and sustained oscillations are expected, what actually happens is
that a steady state is not attained (crossed symbols). High peaks give rise to
a locally low diffusion coefficient, then the fluxes outwards the peaks are
not able to balance the growth of the population and this effect is reinforced
the higher the peak, as can be seen in Fig. 7 where a region of the profiles
is amplified.

Notice that such feedback does not take place for density-independent
diffusivity, in which case a steady state was always achieved. In comparison
to the space-dependent case, we also highlight that for $D_{1}(u)$, when
stationary oscillations are formed in the refuge, they are strongly damped
towards the center. Moreover, the population outside the refuge can attain a
moderate level at the interface and decay slowly with the distance to the
interface.

Figure 7: Detail of the time evolution if $u(x)$ and $D_{1}(u)$ for
$\sigma=0.2$ and $d=1.6$ (corresponding to a crossed circle in Fig 6, and
$r_{in}=1$, $r_{out}=-0.02$.

### Enhanced reaction to overcrowding: $D_{2}(u)=d[1-\exp(-u/\sigma)]$

In this section we discuss the effects of the diffusivity $D_{2}(u)$ defined
in Eq. (5), which presents a higher level in response to overcrowding. For
large enough $d$, the high internal diffusivity level forbids pattern
formation, then we considered low values of $d$ only. Typical plots of
population density and diffusivity profiles are shown in Fig. 8.

Figure 8: Stationary profiles of $u$, $D(u)$ and flux $J$, generated by
$D(u)=D_{2}(u)=d[1-\exp(-u/\sigma)]$, for different values of the decay
parameter $\sigma$, with fixed $d=0.03$, $r_{in}=1$, $r_{out}=-0.02$.

Inside the refuge, the diffusivity is nearly constant, $d$, when $\sigma\ll
u_{0}$. As $\sigma$ increases, the average level of the population decreases,
favoring oscillations, which tend to become sustained when the diffusivity
profile becomes bounded by $D_{c1}$. Furthermore, due to the particular
dependency on $u$, the diffusivity is higher at the peaks of $u$ and very
small at the valleys, favoring a steady state with nearly sustained
oscillations, differently to what we observed in the opposite case of reaction
to sparsity. Moreover, this kind of oscillations can lead to fragmentation of
the population inside the refuge, for sufficiently large $\sigma$, a
phenomenon observed for the space-dependent profile but not for $D_{1}(u)$.

With regard to the distribution of the population outside the refuge, notice
the fast decay to a low level, even when $r_{out}$ is only slightly negative.
This is another difference with the profile $D_{1}(u)$.

Figure 9:  Phase diagrams in the plane $\sigma-d$ of $D_{2}(u)$, as in Fig. 6.
Patterns are formed except in the white region in (a). Fixed parameters
$r_{in}=1$, $r_{out}=-0.02$.

The effect of the parameters is summarized in the phase diagrams shown in Fig.
9. In general there is a good agreement between the classification of steady
states (symbols) and the predicted regions. However, where the linear
approximation predicts sustained oscillations (white region in panel (b)),
actually, a slight damping occurs, as can be seen in Fig. 8 for the case
$\sigma=10$ and $d=0.03$. The theoretical prediction based on the
approximation $u\sim u_{0}$ for the inner population can fail when large
amplitude oscillations develop.

##  VI Final remarks

We have discussed the impact of heterogeneous diffusion in a refuge immersed
in a hostile environment. We considered two forms of state-dependent
diffusivity $D(u)$ that mimic opposite reactions to the density, and also
space-dependent diffusivity $D(x)$ where the heterogeneity is anchored on the
quality of the environment, independently of the population.

In all the analyzed cases, when the average diffusivity level is sufficiently
low inside the refuge, for the other parameters fixed, internal patterns
emerge. Regardless of the details of the diffusivity profile inside the
refuge, its average level allows a good prediction of the frontiers for
patterning through the linear stability analysis. This indicates a certain
robustness of the process of pattern formation against heterogeneities in the
diffusivity profile.

Nevertheless, the type of heterogeneity has effects on the shape of the
patterns, not captured by the linear analysis. In the case of overcrowding
reaction, it is possible to have more uniform peaks and fragmentation can
occur, similar to what happens with the two-level profile. In the case of
mobility enhanced by sparsity through $D_{1}$, the peaks rise higher the
closer to the interface and a steady state may be unattainable due to
progressive condensation in these peaks feedback by the dependence on density.

Outside the refuge, the population decays from the interface more slowly with
higher levels of the external diffusivity and within the setting represented
by $D_{1}$, where the diffusivity is low for high concentrations. The
possibility of distinct diffusivity levels inside and outside could explain
profiles observed experimentally, which present tails from the interfaces with
moderate density of individuals even under very adverse conditions.

Our numerical results and theoretical considerations, exploring how density-
dependent diffusivity affects the organization of the population inside and
outside the refuge, reveal macroscopic signatures that could provide insights
about hidden mechanisms. Moreover, this knowledge can be used to refine
mathematical models. Within the particular context of the FKPP equation, in
previous attempts to explain the distribution of populations of bacteria in an
environment subject to UV light with a refuge perry2005experimental , the
diffusivity, measured for low bacterial density, is too high to allow
destabilization of the homogeneous state empirically observed. But if the
mobility is affected by density, being lower in the crowded refuge, or
distinct in hostile conditions, oscillations inside and tails outside the
refuge would be explained. In fact some of the profiles obtained in this work
remind those observed in Perry’s experiments, but a deeper study is necessary
in this direction, and this work may represent a starting point.

Acknowledgements:  We are grateful to Eduardo H.F. Colombo for insightful
discussions. We acknowledge partial financial support by the Coordenação de
Aperfeiçoamento de Pessoal de Nível Superior - Brazil (CAPES) - Finance Code
001. CA also acknowledges partial support received from Brazilian agencies
CNPq and Faperj.

## References

  * (1) J. von Hardenberg, E. Meron, M. Shachak, and Y. Zarmi, Diversity of vegetation patterns and desertification,  Phys. Rev. Lett., vol. 87, no. 19, p. 198101, 2001.
  * (2) M. Rietkerk and J. Van de Koppel, Regular pattern formation in real ecosystems,  Trends in Ecology & Evolution, vol. 23, no. 3, pp. 169–175, 2008.
  * (3) M. Rietkerk, M. C. Boerlijst, F. van Langevelde, R. HilleRisLambers, J. van de Koppel, L. Kumar, H. H. T. Prins, and A. M. de Roos, Self-organization of vegetation in arid ecosystems,  The American Naturalist, vol. 160, no. 4, pp. 524–530, 2002.
  * (4) J. D. Murray,  Mathematical Biology: I. An Introduction.  Interdisciplinary Applied Mathematics, Springer, 2002.
  * (5) M. Fuentes, M. Kuperman, and V. Kenkre, Nonlocal interaction effects on pattern formation in population dynamics,  Phys. Rev. Lett., vol. 91, no. 15, p. 158104, 2003.
  * (6) M. Fuentes, M. Kuperman, and V. Kenkre, Analytical considerations in the study of spatial patterns arising from nonlocal interaction effects,  J. Phys. Chem. B, vol. 108, no. 29, pp. 10505–10508, 2004.
  * (7) E. Hernández-García and C. López, Clustering, advection, and patterns in a model of population dynamics with neighborhood-dependent rates,  Phys. Rev. E, vol. 70, no. 1, p. 016216, 2004.
  * (8) D. Escaff, C. Fernandez-Oto, M. Clerc, and M. Tlidi, Localized vegetation patterns, fairy circles, and localized patches in arid landscapes,  Phys. Rev. E, vol. 91, no. 2, p. 022924, 2015.
  * (9) C. Fernandez-Oto, O. Tzuk, and E. Meron, Front instabilities can reverse desertification,  Phys. Rev. Lett., vol. 122, no. 4, p. 048101, 2019\.
  * (10) G. Piva, E. Colombo, and C. Anteneodo, Interplay between scales in the nonlocal FKPP equation,  Chaos, Solitons & Fractals, vol. 153, p. 111609, 2021.
  * (11) J. M. Smith, D. Harper,  et al.,  Animal signals.  Oxford University Press, 2003.
  * (12) L. Corley, T. Staley, L. Bush, and E. Jones, Influence of colostrum on transepithelial movement of Escherichia Coli 055,  Journal of Dairy Science, vol. 60, no. 9, pp. 1416–1421, 1977.
  * (13) C. E. Tarnita, J. A. Bonachela, E. Sheffer, J. A. Guyton, T. C. Coverdale, R. A. Long, and R. M. Pringle, A theoretical foundation for multi-scale regular vegetation patterns,  Nature, vol. 541, no. 7637, pp. 398–401, 2017.
  * (14) S. L. Pimm and C. N. Jenkins, Connecting habitats to prevent species extinctions conservation biologists are creating links between forest fragments where the most animals with small ranges live,  American Scientist, vol. 107, no. 3, pp. 162–169, 2019.
  * (15) G. A. Maciel and F. Lutscher, How individual movement response to habitat edges affects population persistence and spatial spread,  The American Naturalist, vol. 182, no. 1, pp. 42–52, 2013.
  * (16) B. A. Wintle, H. Kujala, A. Whitehead, A. Cameron, S. Veloz, A. Kukkala, A. Moilanen, A. Gordon, P. E. Lentini, N. C. Cadenhead,  et al., Global synthesis of conservation studies reveals the importance of small habitat patches for biodiversity,  Proceedings of the National Academy of Sciences, vol. 116, no. 3, pp. 909–914, 2019.
  * (17) N. Perry, Experimental validation of a critical domain size in reaction–diffusion systems with escherichia coli populations,  Journal of The Royal Society Interface, vol. 2, no. 4, pp. 379–387, 2005.
  * (18) V. Dornelas, E. H. Colombo, C. López, E. Hernández-García, and C. Anteneodo, Landscape-induced spatial oscillations in population dynamics,  Scientific Reports, vol. 11, no. 1, pp. 1–11, 2021.
  * (19) W. R. DiLuzio, L. Turner, M. Mayer, P. Garstecki, D. B. Weibel, H. C. Berg, and G. M. Whitesides, Escherichia Coli swim on the right-hand side,  Nature, vol. 435, no. 7046, pp. 1271–1274, 2005.
  * (20) J. E. Sosa-Hernández, M. Santillán, and J. Santana-Solano, Motility of Escherichia Coli in a quasi-two-dimensional porous medium,  Phys. Rev. E, vol. 95, no. 3, p. 032404, 2017.
  * (21) M. Muskat, The Flow of Homogeneous Fluids through Porous Media, McGraw-Hill, New York, 1937.
  * (22) J. Arlt, V. A. Martinez, A. Dawson, T. Pilizota, and W. C. Poon, Dynamics-dependent density distribution in active suspensions,  Nature communications, vol. 10, no. 1, pp. 1–7, 2019.
  * (23) E. H. Colombo and C. Anteneodo, Nonlinear diffusion effects on biological population spatial patterns,  Phys. Rev. E, vol. 86, no. 3,p. 036215, 2012.
  * (24) M. E. Cates, D. Marenduzzo, I. Pagonabarraga, and J. Tailleur, Arrested phase separation in reproducing bacteria creates a generic route to pattern formation,  Proceedings of the National Academy of Sciences, vol. 107, pp. 11715–11720, June 2010.
  * (25) C. López, Macroscopic description of particle systems with nonlocal density-dependent diffusivity,  Phys. Rev. E, vol. 74, no. 1, p. 012102, 2006.
  * (26) V. Kenkre and N. Kumar, Nonlinearity in bacterial population dynamics: Proposal for experiments for the observation of abrupt transitions in patches,  Proceedings of the National Academy of Sciences, vol. 105, no. 48, pp. 18752–18757, 2008.
  * (27) E. Lenzi, C. Anteneodo, and L. Borland, Escape time in anomalous diffusive media,  Phys. Rev. E, vol. 63, no. 5, p. 051109, 2001.
  * (28) C. Anteneodo, Non-extensive random walks,  Physica A: Statistical Mechanics and its Applications, vol. 358, no. 2-4, pp. 289–298, 2005.
  * (29) C. Anteneodo, Brownian motors in nonlinear diffusive media,  Phys. Rev. E, vol. 76, no. 2, p. 021102, 2007.
  * (30) W. I. Newman, Some exact solutions to a non-linear diffusion problem in population genetics and combustion,  J. Theor. Biol. , vol. 85, no. 2, pp. 325–334, 1980.
  * (31) M. E. Gurtin and R. C. MacCamy, On the diffusion of biological populations,  Math. Biosciences, vol. 33, no. 1-2, pp. 35–49, 1977.
  * (32) P. Kareiva, Local movement in herbivorous insects: applying a passive diffusion model to mark-recapture field experiments,  Oecologia, vol. 57, no. 3, pp. 322–327, 1983.
  * (33) G. Birzu, S. Matin, O. Hallatschek, and K. S. Korolev, Genetic drift in range expansions is very sensitive to density feedback in dispersal and growth,  arXiv preprint arXiv:1903.11627, 2019.
  * (34) J. Tailleur and M. Cates, Statistical mechanics of interacting run-and-tumble bacteria, Phys. Rev. Lett., vol. 100, no. 21, p. 218103, 2008.
  * (35) R. Cantrell and C. Cosner, Diffusion models for population dynamics incorporating individual behavior at boundaries: applications to refuge design,  Theoretical population biology, vol. 55, no. 2, pp. 189–207, 1999.
  * (36) M. Dos Santos, V. Dornelas, E. Colombo, and C. Anteneodo, Critical patch size reduction by heterogeneous diffusion,  Phys. Rev. E, vol. 102, no. 4, p. 042139, 2020.
  * (37) E. H. Colombo and C. Anteneodo, Population dynamics in an intermittent refuge,  Phys. Rev. E, vol. 94, no. 4, p. 042413, 2016.
  * (38) E. Colombo and C. Anteneodo, Nonlinear population dynamics in a bounded habitat,  J. Theor. Biol., vol. 446, pp. 11–18, 2018.
  * (39) D. Ludwig, D. Aronson, and H. Weinberger, Spatial patterning of the spruce budworm,  J. Math. Biol., vol. 8, no. 3, pp. 217–258, 1979.
  * (40) J. Skellam, Random dispersal in theoretical populations,  Bulletin of Mathematical Biology, vol. 53, no. 1-2, pp. 135–165, 1991.
  * (41) V. Kenkre and M. Kuperman, Applicability of the fisher equation to bacterial population dynamics,  Phys. Rev. E, vol. 67, no. 5, p. 051921, 2003\.
  * (42) V. Dornelas, E. H. Colombo, and C. Anteneodo, Single-species fragmentation: The role of density-dependent feedback, Phys. Review E, vol. 99, no. 6, p. 062225, 2019.


# Shortcuts to Adiabaticity with Inherent Robustness and without Auxiliary
Control

Yiyao Liu  Guangdong Provincial Key Laboratory of Quantum Engineering and
Quantum Materials, School of Physics and Telecommunication Engineering, South
China Normal University, Guangzhou 510006, China  Zhen-Yu Wang
<EMAIL_ADDRESS>Guangdong Provincial Key Laboratory of Quantum
Engineering and Quantum Materials, School of Physics and Telecommunication
Engineering, South China Normal University, Guangzhou 510006, China  Frontier
Research Institute for Physics, South China Normal University, Guangzhou
510006, China

###### Abstract

We develop a framework to construct a new class of shortcuts to adiabaticity
for realization of adiabatic processes in short times. In contrast to previous
strategies for accelerating adiabatic protocols, our approach speeds up the
adiabatic process by applying the adiabatic Hamiltonian at discrete points of
the parameter space, thereby avoiding the use of any counterdiabatic control
that can be experimentally unfeasible, skipping undesired points of the
adiabatic path, and retaining the intrinsic robustness of adiabatic control.
We exemplify our general theory by construction of fast and noise-resilient
control protocols for bosonic modes, three-level systems, and interacting
qubits. Our framework offers a new route to design robust and fast control
methods for general quantum systems.

## Introduction. —

Due to its intrinsic robustness against control imperfections, quantum
adiabatic control has a board range of applications from quantum information
processing to quantum sensing, such as stimulated Raman adiabatic passage
(STIRAP) for state transfer in three-level systems Vitanov  _et al._ (2017);
Vasilev  _et al._ (2009); Gaubatz  _et al._ (1990); Bergmann  _et al._ (1998)
and adiabatic quantum computing Albash and Lidar (2018). However, these
quantum adiabatic control methods are extremely slow following the traditional
adiabatic condition. This slowness can be a big problem because decoherence
and noise could have enough time to spoil the desired states.

Shortcuts to adiabaticity (STA) Demirplak and Rice (2003); Berry (2009); Chen
_et al._ (2010) are alternative fast processes to realize the same final state
evolution and have attracted much attention due to practical and fundamental
interests Torrontegui  _et al._ (2013); Guéry-Odelin  _et al._ (2019); del
Campo and Kim (2019); Santos and Sarandy (2017); Hu  _et al._ (2018); Santos
_et al._ (2020); Takahashi (2019); Abah  _et al._ (2020); Kölbl  _et al._
(2019); Chen  _et al._ (2021); Zheng  _et al._ (2022a). In existing STA
methods, the total Hamiltonian $H^{\prime}=H+H_{{\rm CD}}$ is different from
the original time-dependent Hamiltonian $H$ of adiabatic control. The addition
of the counterdiabatic Hamiltonian $H_{{\rm CD}}$ ensures that an eigenstate
of $H(t=0)$ will evolve to the corresponding eigenstate of $H(T)$ at a later
time $T$. However, implementation of $H_{{\rm CD}}$ can be challenging as it
may require forbidden transitions or unavailable experimental resources
Torrontegui  _et al._ (2013). While for some cases it is possible to derive
experimentally feasible STA control using a suitable unitary transformation
Ibáñez  _et al._ (2012), the procedure is complicated even for the case of
two- and three-level systems Ibáñez  _et al._ (2012); Bason  _et al._ (2012);
Martínez-Garaot  _et al._ (2014); Li and Chen (2016); Baksic  _et al._ (2016);
Du  _et al._ (2016); Zhou  _et al._ (2017). Moreover, although a larger
$H_{{\rm CD}}$ lessens the effects of dissipation and decoherence by
shortening the evolution time, $H_{{\rm CD}}$ generally adds additional
sources of control errors and changes the eigenstates of the total Hamiltonian
$H^{\prime}$. The latter point implies that the evolution of an initial ground
state does not follow the ground-state evolution of the total Hamiltonian.
Therefore, the robustness of STA is not guaranteed. Indeed, a recent
comparative study Torosov  _et al._ (2021) on the control of two-level systems
showed that STA performs worse to typical control errors such as amplitude and
detuning errors. See Fig. 1(a). Designing STA protocols robust against control
errors is a non-trivial problem Song  _et al._ (2021); Ruschhaupt  _et al._
(2012); Chen  _et al._ (2010); Zheng  _et al._ (2022a). In Refs. Xu  _et al._
(2019); Zheng  _et al._ (2022b), it was found that for the case of two-level
systems one can, without the use of any $H_{{\rm CD}}$, speed up the quantum
adiabatic evolution and retain the inherent robustness of adiabatic process,
by utilizing the necessary and sufficient quantum adiabatic condition Wang and
Plenio (2016). This shed new light on designing STA protocols in a new way for
general systems.

Figure 1: (a) Schematic trajectories (defined by the eigenstates of the total
Hamiltonian) to the target state (indicated by the blue dashed circle) for an
infinitely slow adiabatic process (solid blue), an ideal STA (solid green), a
STA with control errors (dashed light green), and a STAM control (orange
crosses). STA can be sensitive to control errors $\delta$. STAM and adiabatic
control are robust and have the same instantaneous eigenstates of the total
Hamiltonian. (b) A sketch of the function $F(\lambda)$. The integral for the
regions of the same color is zero. (c) Nodes $n$ and $m$ that are not
connected have $g_{n,m}=0$ in STAM. The colors denote the types of the nodes.
(d) An example of (c) for a bosonic mode.

In this Letter, we develop a general theory to construct a new type of STA
control. This STA by modulation (STAM) method is achieved by dynamically
modulating the functional forms of the parameters in the Hamiltonian for
quantum adiabatic evolution. The resulting protocols are simple, do not
require additional control, retain the robustness of adiabatic processes, and
can avoid the control and interactions that are not feasible in experiments.
We illustrate the idea and show superior performances of STAM by applying it
to bosonic systems, three-level systems, and interacting qubits.

## STAM. —

We write the time-dependent Hamiltonian in its diagonal form ($\hbar=1$)

$H(\lambda)=\sum_{n}E_{n}|n_{\lambda}\rangle\langle n_{\lambda}|,$ (1)

where $E_{n}$ are the eigenenergies, and the instantaneous eigenstates
$\\{|n_{\lambda}\rangle\\}$ depend on a monotonic function
$\lambda=\lambda(t)$ which starts from an initial path point $\lambda(0)=0$.
In the adiabatic limit of an infinitely slow change of $H$, the time-ordered
evolution $U=\mathcal{T}e^{-i\int_{0}^{t}Hdt^{\prime}}$ becomes the adiabatic
evolution operator $U_{{\rm
adia}}=\sum_{n}e^{-i\varphi_{n}}|n_{\lambda}\rangle\langle n_{0}|$. Since the
choice of gauge does not change the Hamiltonian and its evolution operator
Wang and Plenio (2016), for convenience we use the Born-Fock gauge
$g_{n,n}(\lambda)=0$ such that
$\int_{0}^{\lambda}g_{n,n}(\lambda^{\prime})d\lambda^{\prime}=0$, where

$g_{n,m}(\lambda)\equiv\langle
n_{\lambda}|i\frac{d}{d\lambda}|m_{\lambda}\rangle.$ (2)

In this gauge, $\varphi_{n}=\int_{0}^{t}E_{n}(t^{\prime})dt^{\prime}$ equals
the dynamic phase.

To retain the intrinsic robustness of adiabatic process we keep using the
Hamiltonian Eq. (1). We realize the target evolution $U_{{\rm adia}}$ within a
finite time by choosing a proper functional form of varying $H(\lambda)$. In
the Supplemental Material SM  we obtain the fidelity Wang  _et al._ (2008)
between the two evolution operators $U$ and $U_{{\rm adia}}$ as
$\mathcal{F}=|{\rm Tr}U_{D}|/L_{{\rm dim}}\geq 1-|D|$, where the normalization
factor $L_{{\rm dim}}={\rm Tr}I$ equals the trace of the identity operator
$I$; the unitary operator

$U_{D}(\lambda)\equiv U_{{\rm
adia}}^{\dagger}U=\mathcal{P}e^{i\int_{0}^{\lambda}W(\lambda^{\prime})d\lambda^{\prime}},$
(3)

where $W(\lambda)=\sum_{n\neq
m}e^{i[\varphi_{n}(\lambda)-\varphi_{m}(\lambda)]}g_{n,m}(\lambda)|n_{0}\rangle\langle
m_{0}|$ and $\mathcal{P}$ indicates ordering with respect to $\lambda$; the
deviation $D$ from the target evolution has a universal upper bound

$|D(\lambda)|<2\lambda L_{{\rm dim}}\sqrt{\epsilon_{{\rm ave}}\left(g_{{\rm
max}}^{2}L_{{\rm dim}}+g_{{\rm max}}^{\prime}\right)g_{{\rm max}}},$ (4)

which is valid for any Hamiltonian. Here the least upper bounds $g_{{\rm
max}}=\sup_{{\rm\lambda^{\prime}\in(0,\lambda)}}|g_{n,m}(\lambda^{\prime})|$
and $g_{{\rm
max}}^{\prime}=\sup_{{\rm\lambda^{\prime}\in(0,\lambda)}}|\frac{d}{d\lambda^{\prime}}g_{n,m}(\lambda^{\prime})|$
are fixed for a given evolution path. $|D(\lambda)|$ is small when the factor
$\epsilon_{{\rm
ave}}=\sup_{{\rm\lambda^{\prime}\in(0,\lambda)}}|\int_{0}^{\lambda^{\prime}}e^{i[\varphi_{n}(\lambda^{\prime\prime})-\varphi_{m}(\lambda^{\prime\prime})]}d\lambda^{\prime\prime}|$
is small. Note that if $W(\lambda)$ commutes with itself at different values
of $\lambda$ and $\int_{0}^{\lambda}W(\lambda^{\prime})d\lambda^{\prime}=0$,
Eq. (3) gives $U_{D}(\lambda)=I$ and we realize the target evolution $U_{{\rm
adia}}$.

To have $W(\lambda)$ commutes with itself, the first step is to choose

$|n_{\lambda}\rangle=e^{-iG\lambda}|n_{0}\rangle,$ (5)

for the eigenstates $\\{|n_{\lambda}\rangle\\}$ of the Hamiltonian
$H(\lambda)$ using a constant Hermitian operator $G$. In this manner each
$g_{n,m}=\langle
n_{0}|e^{iG\lambda}i\frac{d}{d\lambda}e^{-iG\lambda}|m_{0}\rangle=\langle
n_{0}|G|m_{0}\rangle$ becomes a constant, by using Eq. (2). This also reduces
the bound of the non-adiabatic error $|D|$ in Eq. (4) since $g_{{\rm
max}}^{\prime}=0$. The second step is to modulate the functional forms of the
dynamic phases with respect to $\lambda$ for all pairs of $\\{n,m\\}$ with
$g_{n,m}\neq 0$ such that
$e^{i[\varphi_{n}(\lambda)-\varphi_{m}(\lambda)]}=F(\lambda)$, where
$F(\lambda)\equiv e^{ij\pi}$ when $\lambda\in[\lambda_{j},\lambda_{j+1})$, see
Fig. 1(b). For other pairs of $\\{n,m\\}$ that
$e^{i[\varphi_{n}(\lambda)-\varphi_{m}(\lambda)]}\neq F(\lambda)$, we set
$g_{n,m}=0$. In this way, $W(\lambda)=F(\lambda)G$ commutes with itself, and
from Eq. (3)
$U_{D}(\Theta_{j})=e^{i\int_{0}^{\Theta_{j}}W(\lambda)d\lambda}=I$, where
$\Theta_{j}=2\sum_{k=1}^{j}(-1)^{j+k}\lambda_{k},\ (j=1,2,\ldots,N)$, see Fig.
1(b). Therefore, the controlled evolution $U=U_{{\rm adia}}$ perfectly
realizes the target evolution at $\Theta_{j}$. For simplicity we assume
equally spaced points $\lambda_{j}=\frac{\Theta_{N}}{2N}(2j-1)$ and
$\Theta_{j}=j\Theta_{N}/N$.

While other ways are possible to achieve the desired functional forms of
$\\{\varphi_{n}(\lambda)\\}$ for $F(\lambda)$, here we apply a sequence of
control $H(\lambda)$ with
$\lambda=\lambda_{1},\lambda_{2},\ldots,\lambda_{N}$. In contrast to
traditional adiabatic control where $H(\lambda)$ varies slowly over a
continuous range of $\lambda$, applying $H(\lambda)$ only at
$\\{\lambda_{j}\\}$ can avoid unwanted points $\lambda\notin\\{\lambda_{j}\\}$
that could be challenging or not feasible in experiments. We note that some of
$\\{|n_{\lambda}\rangle\\}$ could be chosen as dark states because only
relative energies are important.

Any given initial and target states can be connected via Eq. (5) with a
suitable choice of $G$. This can be shown by an explicit example that
$G=\sum_{n=2}^{L_{{\rm dim}}}g_{n,1}|n_{0}\rangle\langle 1_{0}|+{\rm h.c.}$.
In this case, we may choose $\varphi_{1}=E_{1}=0$ and, for $n\neq 1$,
$E_{n}=E$, $\varphi_{n}(\lambda)=j\pi$ when
$\lambda\in[\lambda_{j},\lambda_{j+1})$. The sequence, where each
$H(\lambda_{j})$ is applied for a pulse duration of $\pi/E$, transforms an
arbitrary initial state $|1_{0}\rangle$ to a general superposition state
$|1_{\Theta_{N}}\rangle=\cos(g\Theta_{N})|1_{0}\rangle+i\sin(g\Theta_{N})\sum_{n=2}^{L_{{\rm
dim}}}\tilde{g}_{n,1}|n_{0}\rangle$ with $\tilde{g}_{n,1}=g_{n,1}/g$ and
$g=\sqrt{\sum_{n=2}^{L_{{\rm dim}}}|g_{n,1}|^{2}}$.

To find which $g_{n,m}$ can be nonzero, we consider the clusters where the
nodes marked by state labels $n$ and $m$ are connected (disconnected) if
$g_{n,m}\neq 0$ ($g_{n,m}=0$), as illustrated in Fig. 1(c). Since for
$g_{n,m}\neq 0$ the difference of dynamic phases $\varphi_{n}-\varphi_{m}$ has
to be changed by an amount of $(2k_{nm}+1)\pi$ (with $k_{nm}$ being an
integer) at each point $\lambda_{j}$, each cluster $\mathcal{C}$ of a size
larger than one is required to have exactly two types of nodes,
$\mathcal{A}_{\mathcal{C}}$ and $\mathcal{B}_{\mathcal{C}}$. If
$n\in\mathcal{A}_{\mathcal{C}}$ and $m\in\mathcal{B}_{\mathcal{C}}$,
$e^{i[\varphi_{n}(\lambda)-\varphi_{m}(\lambda)]}=F(\lambda)$, while for the
pair $\\{n,m\\}\in\mathcal{A}_{\mathcal{C}}$ or
$\\{n,m\\}\in\mathcal{B}_{\mathcal{C}}$ the nodes $n$ and $m$ cannot be linked
(i.e., $g_{n,m}=0$).

## Fast preparation of coherent states.—

A simple example of choosing the values of $g_{n,m}$ is illustrated for a
bosonic mode in Fig. 1(d), where we choose $g_{n+1,n}=i\alpha\sqrt{n+1}$ in
the Fock state basis $\\{|n\rangle\\}$ ($n\geq 0$). Therefore $G=i\alpha
a^{\dagger}-i\alpha^{*}a$, where $a$ is the annihilation operator.
$e^{-iG\lambda}$ is a displacement operator which transforms the Hamiltonian
$\omega a^{\dagger}a$ of a harmonic oscillator to
$H(\lambda)=e^{-iG\lambda}\left(\omega a^{\dagger}a\right)e^{iG\lambda}=\omega
a^{\dagger}a-\lambda\omega(\alpha
a^{\dagger}+\alpha^{*}a)+\omega|\lambda\alpha|^{2}$. Applying $H(1/2)$ for a
duration of $\pi/\omega$ transfers the original ground state $|0\rangle$ to
the ground state of $H(1)$, i.e., a coherent state
$|\alpha\rangle=e^{-iG}|0\rangle$. This is much faster than adiabatically
varying the control $H(\lambda)$.

Figure 2: Adaptive control for three-level systems. (a) Indirect coupling of
$|0\rangle$ and $|1\rangle$ qubit states via two driving fields with a single-
photon detuning $\Delta$ in a $\Lambda$-type system. (b) The time envelopes of
Stokes $\Omega_{s}(t)$ (blue line) and pump $\Omega_{p}(t)$ (red dash) pulses
of the STAM protocol have simple pulse shapes. The cases of $N=2,3,4$ are
shown. (c) The gate fidelity (yellow triangles) for $N=1$ and
$\Theta_{N}=\frac{\pi}{2}$ in Eq. (6) within the qubit subspace at different
detuning $\Delta$, by including a relaxation rate
$\Gamma_{e}=1.5\Omega/(2\pi)$ for decay of the intermediate state $|e\rangle$
into the qubit subspace and a dephasing rate $\Gamma_{{\rm
dep}}=0.05\Omega/(2\pi)$ on the qubit states. Blue circles (green rectangles)
are the result of setting $\Gamma_{{\rm dep}}=0$ ($\Gamma_{e}=0$). (d)
Robustness of the transfer efficiency as a function of the amplitude and
detuning errors (measured in terms of $\Omega$) for STAM with different $N$.
(e) As (d) but for comparison between the superadiabatic transitionless
driving (SATD) developed in Refs. Zhou  _et al._ (2017); Baksic  _et al._
(2016) and the STAM. Both protocols have the same maximal value of control
field for $\Omega_{p}$ $\Omega_{s}$, and the same length of control time
$T=2\pi/\Omega$. See SM  for more details of the simulations.

## Intrinsic robustness.—

A general, rigorous theory of the robustness for standard adiabatic control
can be seen from Eq. (4). When the control becomes adiabatic, $\epsilon_{{\rm
ave}}$ is small (due to the averaging by
$e^{i[\varphi_{n}(\lambda)-\varphi_{m}(\lambda)]}$) and hence the fidelity
$\mathcal{F}$ is high Wang and Plenio (2016); Xu  _et al._ (2019). It is
obvious that $\epsilon_{{\rm ave}}$ remains small even in the presence of a
positive error on $|\varphi_{n}(\lambda)-\varphi_{m}(\lambda)|$. Thus a more
adiabatic control is more robust against control amplitude errors. For STAM, a
relative error $\delta_{e}$ between ${\varphi_{n}(\lambda)}$ induces
$e^{i(\varphi_{n}-\varphi_{m})}=e^{ij\pi(1+\delta_{e})}$ when
$\lambda\in[\lambda_{j},\lambda_{j+1})$. In SM , we obtain $\epsilon_{{\rm
ave}}=\frac{\Theta_{N}}{2N}[1+\pi|\delta_{e}|+O(\delta_{e}^{2})]$, which
remains small when $N$ is large. This implies that the robustness is stronger
when the distance $\lambda_{j+1}-\lambda_{j}$ between subsequent points is
smaller. Furthermore, the robustness of STAM is enhanced by the constraint
$g_{{\rm max}}^{\prime}=0$ because it gives a smaller deviation
$|D(\lambda)|$.

For STAM there is another mechanism to enhance the robustness against fast
fluctuating errors. In contrast to traditional adiabatic approach that
$\lambda$ changes in time, in STAM $H(\lambda)$ is applied at a fixed
$\lambda=\lambda_{j}$ for a pulse duration $\tau_{p}$. At the same point
$\lambda_{j}$, the random energy drifts $\beta(t)$ at different moments
accumulate as a pulse area error
$\delta=\int_{t_{j}}^{t_{j}+\tau_{p}}\beta(t^{\prime})dt^{\prime}$. A
fluctuating drift $\beta(t)$ with a zero mean implies that the expectation
value of $\delta$ is $\overline{\delta}=0$. This is the same mechanism of
robustness as the geometric gate in Ref. Zu  _et al._ (2014), because fast
fluctuation at the boundary ($\beta$) hardly changes the area ($\delta$). The
variance
$\overline{\delta^{2}}=\int_{t_{j}}^{t_{j}+\tau_{p}}\int_{t_{j}}^{t_{j}+\tau_{p}}dt_{1}dt_{2}\overline{\beta(t_{1})\beta(t_{2})}$
is negligible when the correlation time of the correlation function
$\overline{\beta(t)\beta(0)}$ is small.

## Adaptive control for $\Lambda$-systems.—

To demonstrate STAM, consider a $\Lambda$-type three-level system where there
is no direct coupling between the two qubit states $|0\rangle$ and
$|1\rangle$, see Fig. 2(a). For this system, the adiabatic jumping protocol in
Xu  _et al._ (2019) and the STA protocols Chen  _et al._ (2010) which require
direct coupling between $|0\rangle$ and $|1\rangle$ cannot be used. To realize
robust and fast state transfer between $|0\rangle$ and $|1\rangle$ by a dark
state
$|1_{\lambda}\rangle=e^{-iG\lambda}|1\rangle=\cos(\lambda)|1\rangle-\sin(\lambda)e^{i\phi}|0\rangle$
with $E_{1}=0$, we choose $g_{1,2}=i\cos\xi$, $g_{1,3}=i\sin\xi$, $g_{2,3}=0$,
and the initial eigenstates $|2_{\lambda=0}\rangle=\cos\xi
e^{i\phi}|0\rangle-\sin\xi|e\rangle$, $|3_{\lambda=0}\rangle=\sin\xi
e^{i\phi}|0\rangle+\cos\xi|e\rangle$. See SM  for more details. The
requirement $\langle 0|H|1\rangle=0$ gives
$(E_{3}-E_{2})\cos(2\xi)=(E_{3}+E_{2})$ and the Hamiltonian of STAM,
$H(\lambda)=\Omega(|B_{\lambda}\rangle\langle e|+|e\rangle\langle
B_{\lambda}|)+\Delta|e\rangle\langle e|,$ where the bright state
$|B_{\lambda}\rangle=\sin\lambda|1\rangle+e^{i\phi}\cos\lambda|0\rangle$, the
detuning $\Delta=E_{2}+E_{3}$, and the overall coupling strength
$\Omega=\sqrt{-E_{2}E_{3}}$. To meet the condition for $F(\lambda)$,
$E_{2}=-(2k_{2}+1)\pi/t_{p}$ and $E_{3}=(2k_{3}+1)\pi/t_{p}$ with the integers
$k_{2},k_{3}\geq 0$. The amplitudes of the Stokes
($\Omega_{s}=2\Omega\cos\lambda$) and pump ($\Omega_{p}=2\Omega\sin\lambda$)
pulses are illustrated in Fig. 2(b). Applying $H(\lambda)$ according to the
sequence $\lambda=\lambda_{1},\lambda_{2},\ldots,\lambda_{N}$ realizes a
universal single-qubit quantum gate,

$U=\left(\begin{array}[]{cc}\cos\Theta_{N}&(-1)^{N}e^{-i\phi}\sin\Theta_{N}\\\
-e^{i\phi}\sin\Theta_{N}&(-1)^{N}\cos\Theta_{N}\end{array}\right),$ (6)

in the basis $\\{|1\rangle,|0\rangle\\}$ of the qubit subspace. For even $N$
the gate is a rotation of $2\Theta_{N}$ along a direction in the $x-y$ plane
of the Bloch sphere, while for odd $N$ it is a $\pi$ rotation along an
arbitrary direction.

It is interesting that the value of $\Delta$ is adjustable for optimal
performance. For the resonant case $\Delta=0$, $t_{p}=\pi/\Omega$ reaches the
quantum speed limit Pires  _et al._ (2016); Marvian  _et al._ (2016) and the
protocol for $N=1$ coincides with the geometric gate in Ref. Zu  _et al._
(2014). On the other hand, a larger $\Delta$ can be used to reduce the effect
of dissipation on the intermediate state $|e\rangle$. The traditional method
to suppress this dissipation is the adiabatic elimination of $|e\rangle$ via a
slow, second-order control in the regime $\Omega\ll\Delta$. This slowness is a
large obstacle in the presence of dephasing noise since STIRAP is very
sensitive to the two-photon detuning (equivalent to dephasing noise on the
qubit subspace) Guéry-Odelin  _et al._ (2019). In contrast, STAM can
accelerate the control in the intermediate regime $\Delta\sim\Omega$ to
overcome this obstacle, offering a higher gate fidelity in the presence of
both dephasing and dissipation [see Fig. 2(c)].

Comparing with existing robust STA for $\Lambda$ systems Zhou  _et al._
(2017); Baksic  _et al._ (2016), the robustness of STAM is much better and can
be further enhanced by using a larger $N$ [see Fig. 2(d),(e)].

Figure 3: (a) The energy levels of Eq. (7) are degenerate at the point $s=1$,
where the target ground state encodes the solution. (b) The fidelity (blue
solid line) of standard adiabatic evolution to the target ground state at
$s=1$ as a function of the total evolution time $T$. The fidelity of the STAM
protocol with $N=1$ is shown by the orange dashed line and the minimal
required time $T$ is indicated by an orange circle. (c) The fidelity as a
function of the strength of a local perturbation
$\epsilon_{x}\sigma_{x}^{(1)}$ for standard adiabatic control (blue solid and
green dash-dotted lines) and STAM (orange dashed line). (d) As (b) but with
the addition of a weak local perturbation $0.05E\sigma_{x}^{(1)}$.

## Application to coupled qubits.—

As another example, consider the Hamiltonian for adiabatic quantum
computation,

$H=(1-s)H_{0}+sH_{1},$ (7)

where $H_{0}=-\frac{E}{2}(\sigma_{z}^{(1)}+\sigma_{z}^{(2)})$ and
$H_{1}=-E\sigma_{x}^{(1)}\sigma_{x}^{(2)}$ ($E>0$) Xu  _et al._ (2017).
Starting from the ground state of $H_{0}$,
$|11\rangle\equiv|1\rangle\otimes|1\rangle$, slowly increasing the value of
$s=s(t)$ from $s=0$ to $s=1$ is supposed to adiabatically prepare the target
solution state $|\psi_{+}\rangle=\frac{1}{\sqrt{2}}(|00\rangle-|11\rangle)$.
However, $H$ has a gap closing at the final point $s=1$, with another ground
state $|\psi_{-}\rangle=\frac{1}{\sqrt{2}}(|01\rangle-|10\rangle)$, see Fig.
3(a). Therefore, a local fluctuation (e.g., a weak perturbation
$\epsilon_{x}\sigma_{x}^{(1)}$ with $\epsilon_{x}>0$) breaks the parity
symmetry of $H$ and leads to a wrong result (e.g.,
$|\psi_{x}\rangle=\frac{1}{\sqrt{2}}(|\psi_{+}\rangle+|\psi_{-}\rangle$) in
the adiabatic limit, as shown in Fig. 3.

STAM avoids the use of the Hamiltonian at the degenerate point $s=1$ and is
much faster. For this case we use $g_{1,2}=g_{2,1}=-1$ (other $g_{n,m}=0$) and
$|2_{0}\rangle=i|00\rangle$ such that $|1_{\lambda}\rangle$ connects the
initial ground state $|1_{0}\rangle=|11\rangle$ and the target ground state
$|\psi_{+}\rangle$. The other initial eigenstates are the general
superposition states $|3_{0}\rangle=\cos\beta
e^{i\xi}|01\rangle-\sin\beta|10\rangle$ and $|4_{0}\rangle=\sin\beta
e^{i\xi}|01\rangle+\cos\beta|10\rangle$. Setting $\xi=0$, $\beta=\pi/4$,
$E_{1}=-E_{2}=E$, $E_{3}=-E_{4}=E\sin(2\lambda)$, we obtain the Hamiltonian
$H(\lambda)=\cos(2\lambda)H_{0}+\sin(2\lambda)H_{1}$ as Eq. (7) but with a
different interpolation. As shown in Fig. 3, compared with standard adiabatic
control with a linear interpolation of $s=t/T$, STAM with $N=1$ achieves a
higher state fidelity using a shorter time $\pi/(2E)$ and is much more robust
to local fluctuations because the degenerate point $\lambda=\pi/4$ (i.e.,
$s=1$) is avoided.

We can also use a different Hamiltonian, say,
$H(\lambda)=\cot(2\lambda)H_{0}+H_{1}$, which has similarities to the model of
Landau-Zener-Stückelberg transitions Shevchenko  _et al._ (2010). In contrast
to the adiabatic control that continuously varies $\lambda$ from $\lambda=0$
to $\lambda=\pi/4$, our control is bounded because the divergent point
$\lambda=0$ is avoided. Because the interaction $H_{1}$ is a constant and
$\cot(2\lambda)H_{0}$ can be varied by the detuning of control field, this
control is faster and is easy to realize in systems such as a nitrogen-vacancy
(NV) center in diamond and nuclear spins Taminiau  _et al._ (2012); Wang  _et
al._ (2016, 2017); Bradley  _et al._ (2019).

## Conclusion.—

We developed a general theory for general quantum systems to construct a new
kind of STA protocols that speed up the adiabatic evolution only by the use of
interactions in the original adiabatic Hamiltonian. Without the requirement of
counterdiabatic control or its equivalence using unitary transformations, the
resulting STAM protocols are more experimentally feasible and can avoid
obstacles due to degenerate points, unbounded values of control fields, and
unavailable control resources. In addition to the ability to retain the
intrinsic robustness of quantum adiabatic control, STAM protocols have other
mechanisms to suppress the effects of control errors and noise from the
environment.

Our results provide a new route to design robust and high-speed control
methods, as demonstrated by the examples of preparation of coherent states,
reliable evolution to the solution state in an adiabatic quantum computing
model, as well as a three-level system control protocol which has adjustable
values of detuning for optimal performance in the presence of both dephasing
noise and excited state decay. In future work, the pulse-sequence like STAM
could be modified to have the ability of coherence protection as dynamical-
decoupling (DD) pulse sequences Viola  _et al._ (1999); Yang  _et al._ (2011);
Xu  _et al._ (2019) and may be incorporated in DD-based sensing with adiabatic
control Haase  _et al._ (2018); Whaites  _et al._ (2022).

###### Acknowledgements.

This work was supported by National Natural Science Foundation of China (Grant
No. 12074131) and the Natural Science Foundation of Guangdong Province (Grant
No. 2021A1515012030). We thank Dr. Min Ming for useful suggestions and
discussions.

## References

  * 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).
  * Vasilev  _et al._ (2009) G. S. Vasilev, A. Kuhn, and N. V. Vitanov, “Optimum pulse shapes for stimulated raman adiabatic passage,” Phys. Rev. A 80, 013417 (2009).
  * Gaubatz  _et al._ (1990) U. Gaubatz, P. Rudecki, S. Schiemann, and K. Bergmann, “Population transfer between molecular vibrational levels by stimulated raman scattering with partially overlapping laser fields. a new concept and experimental results,” J. Chem. Phys. 92, 5363 (1990).
  * Bergmann  _et al._ (1998) K. Bergmann, H. Theuer, and B. W. Shore, “Coherent population transfer among quantum states of atoms and molecules,” Rev. Mod. Phys. 70, 1003 (1998).
  * Albash and Lidar (2018) T. Albash and D. A. Lidar, “Adiabatic quantum computation,” Rev. Mod. Phys. 90, 015002 (2018).
  * Demirplak and Rice (2003) M. Demirplak and S. A. Rice, “Adiabatic population transfer with control fields,” J. Phys. Chem. A 107, 9937 (2003).
  * Berry (2009) M. V. Berry, “Transitionless quantum driving,” J. Phys. A 42, 365303 (2009).
  * Chen  _et al._ (2010) X. Chen, I. Lizuain, A. Ruschhaupt, D. Guéry-Odelin, and J. G. Muga, “Shortcut to adiabatic passage in two- and three-level atoms,” Phys. Rev. Lett. 105, 123003 (2010).
  * 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,” Adv. At. Mol. Opt. Phys. 62, 117 (2013).
  * 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).
  * del Campo and Kim (2019) A. del Campo and K. Kim, “Focus on shortcuts to adiabaticity,” New J. Phys. 21, 050201 (2019).
  * Santos and Sarandy (2017) A. C. Santos and M. S. Sarandy, “Generalized shortcuts to adiabaticity and enhanced robustness against decoherence,” J. Phys. A 51, 025301 (2017).
  * Hu  _et al._ (2018) C.-K. Hu, J.-M. Cui, A. C. Santos, Y.-F. Huang, M. S. Sarandy, C.-F. Li, and G.-C. Guo, “Experimental implementation of generalized transitionless quantum driving,” Opt. Lett. 43, 3136 (2018).
  * Santos  _et al._ (2020) A. C. Santos, A. Nicotina, A. M. Souza, R. S. Sarthour, I. S. Oliveira, and M. S. Sarandy, “Optimizing nmr quantum information processing via generalized transitionless quantum driving,” EPL (Europhys. Lett.) 129, 30008 (2020).
  * Takahashi (2019) K. Takahashi, “Hamiltonian engineering for adiabatic quantum computation: Lessons from shortcuts to adiabaticity,” J. Phys. Soc. Jpn. 88, 061002 (2019).
  * Abah  _et al._ (2020) O. Abah, R. Puebla, and M. Paternostro, “Quantum state engineering by shortcuts to adiabaticity in interacting spin-boson systems,” Phys. Rev. Lett. 124, 180401 (2020).
  * Kölbl  _et al._ (2019) J. Kölbl, A. Barfuss, M. S. Kasperczyk, L. Thiel, A. A. Clerk, H. Ribeiro, and P. Maletinsky, “Initialization of single spin dressed states using shortcuts to adiabaticity,” Phys. Rev. Lett. 122, 090502 (2019).
  * Chen  _et al._ (2021) Y.-H. Chen, W. Qin, X. Wang, A. Miranowicz, and F. Nori, “Shortcuts to adiabaticity for the quantum rabi model: Efficient generation of giant entangled cat states via parametric amplification,” Phys. Rev. Lett. 126, 023602 (2021).
  * Zheng  _et al._ (2022a) W. Zheng, Y. Zhang, Y. Dong, J. Xu, Z. Wang, X. Wang, Y. Li, D. Lan, J. Zhao, S. Li, X. Tan, and Yu. Y., “Optimal control of stimulated raman adiabatic passage in a superconducting qudit,” npj Quantum Inf. 8, 9 (2022a).
  * 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).
  * 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 (2012).
  * 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).
  * Li and Chen (2016) Y.-C. Li and X. Chen, “Shortcut to adiabatic population transfer in quantum three-level systems: Effective two-level problems and feasible counterdiabatic driving,” Phys. Rev. A 94, 063411 (2016).
  * Baksic  _et al._ (2016) A. Baksic, H. Ribeiro, and A. A. Clerk, “Speeding up adiabatic quantum state transfer by using dressed states,” Phys. Rev. Lett. 116, 230503 (2016).
  * Du  _et al._ (2016) Y.-X. Du, Z.-T. Liang, Y.-C. Li, X.-X. Yue, Q.-X. Lv, W. Huang, X. Chen, H. Yan, and S.-L. Zhu, “Experimental realization of stimulated raman shortcut-to-adiabatic passage with cold atoms,” Nat. Commun. 7, 12479 (2016).
  * Zhou  _et al._ (2017) B. B. Zhou, A. Baksic, H. Ribeiro, C. G. Yale, F. J. Heremans, P. C. Jerger, A. Auer, G. Burkard, A. A. Clerk, and D. D. Awschalom, “Accelerated quantum control using superadiabatic dynamics in a solid-state lambda system,” Nat. Phys. 13, 330 (2017).
  * Torosov  _et al._ (2021) B. T. Torosov, B. W. Shore, and N. V. Vitanov, “Coherent control techniques for two-state quantum systems: A comparative study,” Phys. Rev. A 103, 033110 (2021).
  * Song  _et al._ (2021) X.-K. Song, F. Meng, B.-J. Liu, D. Wang, L. Ye, and M.-H. Yung, “Robust stimulated raman shortcut-to-adiabatic passage with invariant-based optimal control,” Opt. Express 29, 7998 (2021).
  * Ruschhaupt  _et al._ (2012) A. Ruschhaupt, X. Chen, D. Alonso, and J. G. Muga, “Optimally robust shortcuts to population inversion in two-level quantum systems,” New J. Phys. 14, 093040 (2012).
  * Xu  _et al._ (2019) K. Xu, T. Xie, F. Shi, Z.-Y. Wang, X. Xu, P. Wang, Y. Wang, M. B. Plenio, and J. Du, “Breaking the quantum adiabatic speed limit by jumping along geodesics,” Sci. Adv. 5, eaax3800 (2019).
  * Zheng  _et al._ (2022b) W. Zheng, J. Xu, Z. Wang, Y. Dong, D. Lan, X. Tan, and Y. Yu, “Accelerated quantum adiabatic transfer in superconducting qubits,” Phys. Rev. Applied 18, 044014 (2022b).
  * Wang and Plenio (2016) Z.-Y. Wang and M. B. Plenio, “Necessary and sufficient condition for quantum adiabatic evolution by unitary control fields,” Phys. Rev. A 93, 052107 (2016).
  * (33) See Supplemental Material [url] for more details.
  * Wang  _et al._ (2008) X. Wang, C.-S. Yu, and X. X. Yi, “An alternative quantum fidelity for mixed states of qudits,” Phys. Lett. A 373, 58 (2008).
  * Zu  _et al._ (2014) C. Zu, W.-B. Wang, L. He, W.-G. Zhang, C.-Y. Dai, F. Wang, and L.-M. Duan, “Experimental realization of universal geometric quantum gates with solid-state spins,” Nature 514, 72 (2014).
  * Pires  _et al._ (2016) D. P. Pires, M. Cianciaruso, L. C. Céleri, G. Adesso, and D. O. Soares-Pinto, “Generalized geometric quantum speed limits,” Phys. Rev. X 6, 021031 (2016).
  * Marvian  _et al._ (2016) I. Marvian, R. W. Spekkens, and P. Zanardi, “Quantum speed limits, coherence, and asymmetry,” Phys. Rev. A 93, 052331 (2016).
  * Xu  _et al._ (2017) K. Xu, T. Xie, Z. Li, X. Xu, M. Wang, X. Ye, F. Kong, J. Geng, C. Duan, F. Shi, and J. Du, “Experimental adiabatic quantum factorization under ambient conditions based on a solid-state single spin system,” Phys. Rev. Lett. 118, 130504 (2017).
  * Shevchenko  _et al._ (2010) S. N. Shevchenko, S. Ashhab, and F. Nori, “Landau–Zener–Stückelberg interferometry,” Phys. Rep. 492, 1 (2010).
  * Taminiau  _et al._ (2012) T. H. Taminiau, J. J. T. Wagenaar, T. van der Sar, F. Jelezko, V. V. Dobrovitski, and R. Hanson, “Detection and control of individual nuclear spins using a weakly coupled electron spin,” Phys. Rev. Lett. 109, 137602 (2012).
  * Wang  _et al._ (2016) Z.-Y. Wang, J. F. Haase, J. Casanova, and M. B. Plenio, “Positioning nuclear spins in interacting clusters for quantum technologies and bioimaging,” Phys. Rev. B 93, 174104 (2016).
  * Wang  _et al._ (2017) Z.-Y. Wang, J. Casanova, and M. B. Plenio, “Delayed entanglement echo for individual control of a large number of nuclear spins,” Nat. Commun. 8, 14660 (2017).
  * Bradley  _et al._ (2019) C. E. Bradley, J. Randall, M. H. Abobeih, R. C. Berrevoets, M. J. Degen, M. A. Bakker, M. Markham, D. J. Twitchen, and T. H. Taminiau, “A ten-qubit solid-state spin register with quantum memory up to one minute,” Phys. Rev. X 9, 031045 (2019).
  * Viola  _et al._ (1999) L. Viola, E. Knill, and S. Lloyd, “Dynamical decoupling of open quantum systems,” Phys. Rev. Lett. 82, 2417 (1999).
  * Yang  _et al._ (2011) W. Yang, Z.-Y. Wang, and R.-B. Liu, “Preserving qubit coherence by dynamical decoupling,” Front. Phys. 6, 2 (2011).
  * Haase  _et al._ (2018) J. F. Haase, Z.-Y. Wang, J. Casanova, and M. B. Plenio, “Soft quantum control for highly selective interactions among joint quantum systems,” Phys. Rev. Lett. 121, 050402 (2018).
  * Whaites  _et al._ (2022) O. T. Whaites, J. Randall, T. H. Taminiau, and T. S. Monteiro, “Adiabatic dynamical-decoupling-based control of nuclear spin registers,” Phys. Rev. Research 4, 013214 (2022).