# MonarchNet: Differentiating Monarch Butterflies from Butterflies Species
with Similar Phenotypes

Thomas Y. Chen
The Academy for Mathematics, Science, and Engineering
<EMAIL_ADDRESS>

###### Abstract

In recent years, the monarch butterfly’s iconic migration patterns have come
under threat from a number of factors, from climate change to pesticide use.
To track trends in their populations, scientists as well as citizen scientists
must identify individuals accurately. This is uniquely key for the study of
monarch butterflies because there exist other species of butterfly, such as
viceroy butterflies, that are "look-alikes" (coined by the Convention on
International Trade in Endangered Species of Wild Fauna and Flora), having
similar phenotypes. To tackle this problem and to aid in more efficient
identification, we present MonarchNet, the first comprehensive dataset
consisting of butterfly imagery for monarchs and five look-alike species. We
train a baseline deep-learning classification model to serve as a tool for
differentiating monarch butterflies and its various look-alikes. We seek to
contribute to the study of biodiversity and butterfly ecology by providing a
novel method for computational classification of these particular butterfly
species. The ultimate aim is to help scientists track monarch butterfly
population and migration trends in the most precise and efficient manner
possible.

##  1 Introduction

The famous monarch butterflies (Danaus plexippus) and their migration patterns
are an iconic phenomenon of the natural world. The species’ unique behavior
has warranted a variety of in-depth studies over the years
overwintering1989reproductive ; wassenaar1998natal ; mouritsen2002virtual ;
reppert2010navigational ; heinze2011sun ; kass2020biotic . Monarch butterflies
are important pollinators and crucial members of the ecosystems to which they
contribute landis2014create . However, in recent times, there have been
concerns about dramatic declines in population and disruptions to their
migrations, caused by a variety of reasons, including climate change,
deforestation, droughts, loss of habitat, and loss of milkweed
brower2012decline ; agrawal2019advances ; thogmartin2017monarch . Notably,
citizen science has been crucial to discovering and monitoring these troubling
declines schultz2017citizen . In order to gather data and gain important
insights into monarch butterfly trends, scientists and citizen scientists must
first identify the insects accurately and efficiently.

Monarch butterflies have many counterparts that have similar external
appearances, from viceroy butterflies (Limenitis archippus) to soldier
butterflies (Danaus eresimus). The five species displayed in Figure 1 are of
particular interest. In some cases, particularly for the viceroy butterfly,
these similarities are due to mimicry (specifically Müllerian mimicry)
alexander2019mechanism . In this case, the main visual difference is the black
line across the viceroy’s hind wings. Hence, recognizing monarch butterflies
correctly and not unknowingly observing an entirely different species by
mistake can be a challenge when studying the species. To approach this problem
with machine learning, we must train deep neural networks on butterfly image
data to enable instantaneous and accurate classification.

Currently, there are not any image datasets that incorporate both monarch
butterflies and their close neighbors in physical appearance. Thus, in this
paper, we present a novel butterfly imagery dataset, MonarchNet, for the
purpose of classifying butterflies that have a similar physical phenotype. We
also train a baseline convolutional neural network (CNN) on this dataset for
assessment.

##  2 Dataset Overview

MonarchNet contains labeled images of monarch butterflies (Danaus plexippus)
and of species that have a similar physical appearance: viceroy butterflies
(Limenitis archippus), red admirals (Vanessa atalanta), painted ladies
(Vanessa cardui), queen butterflies (Danaus gilippus), and soldier butterflies
(Danaus eresimus). These categories contain 200192, 98765, 53654, 26578,
12534, and 14213 images, respectively, with each image containing exactly one
butterfly.

(a) Monarch Butterfly (Danaus plexippus)

(b) Viceroy Butterfly (Limenitis archippus)

(c) Red Admiral (Vanessa atalanta)

(d) Painted Lady (Vanessa cardui)

(e) Queen Butterfly (Danaus gilippus)

(f) Soldier Butterfly (Danaus eresimus)

Figure 1: Sample images of a monarch butterfly and members of 5 phenotypically
similar species

##  3 Dataset Curation

Image data for the MonarchNet dataset are sourced and downloaded from Google
Images. This was chosen because it was determined to be the most optimal way
to feasibly gather a large and sufficient volume of data for each of the six
butterfly species categories. Only images with the Creative Commons Licenses
were used and image urls are saved. Queries of both species’ common and
scientific name (i.e. viceroy butterfly as well as Limenitis archippus) were
used to obtain images. For all queries (scientific or colloquial), we label
the downloaded images with the relevant scientific name. Duplicate or near-
duplicate images are discarded using image hashing. We obtain 405,936 images
in total, with a distribution displayed in Figure 2. Additionally, bounding
boxes for the butterflies are annotated and provided as labels.

Figure 2: Distribution of Butterfly Species Labels

However, images sourced from Google Images are noisy and must be filtered. We
manually clean the data and verify that images fit their respective labels.
This is done through re-examining the image source or verifying specific
characteristics relevant to the corresponding species.

##  4 Baseline Model

We train a baseline CNN to evaluate the dataset and its usefulness for the
monarch butterfly versus "look-alike" species classification problem. The
input consists of a labeled butterfly image, with the species and bounding box
coordinates. The output is a digit from 0 to 5 representing the predicted
species of the butterfly in the image. 0, 1, 2, 3, 4, and 5 correspond to the
classes a, b, c, d, e, and f shown in Figure 1, respectively. The model
architecture is ResNet50 he2016deep , pretrained on ImageNet data
deng2009imagenet . The criterion for optimization is the cross-entropy loss
function. We train on a randomly selected 80% of the dataset with a batch size
of 32. The Adam optimizer with a learning rate of 0.01 is utilized. The
network is trained for 100 epochs on NVIDIA Tesla K80 GPUs.

###  4.1 Results and Discussion

Testing on the remaining 20% of the data, our baseline model achieves a
weighted F1 score of 0.824. F1 scores for individual categories are displayed
in Table 1. This measure of evaluation is used due to the imbalanced
distribution of the dataset. For example, if simple accuracy was our metric,
if the model always predicted the monarch butterfly (Danaus plexippus)
category for each image input, we would see an artificially inflated
percentage because images with that particular label make up a
disproportionately large part of the dataset (Fig. 2).

Table 1: Baseline Model F1 Scores Butterfly Species | F1 Score
---|---
Danaus plexippus (0) | 0.795
Limenitis archippus (1) | 0.761
Vanessa atalanta (2) | 0.962
Vanessa cardui (3) | 0.975
Danaus gilippus (4) | 0.823
Danaus eresimus (5) | 0.854

From Table 1, we see that the model performs best, as per F1 score, on Vanessa
atalanta and Vanessa cardui. This is a reasonable result because these two
species are more physically differentiable, for both humans and AI, from the
other species included in MonarchNet. A brief glance at Figure 1 reveals
relatively significant phenotypical differences. The CNN performs the worst on
images with the ground truth label Danaus plexippus and Limenitis archippus.
This is also justifiable because the monarch and viceroy butterflies are
arguably the most similar in characteristics and may cause confusion. The
remaining two species lie in the middle in terms of the model’s performance.

##  5 Conclusion

In this work, we presented the MonarchNet dataset and a corresponding baseline
CNN to provide a machine-learning-based approach to visually differentiating
between monarch butterflies and five other species that appear similarly. This
research contributes to monarch butterfly research as well as the study of the
biodiversity and ecology of relevant ecosystems by introducing a novel method
for tackling this problem. This is crucial to the monitoring of the decline of
monarch butterflies and the assessment and mitigation of threats to
populations. Our AI method is especially helpful for citizen scientists, who
are vital to these efforts. Effective conservation strategies and habitat
restoration efforts can only be developed if accurate data is gathered. Future
work may include curating datasets that incorporate a larger range of
butterfly species, in order to aid in conservation efforts around the globe,
not solely limited to monarchs (for example, in the Amazon rainforest). The
computer vision community is encouraged to strive to improve upon the baseline
model presented here when the dataset and code are released.

## References

  * (1) Anurag A Agrawal.  Advances in understanding the long-term population decline of monarch butterflies.  Proceedings of the National Academy of Sciences, 116(17):8093–8095, 2019.
  * (2) VN Alexander.  The mechanism for mimicry: Instant biosemiotic selection or gradual darwinian fine-tuning selection?  Biosemiotics, 12(1):39–55, 2019.
  * (3) Lincoln P Brower, Orley R Taylor, Ernest H Williams, Daniel A Slayback, Raul R Zubieta, and M Isabel Ramirez.  Decline of monarch butterflies overwintering in mexico: is the migratory phenomenon at risk?  Insect Conservation and Diversity, 5(2):95–100, 2012.
  * (4) Jia Deng, Wei Dong, Richard Socher, Li-Jia Li, Kai Li, and Li Fei-Fei.  Imagenet: A large-scale hierarchical image database.  In 2009 IEEE Conference on Computer Vision and Patter Recognition, pages 248–255. IEEE, 2009.
  * (5) Kaiming He, Xiangyu Zhang, Shaoqing Ren, and Jian Sun.  Deep residual learning for image recognition.  In Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, pages 770–778. IEEE, 2016.
  * (6) Stanley Heinze and Steven M Reppert.  Sun compass integration of skylight cues in migratory monarch butterflies.  Neuron, 69(2):345–358, 2011.
  * (7) William S. Herman, Lincoln P. Brower, and William H. Calvert.  Reproductive tract development in monarch butterflies overwintering in california and mexico.  Journal of the Lepidopterists’ Society, 43(1):50–58, 1989.
  * (8) Jamie M Kass, Robert P Anderson, Alejandro Espinosa-Lucas, Verónica Juárez-Jaimes, Esteban Martínez-Salas, Francisco Botello, Gloria Tavera, José Juan Flores-Martínez, and Víctor Sánchez-Cordero.  Biotic predictors with phenological information improve range estimates for migrating monarch butterflies in mexico.  Ecography, 43(3):341–352, 2020.
  * (9) Thomas D Landis, R Kasten Dumroese, and Matthew E Horning.  Create a pollinator garden at your nursery: an emphasis on monarch butterflies.  U.S. Forest Service - Forest Nursery Notes, 34(1-2):2–13, 2014\.
  * (10) Henrik Mouritsen and Barrie J Frost.  Virtual migration in tethered flying monarch butterflies reveals their orientation mechanisms.  Proceedings of the National Academy of Sciences, 99(15):10162–10166, 2002.
  * (11) Steven M Reppert, Robert J Gegear, and Christine Merlin.  Navigational mechanisms of migrating monarch butterflies.  Trends in neurosciences, 33(9):399–406, 2010.
  * (12) Cheryl B Schultz, Leone M Brown, Emma Pelton, and Elizabeth E Crone.  Citizen science monitoring demonstrates dramatic declines of monarch butterflies in western north america.  Biological Conservation, 214:343–346, 2017.
  * (13) Wayne E Thogmartin, Ruscena Wiederholt, Karen Oberhauser, Ryan G Drum, Jay E Diffendorfer, Sonia Altizer, Orley R Taylor, John Pleasants, Darius Semmens, Brice Semmens, et al.  Monarch butterfly population decline in north america: identifying the threatening processes.  Royal Society open science, 4(9):170760, 2017.
  * (14) Leonard I Wassenaar and Keith A Hobson.  Natal origins of migratory monarch butterflies at wintering colonies in mexico: new isotopic evidence.  Proceedings of the National Academy of Sciences, 95(26):15436–15439, 1998.


# Exponential scaling of clock stability with atom number

T. Rosenband<EMAIL_ADDRESS>D. R. Leibrandt  National Institute of
Standards and Technology, 325 Broadway, Boulder, CO 80305

###### Abstract

In trapped-atom clocks, the primary source of decoherence is often the phase
noise of the oscillator. For this case, we derive theoretical performance
gains by combining several atomic ensembles. For example, $M$ ensembles of $N$
atoms can be combined with a variety of probe periods, to reduce the frequency
variance to $M2^{-M}$ times that of standard Ramsey clocks. A similar
exponential improvement is possible if the atomic phases of some of the
ensembles evolve at reduced frequencies. These ensembles may be constructed
from atoms or molecules with lower-frequency transitions, or generated by
dynamical decoupling. The ensembles with reduced frequency or probe period are
responsible only for counting the integer number of $2\pi$ phase wraps, and do
not affect the clock’s systematic errors. Quantum phase measurement with
Gaussian initial states allows for smaller ensemble sizes than Ramsey
spectroscopy.

Atomic clocks have many technical and scientific applications, and also serve
as model systems for quantum metrology. An important result of quantum
mechanics is that the variance of phase measurements can be reduced if
entangled particles are used rather than an equal number of unentangled ones
Caves1981 ; Giovannetti2004 . This result has been applied to the theoretical
analysis of atomic clocks that operate by repeatedly measuring the phase
difference between atomic qubits and a classical oscillator
WMI1993ProjectionNoise ; Huelga1997 ; DJW1998bible ; Buzek1999Optimal ;
riis2004optimum ; Lukin2004Stability ; Rosenband2011NumTest ;
Mullan2012Semidefinite . However, these analyses have not taken full advantage
of the hybrid nature of passive atomic clocks, where a classical oscillator
with limited coherence reproduces the resonance frequency of more coherent
quantum systems. The present model derives an exponential gain from additional
couplings between atoms and the laboratory system, to digitize their relative
phase evolution. It is related to classical analog-to-digital conversion
techniques kester2005data  and digital computation of the $arctan$ function
Volder1959CORDIC , which have been extended to the use of qubits for the
measurement of classical fields Vaidman2004  and mesoscopic spins Giedke2006 .
Here it is shown that if the information from $M$ ensembles of $N$ atoms is
combined classically, the frequency variance may be reduced to $MD^{1-M}/2$
times the variance of standard Ramsey clocks with $MN$ atoms, where the value
of $D$ is constrained by the ensemble size. The method works by extending the
range of invertible phases from the usual $-\pi\ldots\pi$ interval to
$-D^{M-1}\pi\ldots D^{M-1}\pi$. One ensemble is used for normal phase
measurements, while the other $M-1$ ensembles are measured with free-evolution
periods that are reduced by factors $D^{-1},\ldots,D^{1-M}$, or evolve at
frequencies that are reduced by the same factors. This allows unambiguous
counting of the number of $2\pi$ phase wraps that occur in the first ensemble.

The passive atomic clocks considered here operate by periodically measuring
the differential phase evolution between a classical oscillator and an
ensemble of atoms. These measurements yield corrections to the oscillator,
which is adjusted in a feedback loop to produce an output signal whose
frequency approximates that of the ideal unperturbed atomic resonance. The
stability of atomic clocks, which describes the uncertainty of this
approximation as a function of the averaging duration, is limited by quantum
projection noise WMI1993ProjectionNoise  for short durations. If a clock
operates by repeatedly applying the Ramsey technique Ramsey1956  to measure
the phase difference of $N$ unentangled atoms with respect to the oscillator,
the projection-noise-limited fractional-frequency stability is
$\sigma_{yp}(\tau)=(\omega\sqrt{NT\tau})^{-1}$, where $T$ is the free-
evolution period, $\tau$ is the total duration, and $\omega$ is the angular
frequency of the atoms WMI1993ProjectionNoise ; Huelga1997 . This stability
limit is an important figure of merit, and it connects atomic clocks to
quantum metrology, because the use of entangled states allows
$\sigma_{yp}(\tau)$ to scale more favorably as $N^{-1}$. A clock that operates
for duration $\tau$ has a fractional time-variance of
$var(\tau)/\tau^{2}=\sigma_{yp}(\tau)^{2}+\sigma_{ya}(\tau)^{2}$, where
$\sigma_{ya}(\tau)$ is the accuracy, which does not decrease below the bounds
set by the systematic errors. For short measurement durations, quantum
measurement noise restricts the degree to which atomic clocks can reproduce
the ideal frequency $\omega$. For long durations, the dominant uncertainty
comes from systematic errors due to effects such as unknown environmental
factors.

Proposals to use the available atoms more efficiently, and thereby reduce
$\sigma_{yp}$, include the use of spin-squeezed states Kitagawa1993Squeezed ;
Wineland1994Squeezed ; Lukin2004Stability ; Vuletic2010Squeeze ; more general
quantum measurements are also possible JJB1996Entanglement ; Buzek1999Optimal
; Dobrzanski2011 ; Rosenband2011NumTest ; Mullan2012Semidefinite . Such
approaches are restricted by the fact that the atomic phase evolves as
$\phi=\int^{T}_{0}(\omega-\Omega(t))dt$ with respect to the classical
oscillator, whose frequency is $\Omega(t)$. When the atoms evolve freely,
there is no way to unambiguously resolve phases beyond $\pm\pi$. Because of
this, $T$ is usually constrained to a duration that is short enough so that
the atom-oscillator phase difference has not drifted outside the $-\pi$ to
$\pi$ range of invertibility. Here it is shown that the range of phase
invertibility and therefore $T$ can grow exponentially when multiple atomic
ensembles are available.

Consider a passive atomic clock as described above. Atomic decoherence is
assumed to be negligible compared to oscillator decoherence. This idealization
approximates many current atomic clocks, where the oscillator is a laser whose
output has a $1/f$ power spectrum of frequency fluctuations
Numata2004ThermalNoise . Then the phase difference $\phi$ can be thought of as
a Gaussian random variable with standard-deviation $\sigma_{\phi}=\alpha T$,
where $\alpha$ characterizes the oscillator noise level. Only phases inside
some interval $-\theta\ldots\theta$ can be unambiguously inverted, and it is
reasonable to require that $\sigma_{\phi}\leq\theta/6$, making phase-wrapping
errors six-sigma events. Other choices of the ratio $\sigma_{\phi}/\theta$ are
also possible, to reduce this probability to a level that is considered
insignificant. In this Letter, such a clock is called a standard Ramsey clock.
Here the limit is $\theta=\pi/2$ and $T=\pi/{12\alpha}$, because the clock
measures a single quadrature of $\phi$. For $MN$ atoms the projection-noise-
limited fractional frequency stability is
$\sigma_{yp}(\tau)=\sqrt{12\alpha}/(\omega\sqrt{MN\pi\tau})$.

It has been suggested Buzek1999Optimal ; Lloyd2013Private  that if $M$
oscillators with frequencies $\omega$, $\omega/2$,…,$\omega/2^{M-1}$ are
available, then durations up to $2^{M}\pi/\omega$ can be measured with the
precision of the fastest oscillator. Here we show that when $M$ atomic
ensembles are probed with reduced free-evolution periods or which evolve at
reduced frequencies, the range of unambiguous phase excursions is extended to
$\theta=D^{M-1}\pi$, where $D^{0},D^{1},\ldots,D^{M-1}$ are either probe-time
or frequency-division ratios. Then the Ramsey free-evolution period in the
presence of $1/f$ oscillator noise can grow to $T=D^{M-1}\pi/(6\alpha)$. The
projection-noise limited stability for the the $N$ unperturbed atoms is then
$\sigma_{yp}(\tau)=\sqrt{6\alpha}/(\omega\sqrt{D^{M-1}N\pi\tau})$. This
represents a reduction in the clock’s frequency variance to $MD^{1-M}/2$ times
the variance of a standard Ramsey clock with the same number of atoms.

The hypothetical clock contains $M$ ensembles of $N$ atoms. Atoms are thought
of as two-level systems or qubits. Each ensemble $A_{j}$ is divided into two
sub-ensembles of $N/2$ atoms, labeled $X_{j}$ and $Y_{j}$, where
$j\in{0...M-1}$. The atoms of all $X_{j}$ are simultaneously prepared in the
quantum ground state, and probed by a Ramsey sequence that consists of a
$\pi/2$-pulse about the Bloch-sphere x-axis, followed by evolution for a
period $T$. During this period, the atoms evolve through a reduced phase
$\phi_{j}=\int^{T}_{0}(\omega-\Omega(t))/D^{j}dt$ with respect to the
oscillator. A possible construction by use of dynamical decoupling is
described later. The evolution period ends with a second $\pi/2$-pulse about
the x-axis, and final measurement along the z-axis. The measurement result is
recorded as the fraction of atoms $x_{j}$ in the excited state. The $Y_{j}$
sub-ensembles are simultaneously probed in the same way, except that the
second $\pi/2$-pulse is about the y-axis, and their results are recorded as
$y_{j}$. The expectation values $\langle x_{j}\rangle=(1+\cos{\phi_{j}})/2$
and $\langle y_{j}\rangle=(1+\sin{\phi_{j}})/2$ are quadratures of $\phi_{j}$,
and the value $\beta_{j}=arg((x_{j}-1/2)+i(y_{j}-1/2))$ is a measurement of
$\phi_{j}$, modulo $2\pi$. Here $arg$ is the argument function, whose range is
$-\pi\ldots\pi$.

First consider the phase uncertainty of each ensemble. The measurements
$x_{j}$ and $y_{j}$ yield an estimate of the phase $\phi_{j}$. This can be
written as $\phi_{j}=\beta_{j}+2P_{j}\pi+G_{j}$, where $P_{j}$ is an integer
whose value may be unknown, and $G_{j}$ is a random variable that completes
the equation. The distribution of $G_{j}$ values is determined by the $arg$
function and the scaled binomial random variables $x_{j}$ and $y_{j}$. In the
limit of large $N$, $G_{j}$ approximates a Gaussian random variable with
variance $\sigma^{2}=1/N$, independently of the actual phases
WMI1993ProjectionNoise . To determine the ensemble size $N$ that is required
to reduce the probability of inversion errors to a certain level, we perform
precise numerical calculations and do not use the Gaussian approximation.

Figure 1: Depiction of how phase measurements can be combined for improved
clock stability. In (a) two ensembles with frequency ratio $D$ evolve
synchronously for period $T$, and their phase measurements are then combined
(see Eq. 1) to derive the integer number of $2\pi$ phase wraps $P_{j-1}$ of
the second ensemble. In (b) the first ensemble is measured $D$ times during
each measurement of the second ensemble. This information also yields the
number of phase wraps. Here the second index is an integer with range
$0\ldots(D^{j}-1)$ and refers to the sequence of measurements during one
period $T$.

Results from measuring the $M$ ensembles are combined to yield an estimate of
the total phase $\phi_{0}$. This process can be understood iteratively from
lowest frequency to highest, as shown in Fig. 1. Consider the ensembles
$A_{j}$ and $A_{j-1}$ where the integer $P_{j}$ is assumed to be known, and
$P_{j-1}$ is unknown, but can be derived from the earlier equations involving
$\phi_{j}$ as:

$P_{j-1}=DP_{j}+\frac{D\beta_{j}-\beta_{j-1}}{2\pi}+\frac{DG_{j}-G_{j-1}}{2\pi}$
(1)

In practice, $P_{j-1}$ will be calculated as
$r(DP_{j}+\frac{D\beta_{j}-\beta_{j-1}}{2\pi})$ where $r()$ rounds its
argument to the nearest integer. The probability of assigning an incorrect
value to $P_{j-1}$ is equal to the probability that $|DG_{j}-G_{j-1}|\geq\pi$,
which can be made arbitrarily small by choice of ensemble size $N$ (see Fig.
2). Therefore, if the integer $P_{j}$ is known, then the integer $P_{j-1}$ can
be derived from the measurements. By induction, it follows that the number of
phase wraps $P_{0}$ is known, because the evolution time was chosen such that
$P_{M-1}$ is known to be zero. If errors in the determination of $P_{j-1}$
from $P_{j}$ are equivalent to a six-sigma Gaussian events ($2\times 10^{-9}$
probability), then $N=46$ is the minimum ensemble size when $D=2$ (see Fig.
2). The total error probability during each period $T$ is equal to $M-1$ (the
number of applications of Eq. 1) times the probability shown in Fig 2.

Reduced-frequency clocks might be constructed from atoms or molecules with
lower-frequency resonances, where the oscillator is derived by frequency
division of the original oscillator SAD2000Comb ; Schnatz1996CohLink . In this
case, the frequency ratios between ensembles $\phi_{j-1}/\phi_{j}$ will be
given by the natural frequency ratios of the transitions, rather than the
equal ratios $D$ considered so far, but the previous arguments still apply,
and we can derive a clock variance that is reduced by the product of the
transition-frequency ratios.

A sequence of $\pi$-pulses, known as dynamical-decoupling, can also reduce the
rate of phase evolution. Optimized pulse sequences have been developed to
extend the coherence time of qubits, with the goal of making their phase
evolution with respect to the oscillator as close to zero as possible. The
same approach that rephases qubits for frequency drifts of polynomial order
$n$, by $n+1$ $\pi$-pulses Uhrig2007 , can reduce the phase evolution to
$\phi_{j}=\phi_{0}D^{-j}$. For example, the reduced phases can be generated in
the face of linear drift of the oscillator frequency
$\Omega(t)=\Omega_{0}+t\Omega_{1}$ during a period $T_{D}$ by one $\pi$-pulse
at $T_{A}=(T_{D}/4)(1+D^{-j})$, and a second $\pi$-pulse at
$T_{B}=(T_{D}/4)(3-D^{-j})$. In this case,
$\phi_{j}=\int_{0}^{T_{A}}(\omega-\Omega(t))dt-\int_{T_{A}}^{T_{B}}(\omega-\Omega(t))dt+\int_{T_{B}}^{T_{D}}(\omega-\Omega(t))dt=\int_{0}^{T_{D}}(\omega-\Omega(t))D^{-j}dt$,
as is required. Higher-order sequences can be formed by adapting the recipe of
ref. Uhrig2007 . The total evolution period $T$ can contain many decoupling
sequences of length $T_{D}$ to closely approximate the ideal value of
$\phi_{j}$. In real clocks, generating dynamical decoupling sequences with
sufficient fidelity is likely to be difficult.

As an alternative to the quadrature Ramsey measurements discussed above, one
can consider more general quantum measurements on each ensemble to derive the
phase estimates $\beta_{j}$. Optimal measurements will minimize the
probability of large phase errors. First prepare an initial state $\ket{a}$
and after evolution through a phase $\phi_{j}$, measure in some basis. Bužek
et al. Buzek1999Optimal  suggest a measurement basis spanned by the states
peres1980measurement  $\ket{k}=\frac{1}{\sqrt{N+1}}\sum_{m=0}^{N}e^{im2\pi
k/(N+1)}\ket{N,m}$ where $\ket{N,m}$ are the fully symmetrized states with $N$
qubits and $m$ excitation (e.g.
$\ket{3,2}=(\ket{011}+\ket{101}+\ket{110})/\sqrt{3}$). The probability
amplitude of measuring $\ket{k}$ is simply the $k$th coefficient of the
discrete Fourier transform of the values $\bra{a}e^{-im\phi_{j}}\ket{N,m}$,
where $m\in 0\ldots N$ (see Fig. 3). For the Gaussian initial states
Lukin2004Stability  where $\ket{a}=(-1)^{m}e^{-(m-N/2)^{2}/(Nc)}\ket{N,m}$
(unnormalized) the error rate is reduced compared to Ramsey measurements with
an equal number of qubits (see Fig. 2). Numerically, we find an optimal value
of $c=0.735$ and in this case $N=24$ is the smallest number of atoms that
yields an error rate below $2\times 10^{-9}$ when reduced frequency ensembles
with $D=2$ are combined. While they are not necessarily optimal, the Gaussian
states were chosen because the tail of a Gaussian function rapidly approaches
zero, for both the original function and its discrete Fourier transform. This
reduces the probability of large phase errors.

Figure 2: Probability of $2\pi$ phase errors $P$ as a function of ensemble
size $N$. Circles show the total probability that the term
$|(2G_{j}-G_{j-1})/(2\pi)|$ in Eq. 1 exceeds $1/2$, based on the binomial
distributions of $x_{j},y_{j},x_{j-1}$ and $y_{j-1}$, averaged uniformly over
the interval $-\pi\leq\phi_{j}\leq\pi$. The calculated probabilities are
approximated by the fit $P=0.07e^{-0.38N}$. Squares represent the probability
that $2\pi$ phase errors occur when reduced probe periods are used. In this
case the fit is $P=0.04e^{-0.48N}$. The error probabilities for Gaussian
initial states measured in the phase-state basis are shown as ($\times$) and
(+), and their fits are $P=0.05e^{-0.72N}$ and $P=0.25e^{-N}$ respectively.

Exponential scaling behavior is also possible if the atoms of all $M$
ensembles have frequency $\omega$, but each ensemble $A_{j}$ freely evolves
for a reduced period $T_{j}=D^{-j}T$ and is then measured and reprepared. In
this case, $D$ is an integer and measurement and preparation are assumed to be
instantaneous. Here $D$ measurements of duration $T_{j}$ from ensemble $A_{j}$
are combined to derive the integer number of phase wraps that occur during a
single measurement of ensemble $A_{j-1}$ (see Fig. 1). Measurements are
iteratively combined as before. Numerical evaluation shows that if $D=2$ and
$2\pi$ phase errors are again allowed to occur with $2\times 10^{-9}$
probability, then $N=36$ and $N=19$ are the minimum ensemble sizes for Ramsey
and Gaussian ($c=0.635$) measurements respectively (see Fig. 2). In this case,
the total error probability during each period $T$ is equal to $2^{M-1}-1$
times the probability shown in Fig 2.

Figure 3: Probabilities of $\ket{N,m}$ and $\ket{k}$ for phase measurement
with Gaussian states. Lines are drawn to guide the eye. The initial state
$\ket{a}$ for $N=20$ and $c=0.7$ in the $\ket{N,m}$ basis is represented by
(o). Symbols $\times$ \+ . $\square$ show the probabilities of measuring
$\ket{k}$ after phase evolution by $-2\delta,0,2\delta,$ and $(5/2)\delta$
respectively, where $\delta=2\pi/(N+1)$. Phase shifts correspond to
translation of $k$. Large phase measurement errors are confined to the low-
probability tails of the Gaussian function. For the discrete Fourier
transform, the product of Gaussian widths in $m$ and $k$ space is proportional
to $N$. The choice of $c$ is a compromise because the width in $k$ space
should be as narrow as possible to avoid large phase errors, but cutting off
the Gaussian in $m$ space adds unwanted Fourier components.

The above techniques show that oscillator phase noise in atomic clocks may be
mitigated by use of a collection of atomic ensembles. The free-evolution
period $T$ can, in principle, be extended far beyond the oscillator’s
decoherence limit, not only for differential frequency measurements
Chwalla2007Corr ; Olmshenk2007Corr ; CWC2011CorrSpec , but also for clocks
that keep time, as described here. This raises the question of what the
ultimate limits to atomic coherence are in clocks. In early work, the
millisecond transit time of particles in a thermal beam posed an artificial
limit to coherence Ramsey1956 . For microwave frequencies, this was extended
to about one second in atomic fountains Kasevich1989 , and 100 seconds in ion
traps Berkeland1998 .

Technical effects associated with fluctuating atomic frequencies pose one
limit to atomic coherence at optical frequencies and above ($\omega>2\pi\times
10^{14}$ Hz), but many of these can be reduced by experimental techniques.
Vacuum conditions pose another limit, because background-gas collisions can
randomize the atomic phase, or eject atoms from the trap. More fundamental is
the excited-state lifetime, although some optical transitions have very long-
lived excited states, such as the octupole transition in Yb+
Huntemann2012Octupole  and 1S${}_{0}-^{3}$P0 transitions in nuclear-spin zero
isotopes Taichenachev2006SpinZero .

Each ensemble or sub-ensemble may be constructed in a separate chamber to
avoid crosstalk, or a large collection of atoms can be subdivided into
separate ensembles. For the types of phase measurement considered here, the
minimum number requirement for performance gain while maintaining six-sigma
reliability is two ensembles of $36$ unentangled atoms or $19$ entangled
atoms. If the reduced-frequency clocks are generated by dynamical decoupling,
the fidelity of $\pi$-pulses is a critical concern Kirchmai2009HiFi .
Preparation and read-out speed may limit the performance of reduced period
clocks. These difficulties are avoided if ensembles with different natural
resonance frequencies can be combined. Technological barriers exist for all
cases, but it is possible that future clock stability can surpass the limits
of oscillator decoherence.

Helpful discussions with S. Lloyd, D. J. Wineland, M. D. Lukin, E. Knill, M.
Mullan, V. Vuletić and J. A. Sherman are gratefully acknowledged. This work is
supported by DARPA QuASAR, ONR, ARO, and AFOSR. Not subject to U.S. copyright.

## References

  * (1) C. M. Caves, Phys. Rev. D 23, 1693 (Apr 1981)
  * (2) V. Giovannetti, S. Lloyd, and L. Maccone, Science 306, 1330 (2004)
  * (3) W. M. Itano, J. C. Bergquist, J. J. Bollinger, J. M. Gilligan, D. J. Heinzen, F. L. Moore, M. G. Raizen, and D. J. Wineland, Phys. Rev. A 47, 3554 (1993)
  * (4) S. F. Huelga, C. Macchiavello, T. Pellizzari, A. K. Ekert, M. B. Plenio, and J. I. Cirac, Phys. Rev. Lett. 79, 3865 (Nov 1997)
  * (5) D. J. Wineland, C. Monroe, W. M. Itano, D. Leibfried, B. E. King, and D. M. Meekhof, J. Res. NIST 103, 259 (1998)
  * (6) V. Bužek, R. Derka, and S. Massar, Phys. Rev. Lett. 82, 2207 (1999)
  * (7) E. Riis and A. G. Sinclair, J. Phys. B 37, 4719 (2004)
  * (8) A. André, A. S. Sørensen, and M. D. Lukin, Phys. Rev. Lett. 92, 230801 (2004)
  * (9) T. Rosenband, in  _Proc. 2011 ICOLS Conf. ArXiv e-prints, quant-ph:1203.0288v2_
  * (10) M. Mullan and E. Knill, Quantum Info. Comput. 12, 553 (Jul. 2012), ISSN 1533-7146
  * (11) W. A. Kester,  _Data conversion handbook_ (Newnes, 2005) http://www.analog.com/library/analogDialogue/archives/39-06/data_conver%sion_handbook.html
  * (12) J. E. Volder, Electronic Computers, IRE Transactions on EC-8, 330 (1959), ISSN 0367-9950
  * (13) L. Vaidman and Z. Mitrani, Phys. Rev. Lett. 92, 217902 (May 2004)
  * (14) G. Giedke, J. M. Taylor, D. D’Alessandro, M. D. Lukin, and A. Imamoğlu, Phys. Rev. A 74, 032316 (Sep 2006)
  * (15) N. F. Ramsey,  _Molecular Beams_ (Oxford University Press, 1956)
  * (16) M. Kitagawa and M. Ueda, Phys. Rev. A 47, 5138 (1993)
  * (17) D. J. Wineland, J. J. Bollinger, W. M. Itano, and D. J. Heinzen, Phys. Rev. A 50, 67 (1994)
  * (18) I. D. Leroux, M. H. Schleier-Smith, and V. Vuletić, Phys. Rev. Lett. 104, 073602 (Feb 2010)
  * (19) J. J. Bollinger, W. M. Itano, D. J. Wineland, and D. J. Heinzen, Phys. Rev. A 54, R4649 (1996)
  * (20) R. Demkowicz-Dobrzański, Phys. Rev. A 83, 061802 (Jun 2011)
  * (21) K. Numata, A. Kemery, and J. Camp, Phys. Rev. Lett. 93, 250602 (2004)
  * (22) S. Lloyd, private communication, 2013
  * (23) S. A. Diddams, D. J. Jones, J. Ye, S. T. Cundiff, J. L. Hall, J. K. Ranka, R. S. Windeler, R. Holzwarth, T. Udem, and T. W. Hänsch, Phys. Rev. Lett. 84, 5102 (2000)
  * (24) H. Schnatz, B. Lipphardt, J. Helmcke, F. Riehle, and G. Zinner, Phys. Rev. Lett. 76, 18 (Jan 1996)
  * (25) G. S. Uhrig, Phys. Rev. Lett. 98, 100504 (Mar 2007)
  * (26) A. Peres, American Journal of Physics 48, 552 (1980)
  * (27) M. Chwalla, K. Kim, T. Monz, P. Schindler, M. Riebe, C. F. Roos, and R. Blatt, Appl. Phys. B 89, 483 (Dec. 2007)
  * (28) S. Olmschenk, K. C. Younge, D. L. Moehring, D. N. Matsukevich, P. Maunz, and C. Monroe, Phys. Rev. A 76, 052314 (Nov 2007)
  * (29) C. W. Chou, D. B. Hume, M. J. Thorpe, D. J. Wineland, and T. Rosenband, Phys. Rev. Lett. 106, 160801 (Apr 2011)
  * (30) M. A. Kasevich, E. Riis, S. Chu, and R. G. DeVoe, Phys. Rev. Lett. 63, 612 (Aug 1989)
  * (31) D. J. Berkeland, J. D. Miller, J. C. Bergquist, W. M. Itano, and D. J. Wineland, Phys. Rev. Lett. 80, 2089 (Mar 1998)
  * (32) N. Huntemann, M. Okhapkin, B. Lipphardt, S. Weyers, C. Tamm, and E. Peik, Phys. Rev. Lett. 108, 090801 (Feb 2012)
  * (33) A. V. Taichenachev, V. I. Yudin, C. W. Oates, C. W. Hoyt, Z. W. Barber, and L. Hollberg, Phys. Rev. Lett. 96, 083001 (Mar 2006)
  * (34) G. Kirchmair, J. Benhelm, F. Zähringer, R. Gerritsma, C. F. Roos, and R. Blatt, Phys. Rev. A 79, 020304 (Feb 2009)


# Site-selective doublon-holon dynamics in a pumped one-dimensional Hubbard
superlattice
with staggered Coulomb interactions

Zhenyu Cheng  Ying Li  College of Physics and Technology, Guangxi Normal
University, Guilin, Guangxi 541004, China  Hantao Lu  Lanzhou Center for
Theoretical Physics, Key Laboratory of Theoretical Physics of Gansu Province,
and Key Laboratory of Quantum Theory and Applications of MoE, Lanzhou
University, Lanzhou, Gansu 730000, China  Xiang Hu  College of Physics and
Technology, Guangxi Normal University, Guilin, Guangxi 541004, China
Zhongbing Huang  Department of Physics, Hubei University, Wuhan 430062, China
Gregory A. Fiete  Department of Physics, Northeastern University, Boston,
Massachusetts 02115, USA  Department of Physics, Massachusetts Institute of
Technology, Cambridge, Massachusetts 02139, USA  Liang Du  Liang Du:
<EMAIL_ADDRESS>College of Physics and Technology, Guangxi Normal
University, Guilin, Guangxi 541004, China

###### Abstract

Doublon-holon dynamics is investigated in a pumped one-dimensional Hubbard
model with a staggered on-site Coulomb interaction at half-filling. When the
system parameters are set to be in the Mott insulating regime the equilibrium
sublattice density of states exhibits several characteristic peaks,
corresponding to the lower and upper Hubbard bands as well as hybridization
bands. We study the linear absorption spectrum and find two main peaks
characterizing the photon frequencies which excite the ground state to an
excited state. For a system driven by a laser pulse with general intensity and
frequency, both the energy absorption and the doublon-holon dynamics exhibit
distinct behaviors as a function of laser amplitude and frequency. Single-
photon processes are observed at low laser intensity where the energy is
absorbed for resonance laser frequencies. For strong laser intensity multi-
photon induced dynamics are observed in the system, which are confirmed by an
evaluation of the Loschmidt amplitude. The contribution of multi-photon
processes to site-specific double occupancy is also characterized by the
generalized Loschmidt amplitude. The site-selective doublon-holon dynamics are
observed in both the one and multi-photon processes and the site-selective
behavior is explained within a quasiparticle picture. Our study suggests
strategies to optically engineer the doublon-holon dynamics in one dimensional
strongly correlated many-body systems.

**footnotetext: These authors contributed equally to this work

##  I INTRODUCTION

Non-equilibrium control of quantum states in strongly correlated systems with
optical techniques and its physical understanding have been attracting
attention from the condensed matter community in the past decade and remains
challenging[1, 2, 3]. Of particular research interests are: (1) thermalization
and prethermalization behavior as the system is driven far from equilibrium[1,
4]; (2) non-thermal states not accessible in equilibrium[5, 6, 3]; (3) non-
equilibrium control of quantum phase transitions [7, 2, 8, 9, 10, 11]; and (4)
disentangling degrees of freedom in the strongly correlated electronic
system[12].

In strongly correlated systems, the dynamics of the low-dimensional (one and
two spatial dimensions) Hubbard model has been studied extensively, with a
number of exotic behaviors observed. For example, in the pumped one-
dimensional Hubbard model, unconventional superconductivity[13], a photo-
induced insulator to metal transition[14], and charge, spin and $\eta-$spin
separation[15] have been revealed. Furthermore, in laser-pumped systems the
evolution of doubly occupied states (doublons) are dependent on the laser
frequency[16]. Such a non-equilibrium system can be described as a generalized
Gibbs ensemble[17]. On the other hand, in a two-dimensional Hubbard cluster,
the system exhibits Rabi-like oscillations[18], where the oscillation
frequency increases with the drive amplitude (intensity). Recently, d-wave
superconductivity has been observed in a driven two-dimensional plaquette
Hubbard model[19]. In contrast to the one-dimensional case, the doublon
dynamics in the Hubbard two-leg ladder exhibits weak coupling with magnetic
excitations[20]. In addition, impact ionization is observed in the two-
dimensional cluster[21, 22, 23, 24] where the double occupancy rises further
after the laser pulse.

The non-equilibrium studies summarized above are focused on systems with
spatial homogeneity in the Hamiltonian. In such systems, doublons and holons
are generated when a pump light is applied, which leads to an increase in the
total energy of the system. After the passing of the pump light through the
system, the system begins to relax through the recombination of doublons and
holons, where prethermalization behavior occurs[17]. By contrast, the non-
equilibrium dynamics in systems without spacial homogeneity (e.g., ionic
Hubbard model, Hubbard superlattice with a staggered Coulomb interaction) have
received relatively little attention. For the purpose of simplifying the
physical picture and the analysis in nonequilibrium studies of doulon-holon
dynamics, we focus our attention on the Hubbard superlattice with a staggered
Coulomb potential rather than the ionic Hubbard model where charge density
wave and spontaneous dimerized state[25] could occur and complicate the
analysis.

In this paper we focus on the non-equilibrium dynamics of the Hubbard
superlattice with spatially staggered Coulomb interactions, as shown in Fig.1.
In equilibrium such a Hubbard superlattice can be experimentally realized in
condensed matter with nano-scale spatial inhomogeneity[26, 27, 28, 29, 30,
31]. The simplest spatially inhomogeneous case is where there are two sites in
the unit cell with on-site Coulomb interactions $U_{a}$ for the A-sublattice
and $U_{b}$ for the B-sublattice. For a non-interacting B sublattice,
$U_{b}=0$, and finite $U_{a}$ in one-dimension, the system exhibits a
correlated-metallic phase at the particle–hole symmetrical point with zero
spin and zero charge gaps[30, 31]. For a B-sublattice with finite Coulomb
interaction strength $U_{b}>0$, the system is an antiferromagnetic Mott
insulator whis with zero spin gap and finite charge gap[30]. Besides
realizations in condensed matter systems, a spatial modulation of the
interaction has also been reported in 174Yb gas systems[32], and the system
undergoes a Mott metal-insulator transition as the Coulomb interaction
increases[33, 34, 35] in the infinite dimensional Bethe lattice, as shown
within the framework of dynamical mean-field theory[36].

In this work, by introducing a spatial inhomogenity with spatially alternating
(staggered) on-site Coulomb interactions, we study the doublon-holon dynamics
as a function of laser intensity and frequency in the one-dimensional Hubbard
superlattice with $U_{b}=18.0t_{h},U_{b}=3.0t_{h}$ (see Eq.(2)), where the
system is a Mott insulator before the pump is applied. For weak laser
intensity where linear response theory applies, we observe a site-selective
doublon-holon dynamics at laser frequency $\hbar\Omega\approx 3.2t_{h}$, where
the double occupancy of the B-sublattice is enhanced by the laser pulse, while
the double occupancy of A-sublattice remains almost unchanged during and after
the laser pulse. With strong laser intensity, multi-photon effects are
observed, and the site-selected doublon enhanced for laser frequency at
$\hbar\Omega=9.4t_{h}$, with site-A enhanced substantially while site-B
remains unchanged, in contrast to the selective behavior observed at
$\hbar\Omega=3.2t_{h}$. This observation opens the possibility through
appropriate laser protocol to optically engineer doublon-holon states in many-
body interacting systems.

Our paper is organized as follows. In Sec.II, we describe the Hamiltonian of
the pumped one-dimensional Hubbard model with modulated (staggered) site-
dependent Coulomb interaction and the time-dependent Lanczos method of
solution. In Sec.III, the equilibrium sublattice site-specific density of
states and the linear absorption spectrum are obtained using exact
diagonalization. In Sec.IV, we study the non-equilibrium dynamics of doublon-
holons as a function of laser frequency and amplitude, and we study the
eigenstate spectrum using the Loschmidt amplitude. In Sec.V, the generalized
Loschmidt amplitude is introduced to exhibit details of the site-specific
doublon-holon dynamics. Finally, in Sec.VI we present the main conclusions of
the paper.

Figure 1: (Color online) An example of a lattice with 10 sites, where the
onsite Coulomb interaction strength is modulated alternatively (staggered)
with $U_{a}\neq U_{b}$ for different sublattices.

##  II model and method

In equilibrium the one dimensional Hubbard model on the superlattice (A-B
sublattice) is written as [30, 31],

$\displaystyle H$
$\displaystyle=-t_{\mathrm{h}}\sum_{i\sigma}\left(c_{i\sigma}^{\dagger}c_{i+1\sigma}+h.c.\right)+\frac{\epsilon}{2}\sum_{i,\sigma}(n_{2i,\sigma}-n_{2i+1,\sigma})$
$\displaystyle+U_{a}\sum_{i}n_{2i+1,\uparrow}n_{2i+1\downarrow}+U_{b}\sum_{i}n_{2i,\uparrow}n_{2i\downarrow},$
(1)

where $c_{i\sigma}^{\dagger}$ ($c_{i\sigma}$) create (annihilate) a fermionic
particle at site $i$ with spin projection $\sigma$ and
$n_{i\sigma}=c_{i\sigma}^{\dagger}c_{i\sigma}$ for
$\sigma=\uparrow,\downarrow$. Here $t_{\mathrm{h}}$ is the anisotropic hopping
amplitude connecting nearest neighbour sites and $U_{a}(U_{b})$ is the onsite
Coulomb interaction between $\uparrow$ and $\downarrow$ spin electrons for odd
(even) sites in the one-dimensional chain. The energy $\mp\epsilon/2$ is the
onsite orbit energy for the A-B sublattice. The Hamiltonian is particle-hole
symmetric at $\epsilon=(U_{a}-U_{b})/2$ and can be rewritten as,

$\displaystyle H$
$\displaystyle=-t_{\mathrm{h}}\sum_{i\sigma}\left(c_{i\sigma}^{\dagger}c_{i+1\sigma}+h.c.\right)$
$\displaystyle+\sum_{i}U_{i}\left(n_{i\uparrow}-\frac{1}{2}\right)\left(n_{i\downarrow}-\frac{1}{2}\right),$
(2)

where $U_{i}$ is the site dependent Coulomb interaction, which we choose
$U_{i}=U_{a}(U_{b})$ for odd (even) sites in the one dimensional chain.

In this paper, the non-equilibrium calculations are done in the particle-hole
symmetric limit. We set $t_{\mathrm{h}}=1$ as the energy unit and
correspondingly, the unit of time is then the inverse of energy,
$t_{\mathrm{h}}^{-1}$. The site-specific density of states is defined as,

$\displaystyle\rho_{\alpha}(\omega)=\sum_{i\in\alpha,\sigma}\sum_{\phi}$
$\displaystyle|\langle\phi|c_{i\sigma}^{\dagger}|\psi_{0}\rangle|^{2}\delta(\omega-
E_{\phi}+E_{0})$ $\displaystyle+$
$\displaystyle|\langle\phi|c_{i\sigma}|\psi_{0}\rangle|^{2}\delta(\omega+E_{\phi}-E_{0})$
(3)

where $\\{|\phi\rangle\\}$ are the eigenstates of the equilibrium Hamiltonian
in Eq.(2) with respect to energy eigenvalues $E_{\phi}$, and
$|\psi_{0}\rangle$ is the ground state with energy $E_{0}$. The definition of
the linear absorption spectrum is [37, 13, 14],

$\displaystyle\alpha(\omega)=-\frac{1}{\pi}\mathrm{Im}\langle\psi_{0}|\hat{j}\frac{1}{\omega-(H-E_{0})}\hat{j}|\psi_{0}\rangle,$
(4)

and the current density operator is defined as,

$\displaystyle\hat{j}=it_{\mathrm{h}}\sum_{i\sigma}(c_{i\sigma}^{\dagger}c_{i+1,\sigma}-c_{i+1\sigma}^{\dagger}c_{i,\sigma}).$
(5)

The equation for linear absorption, Eq.(4), only works in low laser intensity
regime, where linear response theory is applicable.

We consider a system exposed to an external laser pulse with vector potential
(directed along the chain direction),

$\displaystyle
A(t)=A_{0}\exp[-(t-t_{p})^{2}/2t_{d}^{2}]\cos[\Omega(t-t_{p})],$ (6)

where $A_{0}$ is the laser intensity, $\Omega$ is the laser frequency and the
laser pulse is peaked at $t_{p}$ with $t_{d}$ characterizing the duration time
(pulse width) of light. The time-dependent Hamiltonian is written using the
Peierls substitution,

$\displaystyle H(t)$
$\displaystyle=-t_{\mathrm{h}}\sum_{i\sigma}\left(\exp[iA(t)]c_{i\sigma}^{\dagger}c_{i+1\sigma}+h.c.\right)$
$\displaystyle+\sum_{i}U_{i}\left(n_{i\uparrow}-\frac{1}{2}\right)\left(n_{i\downarrow}-\frac{1}{2}\right).$
(7)

We choose the chain size to be $L=10$ with periodic boundary conditions and a
coarse-grained time $\delta t=0.005t_{\mathrm{h}}^{-1}$. In the following, we
restrict ourselves to the case of half-filling with periodic boundary
conditions, where the total number of electrons $N$ is equal to the number of
sites in the chain $L$. Furthermore, we assume the total magnetization in the
system vanishes, which means the number of up-spin electrons $N_{\uparrow}$ is
equal to the down-spin electrons $N_{\downarrow}$.

The exact diagonalization method (a standard Lanczos procedure) is employed to
numerically calculate the ground state of the Hamiltonian at time $t=0^{-}$,
where the lase pulse is not yet applied to the system. The ground state is
used as an initial state for the time dependent Schrödinger equation
$i\partial_{t}|\Psi(t)\rangle=H(t)|\Psi(t)\rangle$. The time evolution is
implemented step-by-step based on the time-dependent Lanczos method [38, 39,
40, 41, 16],

$|\Psi(t+\delta t)\rangle\approx e^{-iH(t)\delta
t}|\Psi(t)\rangle\approx\sum_{l=1}^{M}e^{-i\epsilon_{l}\delta
t}|\Phi_{l}\rangle\langle\Phi_{l}|\Psi(t)\rangle,$ (8)

where $\epsilon_{l}$ ($\Phi_{l}$) are the eigenvalues (eigenvectors) of the
tri-diagonal matrix generated by Lanczos iteration with $M\leq 100$. (In
general, the desired M for achieving a given accuracy depends on the setups
and chosen time step sizes [42].) We set the time step size $\delta
t=0.005t_{h}^{-1}$ in our calculation of the time evolution. The physical
observable are computed as,

$\displaystyle\langle O(t)\rangle=\langle\Psi(t)|O|\Psi(t)\rangle.$ (9)

For the case of high laser intensity, multi-photon processes appear or even
dominate the optical excitation in the system. To validate the multi-photon
processes and provide more details on the double-holon dynamics, the Loschmidt
amplitude is used to calculate the spectral density[43, 22, 18],

$L(\omega,t)=\sum_{n}|\langle
E_{n}|\Psi(t)\rangle|^{2}\delta[\omega-(E_{n}-E_{g})]$ (10)

where $|E_{n}\rangle$ and $E_{n}$ are the eigenstates and eigenenergies
($E_{g}$ is the many-body ground state energy) of the unperturbed Hamiltonian
(after pulse Hamiltonian), and $|\Psi(t)\rangle$ is the time evolved state at
time $t$. Note that, the eigenstate spectrum is kept unchanged as a function
of time after the pump pulse has decayed, because $H(t)=H(0)$ and as a result
$|\langle E_{n}|\Psi(t)\rangle|$ is unchanged.

Figure 2: (Color online) For the Hubbard super-lattice with fixed site-
dependent Coulomb interaction strength $U_{b}=3.0$ and $U_{a}=18.0$, (a) The
sub-lattice specified density of states in equilibrium. (b) The linear
absorption spectrum as a function of frequency calculated using Eq.(4), where
linear response theory apply.

Figure 3: (Color online) For the Hubbard super-lattice with fixed site-
dependent Coulomb interaction strength $U_{b}=3.0$ and $U_{a}=18.0$, (a) The
after pulse ($t=15.0$ used) energy as a function of laser frequency and
amplitude. (b) the after pulse energy as a function of laser amplitude with
fixed laser frequency $\Omega=3.2,6.0,9.4,11.4$.

Figure 4: (Color online) For the Hubbard super-lattice with fixed site-
dependent Coulomb interaction strength $U_{b}=3.0$ and $U_{a}=18.0$, The time
evolution of total energy and site-specified double occupancy as a function of
frequency $\Omega$ are plotted in one column with laser intensity $A_{0}=0.1$
(a1-a3), $A_{0}=0.2$ (b1-b3), $A_{0}=0.6$ (c1-c3), $A_{0}=1.0$ (d1-d3).

Figure 5: (Color online) For the Hubbard super-lattice with fixed site-
dependent Coulomb interaction strength $U_{b}=3.0$ and $U_{a}=18.0$, the time
evolution of site dependent double occupancy with the Loschmidt amplitude in
the third row for different laser amplitude $A_{0}$ and frequency $\Omega$.
(a1-a3) $\Omega=3.2$, (b1-b3) $\Omega=6.0$ (c1-c3) $\Omega=9.4$ (d1-d3)
$\Omega=11.4$, respectively.

##  III Density of states and Linear absorption spectrum

In this section, we focus our attention on the Hubbard superlattice in the
Mott insulating phase with finite charge gap and zero spin gap[30], where the
system parameters are $U_{a}=18.0$ and $U_{b}=3.0$ (measured in units of
$t_{h}$). In Fig.2(a), we plot the sublattice-specified density of states in
Eq.(3) at zero temperature using the Lanczos method. Evidently, the particle-
hole symmetries are observed for both A and B sublattices, with each site
singly occupied (half-filling). The energy gap at the Fermi energy
($E_{F}=0.0$) confirms that the system is an insulator. Note the density of
states exhibits four peaks for each sublattice. For the A-sublattice, we
observe the lower and upper Hubbard band at approximately $\pm U_{a}/2=\pm
9.0$ and the two hybridization bands around $\pm U_{b}/2=\pm 1.5$. In
contrast, the upper and lower Hubbard band for the B-sublattice is situated at
$\pm U_{b}/2=\pm 1.5$ with tiny hybridization band at $\pm U_{a}/2=\pm 9.0$.

To study effect of a laser drive on the equilibrium system (taking it out-of-
equilibrium), we first study the energy absorption of the superlattice system
at low laser intensity where linear response theory applies. In Fig.2(b), the
linear absorption spectrum $\alpha(\omega)$ is calculated using Eq.(4), where
Fermi’s golden rule applies[24]. The positions of the peaks in the linear
absorption spectrum (one-photon excited states) are $\omega\approx 3.2,11.4$
(two peaks at $\omega\approx 11.0,11.5$ are merged into a single peak). The
peak at $\omega\approx 3.2$ corresponds to the energy difference between the
lower Hubbard band of the B-sublattice and the upper hybridization peak of the
A-sublattice, or inversely, from lower hybridization peak of A-sublattice to
the upper Hubbard band of B-sublattice. The peak at $\omega\approx 11.4$
corresponds to the energy difference between the lower Hubbard band of the
B-sublattice and the upper Hubbard band of the A-sublattice, or inversely,
from lower Hubbard band of A sublattice to the upper Hubbard band of
B-sublattice. The physical picture can be clearly understood as follows: Since
the hopping between the nearest neighboring A-B sites is dominate, that
process will contribute most to the optical excitation. By contrast, the next-
nearest neighbor A-A or B-B hopping, which is a second order process (compared
with the nearest neighbor A-B first order hopping process) is weaker.

##  IV Non-equilibrium driven behavior in the Mott insulating regime

To gain an overall picture of photon absorption in the system, we plot the
energy following the pulse ($t>15.0$) as a function of laser frequency
$\Omega$ and laser intensity $A_{0}$ in Fig.3(a). Note the energy remains
unchanged after the laser pulse has passed since the time evolution operator
is unitary (no energy is being added or removed from the system). A linear
response behavior is observed for small laser intensity for the one-photon
resonance frequencies $\Omega\approx 3.2,11.4$, which were clarified
previously by the linear absorption spectrum in Fig.2(b). To make the linear
response regime clear, we plot the energy as a function of laser intensity
with fixed laser frequency $\Omega=3.2,6.0,9.4,11.4$ in Fig.3(b). The total
energy exhibits linear energy absorption in the regime $A_{0}<0.30$ for laser
frequencies $\Omega=3.2,11.4$. In contrast, the total energy does not change
until $A_{0}\geq 0.2$ for laser frequencies $\Omega=6.0,9.4$, which suggest
single photon absorption is absent for a lower intensity laser and multi-
photon excitation occur for higher laser intensity.

To provide more details of the optical excitation process, we focus our
attention on the non-equilibrium evolution of total energy and double
occupancy of the system with the super-lattice (staggered Coulomb
interaction). The time-dependent energy as a function of laser frequency are
plotted in the top row of Fig.4(a1, b1, c1, d1), where the laser intensities
are set as $A_{0}=0.1,0.2,0.6,1.0$, respectively. The time evolution of the
site-dependent double occupancy are plotted in Fig.4 (a2-a3) for laser
intensity $A_{0}=0.1$, (b2-b3) for $A_{0}=0.2$ , (c2-c3) for $A_{0}=0.6$ and
(d2-d3) for $A_{0}=1.0$.

For $A_{0}=0.1$ in the linear response regime, the photon energy is absorbed
at frequency $\Omega=11.0$ and a small energy absorption is observed at
frequency $\Omega=3.2$, which are consistent with the equilibrium calculation
of the linear absorption spectrum in Fig.2(b). By inspecting the site-specific
double occupancy as a function of time and laser frequency, we find that an
increase of double occupancy on both site-A and site-B at frequency
$\Omega=11.4$ will contribute to the total energy absorption. While for laser
frequency $\Omega=3.2$, the energy absorption mainly comes from the
enhancement of double occupancy only at site B. The behaviors above can be
explained with the following physical picture: For laser frequency
$\Omega=11.4$, the single photon absorption induced optical excitation will
move electrons from the lower Hubbard band of the A sublattice to the upper
Hubbard band of the B sublattice, which enhances the double occupancy of the
B-site. Correspondingly, the photon can also excite electrons from the lower
Hubbard band of B sub-lattice to the upper Hubbard band of A sublattice,
enhancing the A-site double occupancy simultaneously. By contrast, the physics
is different for optical excitation with laser frequency $\Omega=3.2$. The
single photon process will excite electrons from the lower Hubbard band of B
to the upper hybridization band of A, enhancing the kinetic energy while not
affecting the double occupancy of the A-site. However, the optical excitation
from the lower hybridization A band to the upper Hubbard bands of the B
sublattice will enhance the double occupancy of the B site. As a result, the
system will exhibit site-selective doublon-holon dynamics at $\Omega=3.2$
which results from a single photon excitation. Further increasing the laser
intensity to $A_{0}=0.2$ will not introduce new physics, except it enhances
the effects observed with $A_{0}=0.1$, indicating that $A_{0}=0.2$ is still
situated in the linear response region of optical excitation.

Increasing the laser intensity to $A_{0}=0.6$, besides the enhancement of
energy absorption at $\Omega=3.2,11.4$ observed before, results in energy
absorption at at laser frequencies $\Omega=9.4,6.0,1.5$ respectively. The
extra absorption frequencies can be explained by multi-photon absorption
processes, which indicate $A_{0}=0.6$ is out of the linear response regime.
For $\Omega=6.0,1.5$, the double-photon processes correspond to the excitation
of single photon process of $\Omega=11.4,3.2$, respectively (as the
frequencies are roughly doubled). However, the two photon process for
$\Omega=9.4$ is different from the other two frequencies. For $\Omega=9.4$, an
apparent enhancement of double occupancy is observed for the A-site while not
for the B-site. The short time process is due to an excitation from the lower
Hubbard band of A sublattice to the upper Hubbard band of the B sublattice,
which enhances the double occupancy of B site. At subsequent times, the quasi-
particle is further excite to the upper Hubbard of the A sublattice, which
results in a reduction of the double occupancy of the B site and an increment
of the A site. Thus, the double occupancy of the A sites is enhanced while B
sites do not change in long time regime (which is induced by the two-photon
process in the relative long time regime). Further increasing the laser
intensity up to $A_{0}=1.0$ will enhance the observed doublon-holon dynamics
of $A_{0}=0.6$, and a triple photon absorption is observed at $\Omega=6.5$
(confirmed with Loschmidt amplitude, not shown). In summary, a clear frequency
selection behavior is observed from the energy absorption and site-specific
double occupancy enhancement, while the quasi-particle physical picture for
the selective dynamics can be different for different pump frequencies.

To make the laser frequency selection in the dynamical behavior more clear,
the site-specific double occupation as a function of time with different laser
intensities are plotted in Fig.5, with the Loschmidt amplitude (see also
Sec.V) as a function of $\omega$ plotted in the third row. For laser frequency
$\Omega=3.2$, the double occupancy of the A-site remains nearly unchanged
while the B-site is significantly enhanced with increased laser intensity. For
laser frequency $\Omega=8.0$ (not shown), the double occupancy of both the A-
and the B-site are unchanged with light. The physics can be explained using
the quasiparticle picture of photonic excitation in Ref. [22]. For laser
frequency $\Omega=11.40$, the double occupancy of both A-sites and B-sites are
strongly enhanced. For laser frequency $\Omega=6.0$, we find that the double
occupancy does not change for small laser intensity $A_{0}=0.05,0.10,0.20$,
while for $A_{0}=0.60,1.00$, the double occupancy of both A and B-site
increase in a similar fashion, which can be explained as multiphoton
processes. The multi-photon processes can be identified as double-photon
processes by looking at the Loschmidt amplitude, peaking at
$2\Omega,4\Omega,6\Omega$. Further, the oscillation frequency is independent
of the laser intensity, which is different from the Rabi-like behavior
(oscillation frequency increase with laser intensity) observed in the pumped
two-dimensional Hubbard cluster[18]. At the laser frequency $\Omega=9.4$, the
characterized two-photon process is a one-photon process from the lower
Hubbard band of A to the upper Hubbard band of B, with a subsequent excitation
to upper Hubbard band of A.

Figure 6: (Color online) For the Hubbard super-lattice with fixed site-
dependent Coulomb interaction strength $U_{b}=3.0$ and $U_{a}=18.0$, the
generalized Loschmidt amplitude as a function of energy for different laser
frequency, (a1-a2) $\Omega=3.2$, (b1-b2) $\Omega=6.0$ (c1-c2) $\Omega=9.4$
(d1-d2) $\Omega=11.4$, respectively. The green dash line is the Loschmidt
amplitude, which is a guidance of sum rule in Eq.(13) for eye.

##  V The Generalized Loschmidt Amplitude

In the context of site selective doublon-holon dynamics, it is interesting to
know which energies will contribute significantly to the double occupancy
enhancement of the A-site or the B-site after the laser pulse has passed
($t=40.0$ used here). To this end, we followed previous work[24] and defined
the generalized Loschmidt amplitude with respect to energy and site-specified
double occupancy,

$\displaystyle L_{\hat{H}\\!\hat{D}_{\alpha}}(\omega,D,t)=$
$\displaystyle\sum_{m,n}\langle\psi(t)|E_{n}\rangle\langle
E_{n}|D_{m}\rangle\langle D_{m}|\psi(t)\rangle$
$\displaystyle\times\delta(\omega-E_{n})\delta(D-D_{m}^{\alpha}).$ (11)

where
$\hat{D}_{\alpha}=\sum_{i\in\alpha}\hat{n}_{i\uparrow}\hat{n}_{i\downarrow}$
($\alpha=a,b$ for A, B sub-lattice) and
$\hat{D}_{\alpha}|D_{m}\rangle=D_{m}^{\alpha}|D_{m}\rangle$ with
$|D_{m}\rangle$ the many-body Fock basis with double occupancy
$D_{m}=D_{m}^{a}+D_{m}^{b}$. Here $|E_{n}\rangle$ are energy eigenstates of
the Hamiltonian with eigenvalue $E_{n}$. Due to the fact that the
eigenspectrum and eigenstates of double occupancy are a priori known, a
computational cheaper way to obtain the same information as Eq.(11) is using
projectors[24],

$\displaystyle
L_{H}(\omega,\hat{P}_{m})=\sum_{n}\langle\psi(t)|E_{n}\rangle\langle
E_{n}|\hat{P}_{m}|\psi(t)\rangle\delta(\omega-E_{n})$ (12)

where $\hat{P}_{m}=|D_{m}\rangle\langle D_{m}|$ is the projector operator onto
states with double occupancy, $D_{m}\in\\{0,1,2,\cdots,N/2\\}$,

$\displaystyle\sum_{D_{m}=0}^{N/2}L_{H}(\omega,\hat{P}_{m})=L(\omega).$ (13)

The total double occupancy is expressed as,

$\displaystyle\sum_{D_{m}=0}^{N/2}D_{m}\int d\omega
L_{H}(\omega,\hat{P}_{m})=\langle\hat{D}\rangle.$ (14)

The site-specific double occupancy is expressed as,

$\displaystyle\sum_{D_{m}=0}^{N/2}D_{m}^{\alpha}\int d\omega
L_{H}(\omega,\hat{P}_{m})=\langle\hat{D}_{\alpha}\rangle,$ (15)

where $\alpha=a,b$.

In Fig.6, we plot the generalized Loschmidt amplitude as a function laser
frequency, where the laser amplitude is fixed at $A_{0}=0.60$. For laser
frequency $\Omega=3.2$, the single and double photon processes will induce a
hopping within the $D_{a}=0$ subspace, while the single photon process will
include the state with one double occupancy in the B sublattice, and two
photon processes will induce two double occupancies on the B sites. As a
result, the A site double occupancy do not change and the B site double
occupancy is enhanced.

For laser frequency $\Omega=6.0$, the single photon processes are absent and
only double-photon process will enhance the $D_{a}=1$ subspace. Note that the
four photon process is a subsequent two-photon process (one following the
other). For laser frequency $\Omega=9.4$, the single photon processes are
absent and only a double-photon process will enhance the $D_{b}=1$ subspace,
and the the $D_{a}=0$ subspace. Note that the four photon process is a
subsequent two-photon process (one following the other). For laser frequency
$\Omega=11.4$, the single photon processes are observed and a single-photon
will enhance the $D_{b}=1$ subspace, and the $D_{a}=1$ subspace. Note that the
four photon process is a subsequent two-photon process (one following the
other).

##  VI Conclusion and Discussion

In this paper, we study the doubon-holon dynamics in the Hubbard superlattice
with alternatively spatially modulated (staggered) on-site Coulomb interaction
driven by a laser pulse. The Coulomb interaction parameters are set as
$U_{a}=18.0$ and $U_{b}=3.0$ for odd and even sites of the one-dimensional
chain with periodic boundary conditions adopted. In equilibrium, the site-
specific density of states is studied and shown to exhibit four peaks with
lower and upper Hubbard bands and two hybridization peaks. Within linear
response theory, we find two main resonance frequency–$\Omega=3.2$ and
$\Omega=11.4$–by studying the linear absorption spectrum. Focusing on the
dynamics of site-specific double occupancy, we study the time evolution of
double occupancy at the A and B sites, in the plane of laser intensity and
frequency. For small laser intensity, the system is in the linear response
theory regime, where the resonance frequencies are $\Omega=3.2$ and
$\Omega=11.4$. At $\Omega=3.2$, the double occupancy of the A-sites remain
approximately unchanged, while the double occupancy of the B-sites is largely
enhanced. The site-selective doublon-holon dynamics observed above is due to
the fact that the upper hybridization of the A band is mainly singly occupied,
while the upper Hubbard band of B is mainly doubly occupied. At $\Omega=11.4$,
the double occupancy of both A-sites and B-sites are increased simultaneously
by the laser, which are explained as the resonance excitation from the lower
Hubbard of A(B) to upper the Hubbard band of B(A). For higher laser intensity
beyond the linear response regime, multi-photon processes will play an
important role in the doublon-holon dynamics. At $\Omega=6.0$, the double
occupancy of the A- and B-sites remain almost unchanged (oscillating around
their equilibrium value at time $t=0.0$) for laser intensity
$A_{0}=0.05,0.10,0.20$, while for laser intensity $A_{0}=0.6,1.0$, the double
occupancy of both A, B-sites increase simultaneously. We conclude the double
photon processes at $\Omega=6.0$ for high laser intensity will induced hopping
events observed for single photon processes around $\Omega=11.4$. Furthermore,
at $\Omega=9.4$, the double occupancy of the A-sites are increased
significantly while that of the B-sites remain almost unchanged, which is what
is observed in site-selective doublon-holon dynamics. The physical picture for
$\Omega=9.4$ is different from the one observed at $\Omega=3.2$, where one
hopping event is introduced to explain the dynamics. The phenomena for
$\Omega=3.2$ is explained with two subsequent hopping events: the hopping from
the lower Hubbard band of A to the upper Hubbard band of B, followed by a
second hopping from the upper Hubbard band of B to the upper Hubbard band of
A. The physical pictures described above have been confirmed by direct
evaluation of the Loschmidt amplitudes. Our theoretical study suggests
strategies to engineer the doublon-holon dynamics with specified laser
parameters.

##  Appendix A Particle-Hole symmetry of the Hubbard superlattice

Introducing the particle-hole transformation in a bipartite lattice (A-B
sublattice),

$\displaystyle c_{i\sigma}^{\dagger}\rightarrow\eta_{i}d_{i\sigma},\quad
c_{i\sigma}\rightarrow\eta_{i}d_{i\sigma}^{\dagger},$ (A1)

where

$\displaystyle\eta_{i}=\left\\{\begin{array}[]{rl}1&\mathrm{if}\
i\in\mathrm{A},\\\ -1&\mathrm{if}\ i\in\mathrm{B}.\end{array}\right.$ (A4)

Consider the local terms in one unit cell,

$\displaystyle H_{\mathrm{loc}}$
$\displaystyle=U_{a}n_{a\uparrow}n_{a\downarrow}+U_{b}n_{b\uparrow}n_{b\downarrow}$
$\displaystyle+\epsilon(n_{a}-n_{b})-\mu(n_{a}+n_{b}).$ (A5)

The particle-hole transformed Hamiltonian is,

$\displaystyle\tilde{H}_{\mathrm{loc}}$
$\displaystyle=U_{a}(1-n_{a\uparrow})(1-n_{a\downarrow})+U_{b}(1-n_{b\uparrow})(1-n_{b\downarrow})$
$\displaystyle+\epsilon(n_{b}-n_{a})-\mu(4-n_{a}-n_{b}).$ (A6)

The difference between $H_{\mathrm{loc}}$ and its particle-hole transformation
$\tilde{H}_{\mathrm{loc}}$ is,

$\displaystyle H_{\mathrm{loc}}-\tilde{H}_{\mathrm{loc}}$
$\displaystyle=(U_{a}+2\epsilon-2\mu)n_{a}+(U_{b}-2\epsilon-2\mu)n_{b}$
$\displaystyle+(4\mu-U_{a}-U_{b}).$ (A7)

For the particle-hole symmetry, we have,

$\displaystyle\epsilon=(U_{a}-U_{b})/4,\quad\mu=(U_{a}+U_{b})/4.$ (A8)

##  Appendix B Position of lower and upper Hubbard band of the Hubbard
superlattice

In the homogeneous one band Hubbard model at half-filling, the position of
lower and upper Hubbard bands are situated at $-U/2$ and $U/2$ in equilibrium,
which can be understood by calculating the lesser Green’s function and greater
Green’s function at strong Coulomb interaction limit. The lesser Green’s
function is defined as,

$\displaystyle
G_{\uparrow}^{<}(\omega)=-\sum_{n}\frac{\langle\Psi_{g}|c_{\uparrow}^{\dagger}|n\rangle\langle
n|c_{\uparrow}|\Psi_{g}\rangle}{\omega+i\eta-(E_{0}-E_{n})},$ (B1)

where $E_{0}$ and
$|\Psi_{g}\rangle=|\uparrow,\downarrow,\uparrow,\downarrow,\cdots\rangle$ are
the ground state energy and vector in the subspace
$(N_{\uparrow},N_{\downarrow})$, $E_{n}$ and $|n\rangle$ are the eigen-energy
and eigen-vector in the subspace $(N_{\uparrow}-1,N_{\downarrow})$. As a
result, the $E_{0}-E_{n}$ constitute the position of lower Hubbard band. The
greater Green’s function is defined as,

$\displaystyle
G_{\uparrow}^{>}(\omega)=+\sum_{m}\frac{\langle\Psi_{g}|c_{\uparrow}|m\rangle\langle
m|c_{\uparrow}^{\dagger}|\Psi_{g}\rangle}{\omega+i\eta-(E_{m}-E_{0})},$ (B2)

where $E_{0}$ and $|\Psi_{g}\rangle$ are the ground state energy and vector in
the subspace $(N_{\uparrow},N_{\downarrow})$, $E_{m}$ and $|m\rangle$ are the
eigenenergy and eigenvector in the subspace $(N_{\uparrow}+1,N_{\downarrow})$.
As a result, the energy difference $E_{m}-E_{0}$ constitutes the position of
the upper Hubbard band with one more doubly occupied site.

For the half-filled Hubbard superlattice with staggered Coulomb interactions,
its atomic limit Hamiltonian in one unit cell is written as,

$\displaystyle H_{\mathrm{loc}}$
$\displaystyle=U_{a}\left(n_{a\uparrow}-1/2\right)\left(n_{a\downarrow}-1/2\right)$
$\displaystyle+U_{b}\left(n_{b\uparrow}-1/2\right)\left(n_{b\downarrow}-1/2\right).$
(B3)

The ground state $|\uparrow_{a},\downarrow_{b}\rangle$
($|\downarrow_{a},\uparrow_{b}\rangle$) energy in
$(N_{\uparrow}=1,N_{\downarrow}=1)$ subspace is $-(U_{a}+U_{b})/4$. For the
A-sublattice (where the newly added or deleted electron is from the A site),
in subspace $(N_{\uparrow}=0,N_{\downarrow}=1)$, the ground state energy is
$(U_{a}-U_{b})/4$, which conclude the lower Hubbard band is situated at
$-U_{a}/2$. In subspace $(N_{\uparrow}=2,N_{\downarrow}=1)$, the ground state
energy is $(U_{a}-U_{b})/4$, where the upper Hubbard band is situated at
$U_{a}/2$. For the B-sublattice (where the newly added or deleted electron is
from the B-site), at subspace $(N_{\uparrow}=1,N_{\downarrow}=0)$, the ground
state energy is $(U_{b}-U_{a})/4$, where the lower Hubbard band is situated at
$-U_{b}/2$. In subspace $(N_{\uparrow}=1,N_{\downarrow}=2)$, the ground state
energy is $(U_{b}-U_{a})/4$, where the upper Hubbard band is situated at
$U_{b}/2$.

Figure A1: (Color online) For the Hubbard super-lattice with fixed site-
dependent Coulomb interaction strength $U_{b}=3.0$ and $U_{a}=18.0$, the time
evolution of site dependent double occupancy. (a-b) $\Omega=3.2$ A-site and
B-site, (c-d) $\Omega=11.0$ A-site and B-site, respectively.

##  Appendix C Finite size effects

To determine the finite size effects in the system, we compare our numerical
results of the pumped one-dimensional chain with 32 sites, solved with density
matrix renormalization group (DMRG) theory using Itensor[44, 45]. The system
parameters are $U_{a}=18.0,U_{b}=3.0$. The laser frequency is set as
$\Omega=3.2$ in Fig.A1 (a-b) and $\Omega=11.0$ in Fig.A1 (c-d). The double
occupancy calculated with DMRG and exact diagonalization show good agreement
with each other, which indicate the finite size effects in this work are not
severe.

###### Acknowledgements.

We acknowledge helpful discussions with Tao Li, Huan Li, Xiaojun Zheng. We
gratefully acknowledge funding from the Natural Science Foundation of Guangxi
Province Grant No. 2020GXNSFAA297083, GuiKe AD20297045 and the National
Natural Science Foundation of China (Grant No. 11904143, No. 12174168, No.
12047501, No. 12147211, and No. 12065002), G.A.F. gratefully acknowledges
funding from US National Science Foundation grant DMR-2114825 and the
Alexander von Humboldt Foundation.

## References

  * Aoki  _et al._ [2014] H. Aoki, N. Tsuji, M. Eckstein, M. Kollar, T. Oka, and P. Werner, Nonequilibrium dynamical mean-field theory and its applications, Reviews of Modern Physics 86, 779 (2014).
  * Giannetti  _et al._ [2016] C. Giannetti, M. Capone, D. Fausti, M. Fabrizio, F. Parmigiani, and D. Mihailovic, Ultrafast optical spectroscopy of strongly correlated materials and high-temperature superconductors: A non-equilibrium approach, Advances in Physics 65, 58 (2016).
  * Heyl [2018] M. Heyl, Dynamical quantum phase transitions: A review, Rep. Prog. Phys. 81, 054001 (2018).
  * Mori  _et al._ [2018] T. Mori, T. N. Ikeda, E. Kaminishi, and M. Ueda, Thermalization and prethermalization in isolated quantum systems: A theoretical overview, J. Phys. B: At. Mol. Opt. Phys. 51, 112001 (2018).
  * Berges  _et al._ [2008] J. Berges, A. Rothkopf, and J. Schmidt, Nonthermal Fixed Points: Effective Weak Coupling for Strongly Correlated Systems Far from Equilibrium, Phys. Rev. Lett. 101, 041603 (2008).
  * Chiocchetta  _et al._ [2017] A. Chiocchetta, A. Gambassi, S. Diehl, and J. Marino, Dynamical Crossovers in Prethermal Critical States, Phys. Rev. Lett. 118, 135701 (2017).
  * Fausti  _et al._ [2011] D. Fausti, R. I. Tobey, N. Dean, S. Kaiser, A. Dienst, M. C. Hoffmann, S. Pyon, T. Takayama, H. Takagi, and A. Cavalleri, Light-Induced Superconductivity in a Stripe-Ordered Cuprate, Science 331, 189 (2011).
  * Cavalleri [2018] A. Cavalleri, Photo-induced superconductivity, Contemporary Physics 59, 31 (2018).
  * De La Torre  _et al._ [2021] A. De La Torre, D. M. Kennes, M. Claassen, S. Gerber, J. W. McIver, and M. A. Sentef,  _Colloquium:_ Nonthermal pathways to ultrafast control in quantum materials, Rev. Mod. Phys. 93, 041002 (2021).
  * Budden  _et al._ [2021] M. Budden, T. Gebert, M. Buzzi, G. Jotzu, E. Wang, T. Matsuyama, G. Meier, Y. Laplace, D. Pontiroli, M. Riccò, F. Schlawin, D. Jaksch, and A. Cavalleri, Evidence for metastable photo-induced superconductivity in K3C60, Nat. Phys. 17, 611 (2021).
  * Rowe  _et al._ [2023] E. Rowe, B. Yuan, M. Buzzi, G. Jotzu, Y. Zhu, M. Fechner, M. Först, B. Liu, D. Pontiroli, M. Riccò, and A. Cavalleri, Resonant enhancement of photo-induced superconductivity in K3C60, Nat. Phys. 10.1038/s41567-023-02235-9 (2023).
  * Dal Conte  _et al._ [2012] S. Dal Conte, C. Giannetti, G. Coslovich, F. Cilento, D. Bossini, T. Abebaw, F. Banfi, G. Ferrini, H. Eisaki, M. Greven, A. Damascelli, D. Van Der Marel, and F. Parmigiani, Disentangling the Electronic and Phononic Glue in a High- _T_ c Superconductor, Science 335, 1600 (2012).
  * Kaneko  _et al._ [2019] T. Kaneko, T. Shirakawa, S. Sorella, and S. Yunoki, Photoinduced $\eta$ Pairing in the Hubbard Model, Physical Review Letters 122, 077002 (2019).
  * Ejima  _et al._ [2022] S. Ejima, F. Lange, and H. Fehske, Nonequilibrium dynamics in pumped Mott insulators, Physical Review Research 4, L012012 (2022).
  * Murakami  _et al._ [2023] Y. Murakami, S. Takayoshi, T. Kaneko, A. M. Läuchli, and P. Werner, Spin, Charge, and $\eta$ -Spin Separation in One-Dimensional Photodoped Mott Insulators, Phys. Rev. Lett. 130, 106501 (2023).
  * Innerberger  _et al._ [2020] M. Innerberger, P. Worm, P. Prauhart, and A. Kauch, Electron-light interaction in nonequilibrium: exact diagonalization for time-dependent Hubbard Hamiltonians, European Physical Journal Plus 135, 10.1140/epjp/s13360-020-00919-2 (2020).
  * Murakami  _et al._ [2022] Y. Murakami, S. Takayoshi, T. Kaneko, Z. Sun, D. Golež, A. J. Millis, and P. Werner, Exploring nonequilibrium phases of photo-doped Mott insulators with generalized Gibbs ensembles, Commun Phys 5, 1 (2022).
  * Okamoto and Peronaci [2021] J. Okamoto and F. Peronaci, Floquet prethermalization and Rabi oscillations in optically excited Hubbard clusters, Sci Rep 11, 17994 (2021).
  * Zhang  _et al._ [2023] Y. Zhang, R. Mondaini, and R. T. Scalettar, Photoinduced enhancement of superconductivity in the plaquette Hubbard model, Phys. Rev. B 107, 064309 (2023).
  * Dias Da Silva  _et al._ [2012] L. G. G. V. Dias Da Silva, G. Alvarez, and E. Dagotto, Dynamics of doublon-holon pairs in Hubbard two-leg ladders, Phys. Rev. B 86, 195103 (2012).
  * Kauch  _et al._ [2020] A. Kauch, P. Worm, P. Prauhart, M. Innerberger, C. Watzenböck, and K. Held, Enhancement of impact ionization in Hubbard clusters by disorder and next-nearest-neighbor hopping, Phys. Rev. B 102, 245125 (2020).
  * Maislinger and Evertz [2022] F. Maislinger and H. G. Evertz, Impact ionization and multiple photon absorption in the two-dimensional photoexcited Hubbard model, Phys. Rev. B 105, 045114 (2022).
  * Gazzaneo  _et al._ [2022] P. Gazzaneo, T. M. Mazzocchi, J. Lotze, and E. Arrigoni, Impact ionization processes in a photodriven Mott insulator: Influence of phononic dissipation, Phys. Rev. B 106, 195140 (2022).
  * Watzenböck  _et al._ [2022] C. Watzenböck, M. Wallerberger, L. Ruzicka, P. Worm, K. Held, and A. Kauch, Photoexcitations in the Hubbard model: Generalized Loschmidt amplitude analysis of impact ionization in small clusters, Phys. Rev. B 106, 085135 (2022).
  * Fabrizio  _et al._ [1999] M. Fabrizio, A. O. Gogolin, and A. A. Nersesyan, From Band Insulator to Mott Insulator in One Dimension, Phys. Rev. Lett. 83, 2014 (1999).
  * Paiva and Dos Santos [1996] T. Paiva and R. R. Dos Santos, Electronic Correlations in One-Dimensional Superlattices, Phys. Rev. Lett. 76, 1126 (1996).
  * Paiva and Dos Santos [1998] T. Paiva and R. R. Dos Santos, Metal-insulator transition in one-dimensional Hubbard superlattices, Phys. Rev. B 58, 9607 (1998).
  * Paiva and Dos Santos [2000] T. Paiva and R. R. Dos Santos, Magnetism in one-dimensional Hubbard superlattices, Phys. Rev. B 62, 7007 (2000).
  * Paiva and Dos Santos [2002] T. Paiva and R. R. Dos Santos, Charge-density waves in one-dimensional Hubbard superlattices, Phys. Rev. B 65, 153101 (2002).
  * Duan and Wang [2010] C.-B. Duan and W.-Z. Wang, Incommensurate charge correlation and phase diagram of the one-dimensional superlattice Hubbard model at half-filling, Journal of Physics: Condensed Matter 22, 345601 (2010).
  * Zhang  _et al._ [2015] L.-L. Zhang, J. Huang, C.-B. Duan, and W.-Z. Wang, Structure-dependent metal–insulator transition in one-dimensional Hubbard superlattice, Chinese Phys. B 24, 077101 (2015).
  * Yamazaki  _et al._ [2010] R. Yamazaki, S. Taie, S. Sugawa, and Y. Takahashi, Submicron Spatial Modulation of an Interatomic Interaction in a Bose-Einstein Condensate, Phys. Rev. Lett. 105, 050405 (2010).
  * Saitou  _et al._ [2013] T. Saitou, A. Koga, and A. Yamamoto, Metal–Insulator Transition in Optical Lattice System with Site-Dependent Interactions, J Supercond Nov Magn 26, 1771 (2013).
  * Koga  _et al._ [2013] A. Koga, T. Saitou, and A. Yamamoto, Mott Transitions in the Hubbard Model with Spatially Modulated Interactions, J. Phys. Soc. Jpn. 82, 024401 (2013).
  * Hoang  _et al._ [2016] A.-T. Hoang, T.-H.-Y. Nguyen, T.-T.-T. Tran, and D.-A. Le, Two-component Fermions in Optical Lattice with Spatially Alternating Interactions, J. Phys. Soc. Jpn. 85, 104702 (2016).
  * Georges  _et al._ [1996] A. Georges, G. Kotliar, W. Krauth, and M. J. Rozenberg, Dynamical mean-field theory of strongly correlated fermion systems and the limit of infinite dimensions, Rev. Mod. Phys. 68, 13 (1996).
  * Okamoto [2019] J. Okamoto, Time-dependent spectral properties of a photoexcited one-dimensional ionic Hubbard model: an exact diagonalization study, New Journal of Physics 21, 123040 (2019).
  * Park and Light [1986] T. J. Park and J. C. Light, Unitary quantum time evolution by iterative Lanczos reduction, J. Chem. Phys. 85, 5870 (1986).
  * Mohankumar and Auerbach [2006] N. Mohankumar and S. M. Auerbach, On time-step bounds in unitary quantum evolution using the Lanczos method, Comput. Phys. Commun. 175, 473 (2006).
  * Balzer  _et al._ [2012] M. Balzer, N. Gdaniec, and M. Potthoff, Krylov-space approach to the equilibrium and nonequilibrium single-particle Green’s function, Journal of Physics: Condensed Matter 24, 035603 (2012).
  * Lu  _et al._ [2012] H. Lu, S. Sota, H. Matsueda, J. Bonča, and T. Tohyama, Enhanced Charge Order in a Photoexcited One-Dimensional Strongly Correlated System, Phys. Rev. Lett. 109, 197401 (2012).
  * Moler and Van Loan [2003] C. Moler and C. Van Loan, Nineteen Dubious Ways to Compute the Exponential of a Matrix, Twenty-Five Years Later, SIAM Review 45, 3 (2003).
  * Kennes  _et al._ [2020] D. M. Kennes, C. Karrasch, and A. J. Millis, Loschmidt-amplitude wave function spectroscopy and the physics of dynamically driven phase transitions, Phys. Rev. B 101, 081106 (2020).
  * Fishman  _et al._ [2022] M. Fishman, S. R. White, and E. M. Stoudenmire, The ITensor Software Library for Tensor Network Calculations, SciPost Phys. Codebases , 4 (2022).
  * Yang and White [2020] M. Yang and S. R. White, Time-dependent variational principle with ancillary Krylov subspace, Phys. Rev. B 102, 094315 (2020).