Optimal Hamiltonian simulation for time-periodic
systems
Kaoru Mizuta1 and Keisuke Fujii2,3,1,4
1

RIKEN Center for Quantum Computing (RQC), Hirosawa 2-1, Wako, Saitama 351-0198, Japan

2

Graduate School of Engineering Science, Osaka University, 1-3 Machikaneyama, Toyonaka, Osaka 560-8531,

arXiv:2209.05048v3 [quant-ph] 23 Mar 2023

Japan.
3

Center for Quantum Information and Quantum Biology, Osaka University, Japan.

4

Fujitsu Quantum Computing Joint Research Division at QIQB, Osaka University, 1-2 Machikaneyama, Toyonaka 560-0043, Japan

The implementation of time-evolution operators U (t), called Hamiltonian
simulation, is one of the most promising usage of quantum computers. For
time-independent Hamiltonians, qubitization has recently established efficient
realization of time-evolution U (t) = e−iHt , with achieving the optimal computational resource both in time t and an allowable error ε. In contrast,
those for time-dependent systems require larger cost due to the difficultyof handling time-dependency. In this paper, we establish optimal/nearlyoptimal Hamiltonian simulation for generic time-dependent systems with timeperiodicity, known as Floquet systems. By using a so-called Floquet-Hilbert
space equipped with auxiliary states labeling Fourier indices, we develop a
way to certainly obtain the target time-evolved state without relying on either time-ordered product or Dyson-series expansion. Consequently, the querycomplexity, which measures the cost for implementing the time-evolution, has
optimal and nearly-optimal dependency respectively in time t and inverse errorε, and becomes sufficiently close to that of qubitization. Thus, our protocoltells us that, among generic time-dependent systems, time-periodic systems
provides a class accessible as efficiently as time-independent systems despite
the existence of time-dependency. As we also provide applications to simulation of nonequilibrium phenomena and adiabatic state preparation, our results
will shed light on nonequilibrium phenomena in condensed matter physics andquantum chemistry, and quantum tasks yielding time-dependency in quantum
computation.

1 Introduction
Simulating quantum many-body systems is one of the most promising usages of quantum
computers that can fully exploit their computational powers, as R. Feynman’s proposal in
the beginning of quantum computers [1]. The implementation of a time-evolution operator
under Hamiltonians, called Hamiltonian simulation, has been the most important problem
for efficiently and accurately accomplishing quantum simulation. In fact, the applicationof Hamiltonian simulation nowadays ranges from condensed matter physics to quantum
Kaoru Mizuta: kaoru.mizuta@riken.jp

Accepted in Quantum 2023-03-16, click title to verify. Published under CC-BY 4.0.

1

chemistry, such as quantum dynamics [2, 3, 4] and quantum phase estimation [5, 6, 7, 8].
For long time, Trotterization has been a standard way of Hamiltonian simulation forboth time-independent and time-dependent systems, which can provide simple realization
available in today’s quantum computers [2, 9, 10, 11, 12]. Instead, it requires a huge
number of elementary gates up to poly (1/ε) to achieve an acceptable error ε.
For time-independent systems, various efficient Hamiltonian simulation algorithms
have appeared in the past decade, which yield fewer resources than Trotterization to
implement the time-evolution U (t) = e−iHt [13, 14, 15, 16, 17, 18]. Among them, thequbitization technique [18] achieves the best cost in that its query complexity (a measureof cost) has optimal scaling both in time t and inverse error 1/ε in an additive way. [18].
In contrast, for time-dependent systems, while there exist several efficient Hamiltonian
simulation for constructing time-evolution U (t) [19, 20, 21, 22, 23, 24], they largely relyon discretizing the time. Although the truncated Dyson-series algorithm [19, 20], which
is versatile among them, achieves the query complexity optimal or nearly-optimal both int and 1/ε, its scaling is multiplicative with requiring much larger cost than qubitization.
This originates from the difficulty of efficiently dealing with continuous-time modulation.
It is nontrivial whether we can simulate time-dependent Hamiltonians with much fewer
cost close to that of the qubitization.
Among time-dependent Hamiltonians, one of the most important targets of quantum
simulation is a time-periodic Hamiltonian satisfying H(t + T ) = H(t) with T being aperiod. In fact, quantum systems under time-periodic Hamiltonians are called Floquet
systems, and have been platforms of various nonequilibrium phenomena in condensed
matter physics and quantum chemistry; for instance, they can host nonequilibrium phasesof matter absent in time-independent systems, such as topological phases [25, 26, 27]
and time crystals [28, 29, 30, 31, 32, 33]. Time-periodic Hamiltonians also cover opticalresponses in solids and molecules [34, 35, 36, 37], exemplified by high-harmonic generation and photo-chemical reactions. In addition, adiabatic quantum dynamics such as
Thouless pumping [38, 39, 40] and adiabatic state preparation for quantum computation[41, 6, 42] can be regarded as a part of time-periodic Hamiltonian dynamics under sufficiently large period T . Despite fundamental significance and potential application oftime-periodic Hamiltonians in broad fields, there is no Hamiltonian simulation protocol
which can efficiently handle their time-dependency, while their simulation itself is possibleby the truncated Dyson-series algorithm.
In this paper, we develop an optimal and/or nearly-optimal quantum algorithm fortime-periodic Hamiltonian simulation. The key idea of our formalism is mapping timedependent systems to effective time-independent systems in the infinite-dimensional FloquetHilbert space, obtained by the Fourier series expansion in time. Although the infinitedimensionality prohibits its simulation on quantum computers, we resolve this by formulating the Lieb-Robinson bound [43], amplitude amplification, and the qubitization
technique. Consequently, we can simulate the target time-evolved state with arbitrarysmall error and failure probability with a finite-dimensional Hilbert space by preparing
ancillary states labeling Fourier indices.
The resulting query complexity for time-periodic Hamiltonian simulation has optimal
and nearly-optimal scaling respectively in time t and inverse error 1/ε, and significantly,
these contributions appear in an additive way. This implies that the time-periodic Hamiltonian simulation can be executed with the cost sufficiently close to qubitization, which
provides the theoretically best scaling, despite the existence of time-dependency. We also
note that the auxiliary states for Fourier indices accurately reproduce the exact dynamics
with fewer degrees of freedom than those for discretized time, which is employed for generic

Accepted in Quantum 2023-03-16, click title to verify. Published under CC-BY 4.0.

2

time-dependent systems. This leads to smaller query complexity and simpler oracles of
our algorithm compared to the truncated Dyson-series algorithm. In addition, we provide
simulation of the Fermi-Hubbard model under light and adiabatic state preparation as
potential applications. Thus, our protocol will shed light on the promising usage of quantum computers for nonequilibrium quantum many-body phenomena in condensed matterphysics and quantum chemistry, and the optimal control in quantum computation.
The rest of this paper is organized as follows. In Section 2, we briefly review Floquet
theory and the qubitization for this paper to be self-contained. We provide our main
results in Sections 3-9, with firstly summarizing them in Section 3. In Section 4, we derive
the Lieb-Robinson bound in time-periodic systems, and the way to accurately reproduce
the exact dynamics with the Floquet-Hilbert space. Sections 5 and 6 are devoted toderiving the subroutines of the algorithm. They respectively provide the amplification
protocol that can realize the time-evolved state with sufficiently high success probability(Section 5) and efficient implementation of the time-evolution operator in the FloquetHilbert space (Section 6). Section 7 completes the optimal time-periodic Hamiltonian
simulation, and compares its resource with other algorithms for time-independent / dependent Hamiltonian simulation. We provide some promising applications in Section 8,
and concludes our paper in Section 9.

2 Preliminaries
2.1 Review on Floquet theory
We briefly review Floquet theory to analyze Schrödinger equation under a time-periodic
Hamiltonian,di |ψ(t)i = H(t) |ψ(t)i , H(t + T ) = H(t),(1)
dtand its solution|ψ(t)i = U (t) |ψ(0)i ,U (t) = T exp −i

Z t
0



dt0 H(t0 ) .

(2)Here, the time-evolved state of the system of interest is represented by |ψ(t)i, and weassume that it is defined on a finite-dimensional Hilbert space H. The set of statesdim(H){|ψj i}j=1
denotes a certain choice of the basis of H. Floquet theorem says that thesolution |ψ(t)i is always written in the form ofdim(H)|ψ(t)i =X

cα e−iεα t |φα (t)i ,|φα (t + T )i = |φα (t)i ,

(3)

α=1

similar to Bloch theorem for spatially-periodic systems. Here, εα ∈ [−π/T, π/T ) and|φα (t)i ∈ H are respectively called quasienergy and Floquet state. Time-periodicity ofH(t) and |φα (t)i allows the Fourier series expansion as
H(t) =Xm∈Z

Hm e−imωt ,|φα (t)i =

X

−imωt
|φm
,
α ie

(4)

m∈Z

† . We aswith the frequency ω = 2π/T . The hermiticity of H(t) implies H−m = Hm
P
2sume kH(t)k < ∞ (k·k ; operator norm), which results in m∈Z kHm k < ∞. Introducing an auxiliary degree of freedom {|li}l∈Z (sometimes called photon number) to relate

Accepted in Quantum 2023-03-16, click title to verify. Published under CC-BY 4.0.

3Figure 1: Intuitive understanding of the effective Hamiltonian Heff [Eq. (6)] in the Floquet-Hilbert
space. It can be seen as a single-particle quantum system with the potential H0 − lω and the hopping
H−m on a infinite one-dimensional chain. We should introduce the cutoffs |−lmax + 1i and |lmax i tolmax
embed it on quantum computers, which leads to Heff
[Eq. (18)].

|φlα i ↔ |li |φlα i, every pair of εα , |φα (t)i can be obtained by the time-independent eigenvalue equation,Heff 

X|li |φlα i ,(5)

|li hl + m| ⊗ H−m .

(6)|li |φlα i = εα 

X
l∈Z

l∈Z

with the effective Hamiltonian defined by
Heff =

X|li hl| ⊗ (H0 − lω) +

l∈Z

X
l,m∈Z

Intuitively, this time-independent Hamiltonian Heff describes a static one-dimensional system, where each l-th site has potential energy H0 − lω, as Fig. 1. Off-diagonal terms H−mrepresent hopping by −m sites, which can be recognized as either emission or absorption of |m| photons. We remark that the difficulty of time-dependence does not vanishby this mapping; it is translated into infinite dimensionality of the space spanned by{|li |ψj i | l ∈ Z, j = 1, 2, . . . , dim(H)}, which is called Floquet-Hilbert space or Sambespace [44]. While the effective Hamiltonian Heff is often used for identifying εα , |φα (t)i,
it can be employed also for directly determining the dynamics as follows [45];|ψ(t)i =X

e−ilωt hl|e−iHeff t |0i |ψ(0)i .

(7)

l∈Z

We can extract the dynamics under a time-dependent Hamiltonian H(t) from that under
the time-independent one Heff by preparing auxiliary systems labeled by |li. However,
we note that applying optimal algorithms for time-independent systems to simulating
this extended dynamics is not straightforward, since the Floquet-Hilbert space is infinitedimensional.

2.2 Review on qubitization technique
Here, we briefly review so-called qubitization technique to efficiently implement e−iHt for
time-independent Hamiltonian H on the Hilbert space H [18]. It begins with assuming theexistence of block-encoding, that is, we suppose an na -qubit auxiliary state |Gia (called annaoracle state) and a unitary gate O on the space C2 ⊗ H (called an oracle gate) satisfying
hG|O|Gia =H,αα ≥ kHk.

Accepted in Quantum 2023-03-16, click title to verify. Published under CC-BY 4.0.

(8)4The inequality α ≥ kHk comes from the norm of the unitary gate as kOk = 1, and theparameter α represents a typical energy scale of the whole system. Here, the oracles O and|Gia are supposed to be efficiently implemented; the unitary gates O and G, which realizes
the oracle state from a reference state |0ia as G |0ia = |Gia , require at most C ∈ poly (N )elementary gates.
The explicit construction of the oracles has been explored for certain classes of static
Hamiltonians. For instance, when the Hamiltonian H is given by a linear combination ofunitary (LCU) as
H=

jX
max

αj Uj ,αj ≥ 0,(9)

Uj ; unitary,j=1

a possible choice of the oracles is
O=

jX
max

|ji hj|a ⊗ Uj ,|Gia =

jX
max r
j=1

j=1

αj
|jia ,
α(10)maxwith α = jj=1
αj . The number of ancillary qubits na scales as log jmax . Since jmax ispoly (N ) for typical N -site systems, the number of elementary gates for the oracles amounts
to C ∈ poly (N ). LCUs cover various types of Hamiltonians exemplified by quantum spin
systems composed of local Pauli operators and fermionic systems in condensed matterphysics and quantum chemistry [46]. The block-encoding of other Hamiltonians, such assparse-access matrices and purified density matrices, has been also revealed [18, 47, 48].
Once we find out the oracles O and |Gia , we can organize a unitary gate Wq , implemented by O(1) additional ancillary qubits and O(q)-times usage of O, G, and O(na )elementary gates, which satisfies

P

Wq |0i⊗{na +O(1)} |ψi = |0i⊗{na +O(1)} e−iHt |ψi + O(εq ),(11)

for every state |ψi ∈ H. Here, the error εq comes from approximating e−iHt by a certain
degree-q polynomial of H, and decays as εq ∈ O((αt/q)q ). The scaling of q required formaking the error εq smaller than an acceptable error ε is given by the following formula;log(1/ε)q ∈ Θ αt +.
log(e + (αt)−1 log(1/ε))




(12)

In order to determine the explicit form of Wq , we resort to quantum signal processing [17].
It dictates that only poly (q)-time classical computation is required for this purpose.
To summarize the qubitization technique, we require the following computational resources to implement the time-evolution operator e−iHt with an acceptable error ε (Seealso Table 1);• Number of ancillary qubits; na + O(1)• Query complexity; q [Eq. (12)]• Number of elementary gates;


O

αt +

log(1/ε)(C + na ) .
log(e + (αt)−1 log(1/ε))




(13)

Here, the query complexity is defined by the complexity counted by the number of oracles,
corresponding to q in this case, and gives the scaling of the number of elementary gates
below. Significantly, the query complexity takes an additive form as
αt + o(log(1/ε)) ∈ poly (N ) t + o(log(1/ε)).
Accepted in Quantum 2023-03-16, click title to verify. Published under CC-BY 4.0.

(14)
5

The o(log(1/ε)) term scales as log(1/ε)/ log log(1/ε) under 1/ε → ∞ with fixed time t. Its
scaling is known to be optimal both in t and 1/ε for simulating generic time-independent
Hamiltonians [47].
In our algorithm for time-periodic Hamiltonians, we exploit qubitization as a subroutine to implement time evolution under an effective Hamiltonian in Floquet theory. We
can achieve the query complexity for time-periodic systems in the additive form like Eq.
(14), while the conventional truncated Dyson-series algorithm for generic time-dependent
systems needs the query complexity in a multiplicative form αt × o(log(αt/ε)).

3 Summary of this paper
In this section, we overview our main results on efficient quantum simulation of timeperiodic Hamiltonians. We will provide their detailed derivation in the following Sections
4-7.

3.1 Overview of Algorithm
We first show an overview of the algorithm for time-evolution under a time-periodic Hamiltonian H(t + T ) = H(t). Throughout the main text, we assume that the Fourier components of the Hamiltonian H(t) vanishes at a certain cutoff mmax ∈ O(1);Hm = 0,(15)

if |m| > mmax .

(This discussion can be extended to cases where Hm exponentially decays in |m|. See
Appendix C.).
The central attempt in the algorithm is to efficiently simulate Eq. (7),|ψ(t)i =∞X

e−ilωt hl| e−iHeff t |0i |ψ(0)i ,

(16)

l=−∞

on quantum circuits. The time-dependency is erased by mapping the dynamics to the oneon the Floquet-Hilbert space, which results in the following two benefits. First, we do notdiscretize the time with infinitesimal intervals. We instead rely on Fourier indices l ∈ Z,
which are originally discrete, and actually they lead to much better accuracy with thesame auxiliary degrees of freedom. Second, we can use various Hamiltonian simulation
algorithms for time-independent systems. Since the qubitization technique has achieved
the best query complexity in the time t and the inverse error 1/ε, exploiting it for Heff
accelerate the simulation of time-periodic systems.
On the other hand, in order to simulate the time-evolved state |ψ(t)i via Eq. (16) on
quantum circuits, we also have several problems to be tackled. First, the Floquet-Hilbert
space is infinite dimensional. We have to introduce truncation at a certain Fourier index
lmax , and consider a finite set of indices defined by
Dlmax = {−lmax + 1, −lmax + 2, . . . , lmax } ⊂ Z.

(17)

Then, we focus on the approximate dynamics by
lmax
Heff

=

X

|li hl| ⊗ (H0 − lω) +l∈Dlmax

XX|li hl + m| ⊗ H−m ,l∈Dlmax m6=0;l+m∈Dlmax

(18)
|ψ

lmax

(t)i =

X

e

−ilωt

hl| e

lmax
−iHefft|0i |ψ(0)i .

(19)l∈Dlmax

Accepted in Quantum 2023-03-16, click title to verify. Published under CC-BY 4.0.

6

The finite-dimensional space, spanned by 2lmax × dim(H) states {|li |ψj i | l ∈ Dlmax , j =1, . . . , dim(H)}, is called the truncated Floquet-Hilbert space. While lmax → ∞ reproduces|ψ lmax (t)i → |ψ(t)i, it is nontrivial how we should choose lmax to achieve the accuracy|ψ(t)i − |ψ lmax (t)i ≤ ε. The second problem is the small success probability of postselecting the ancillary state. Even if we succeed in approximation with finite lmax , Eq.
P(19) requires a projection to a unnormalized ancillary state l∈Dlmax |li for the state in
the truncated Floquet-Hilbert space,|Ψlmax (t)i = e−iωt

P
l

lmax
|lihl| −iHeffte|0i |ψ(0)i .

(20)

In the actual computation, we need post-selection to an ancillary state given byX
1|almax i = √|li .
2lmax l∈Dlmax(21)
2

Success probability of the projection is given by halmax |Ψlmax (t)i ' hψ(t)|ψ(t)i /(2lmax ) =1/(2lmax ). As lmax increases to ensure the accuracy, the expected time to successfully obtain |ψ(t)i becomes longer in proportion to it. The final problem is about implementationlmax
of e−iHeff t . While the effective Hamiltonian Heff is time-independent, its structure iscomplicated due to the additional degrees of freedom labeled by {|li}. It is nontrivial
whether the optimal Hamiltonian simulation algorithm, i.e. the qubitization technique,provides the optimality for time-periodic systems.
Our algorithm relying on the truncated Floquet-Hilbert space efficiently simulates thetime-evolved state |ψ(t)i with resolving the above problems. The significant developments
are composed of the following steps;(a) Decision of the Fourier index for truncation, lmax , to achieve an allowable error ε.
(b) Amplification by symmetry and oblivious amplitude amplification to enhance the success probability up to 1 − O(ε).




lmax
(c) Efficient implementation of exp −iHefft by the qubitization technique.

In Step (a), we explicitly derive the upper bound of |ψ(t)i − |ψ lmax (t)i in a similar way
to the Lieb-Robinson bound. We show that the choice of lmax , based onlog(1/ε)lmax ∈ Θ γt +,log(e + (γt)−1 log(1/ε))γ = sup(kH(t) − H0 k),




t

(22)
(23)

is sufficient to make the error smaller than ε. Here, γ gives the scale of the time-dependent
terms in H(t). The second step (b) plays a role in amplifying the success probability of thepost selection from 1/(2lmax ) to 1 − O(ε). The first protocol, which we call amplification
by symmetry, exploits the symmetry of Heff which is always present and inherent in timeperiodic systems. It brings the success probability from 1/(2lmax ) to 1/4 − O(ε) onlywith small cost of O(log lmax ) elementary gates. Following this, we apply the oblivious
amplitude amplification [16], reminiscent of Grover’s search algorithm [49]. The successprobability is further amplified from 1/4 − O(ε) to 1 − O(ε). Exploiting these two kinds ofamplification, we can obtain the target state |ψ(t)i only with a little additional resourcethat does not change the scaling. The final step (c) constructs conversion of the effective
Accepted in Quantum 2023-03-16, click title to verify. Published under CC-BY 4.0.

7





lmax
t can be efficiently implemented. The original effective
Hamiltonian so that exp −iHeff
lmax
is not suitable for the qubitization technique, since the additional
Hamiltonian Heff
degrees of freedom |li makes its block-encoding inefficient. To resolve this problem,
we


lmax
derive an alternative effective Hamiltonian such that the desired evolution exp −iHeff t
is accurately reproduced. The block-encoding for the alternative effective Hamiltonian
requires O(1)-times queries to the oracles for each Fourier component {Hm }, which isimportant for achieving optimal or nearly-optimal dependence in t and 1/ε.

3.2 Main results
Here, we summarize the computational resources for computing the dynamics |ψ(t)i under a time-periodic Hamiltonian H(t). We construct two different quantum algorithms
depending on the time scale of the dynamics. The first case is an adiabatic-like case, where
we are interested in the dynamics over O(1) periods as ωt ∈ O(1). We call it “adiabatic”
since adiabatic dynamics under sufficiently large period T is a typical target, while T is
not required to be large. The other case is generic long-time regime, in which we focus
on dynamics over multiple periods as ωt ∈ Ω(1). The algorithms for both cases follow
Steps (a)-(c) of Section 3.1; the number of ancillary qubits is determined by the truncation
order lmax of Step(a), and the number of elementary gates comes mainly from the cost of
lmax
implementing exp −iHeff t in Step (c).
We assume that the Fourier component Hm becomes zero for |m| > mmax with anO(1) constant mmax . The block-encoding for each Hm is supposed to be given by anoracle unitary gate Om and an na -qubit oracle state |Gm i as
hGm |Om |Gm i =

Hm
,
αm

αm > 0.

(24)

Each oracle state |Gm i is generated by trivial states as |Gm i = Gm |0i. We define theenergy scale of the Hamiltonian H(t) by
α=

X

αm ,

(25)

|m|≤mmax

and suppose that {αm }m can be embedded into an O(1)-qubit quantum system by
Gcoef |0i =X|m|≤mmaxr

αm|mi ∈ C2mmax +1 .
α(26)The query complexity is defined by the number of queries to the oracles Om , Gm , and Gcoef ,where each of them is supposed to require at-most C elementary gates. We summarize thecomputational resources for simulating |ψ(t)i in the adiabatic-like regime and the genericlong-time regime respectively by the following theorems.
Theorem 1. (Resource for adiabatic-like regime)
Suppose that we are interested in the time-evolved state |ψ(t)i from an arbitrary
initial state |ψ(0)i over O(1) periods. The computational resources to obtain it with theacceptable error and the failure probability smaller than O(ε) are summarized as follows;• Number of ancillary qubits;
na + O(log(γt) + log log(1/ε)).
Accepted in Quantum 2023-03-16, click title to verify. Published under CC-BY 4.0.

(27)
8• Scaling of query complexity;
αt +log(1/ε).
log(e + {αt + o(log(1/ε))}−1 log(1/ε))(28)

The o(log(1/ε)) term in the above denominator scales aslog(1/ε).
log(e + (γt)−1 log(1/ε))(29)

O(na + log(γt) + log log(1/ε)).

(30)• Additional gates per query;

The query complexity has optimal scaling in the time t and nearly-optimal scaling in the
inverse error 1/ε.
Theorem 2. (Resource for long-time regime)
Suppose that we are interested in the time-evolved state |ψ(t)i from arbitrary initialstates |ψ(0)i over multiple periods ωt ∈ Ω(1). The computational resources to obtain itwith the acceptable error and the failure probability smaller than O(ε) are summarized
as follows;• Number of ancillary qubits;
na + O(log(γ/ω) + log log(ωt/ε)).

(31)• Scaling of query complexity;
αt +ωt log(ωt/ε).
log(e + {α/ω + o(log(ωt/ε))}−1 log(ωt/ε))(32)

The o(log(ωt/ε)) term in the above denominator scales aslog(ωt/ε).
log(e + (γ/ω)−1 log(ωt/ε))(33)

O(na + log(γ/ω) + log log(ωt/ε))(34)• Additional gates per query;

The scaling of the query complexity is optimal in time t for practical problems up topoly (N )-time, while it is formally nearly-optimal. It is nearly-optimal in the inverse error
1/ε.
The parameters α and γ, defined by Eqs. (25) and (23), respectively represent energy
scales of the overall terms and the time-dependent terms in H(t). They typically scale
as α, γ ∈ poly(N ) with the system size N . In contrast, the frequency ω = 2π/T istypically an O N 0 value much smaller than α and γ (The high-frequency cases where ω
are comparable to or larger than α, γ are trivial, as discussed in Appendix D). While theabove theorems are about time-periodic Hamiltonians with vanishing Fourier components
at mmax ∈ O(1) as Hm = 0 (|m| > mmax ), we obtain similar results for Hamiltonians
with exponentially-decaying Fourier components kHm k . e−O(|m|) . See Appendix C for
its derivation.
Accepted in Quantum 2023-03-16, click title to verify. Published under CC-BY 4.0.

9

4 Building an appropriate truncated Floquet-Hilbert space
This section is devoted to deriving the proper truncation order for the Floquet-Hilbert
space, following (a) in Section 3.1. In order to achieve desirable accuracy for the exact
time-evolved state with the truncated Floquet-Hilbert space, we should obtain an exact
upper bound on its error. Here, we prove the Lieb-Robinson bound in the Floquet-Hilbert
space in Section 4.1, and determine a proper truncation order based on it in Section 4.2.

4.1 Lieb-Robinson bound in Floquet-Hilbert space
The proper choice of the truncation order lmax is determined so that the approximate state|ψ lmax (t)i [See Eq. (19)] can reproduce the time-evolved state |ψ(t)i with satisfying therelation |ψ(t)i − |ψ lmax (t)i ≤ ε. From the explicit formula Eq. (19), we can see that it islmax

important to observe how the transition amplitude hl|e−iHeff t |0i behaves for sufficiently
large l. We first show the upper bound on this transition amplitude in Theorem 3. Since
it is reminiscent of the Lieb-Robinson bound in single-particle quantum systems [50], wecall it the Lieb-Robinson bound in the Floquet-Hilbert space.
Theorem 3. (Bound on transition amplitude)
We assume Hm = 0 for |m| > mmax . Then, for indices l, l0 ∈ Z satisfying |l − l0 | ≥2mmax γt, the transition amplitude is bounded by0lmax

hl| e−iHeff

t

|l0 i

≤ 2

(γt)d|l−l |/mmax e(d|l − l0 |/mmax e)!



≤

γt
0
d|l − l |/mmax e(35)

d|l−l0 |/mmax e

,

(36)

where the parameter γ is defined by Eq. (23).
Proof.— We will omit the superscripts lmax for some operators introduced here since
they are not important. We first employ the interaction picture. With the unitary operation defined byX
U0 (t) = e−iH0 t , H0 =|li hl| ⊗ (H0 − lω),(37)
l∈Dlmax
lmax

the time evolution operator e−iHeff
lmax

e−iHefftt

is represented as= U0 (t)UI (t),UI (t) = T exp −i

(38)

Z t
0



dt0 HI (t0 ) .

(39)Here, the Hamiltonian in the interaction picture is defined bylmax
HI (t) = U0 (t)† (Heff− H0 )U0 (t)

=
I
H−m
(t) =

XX

I
(eimωt |li hl + m| ⊗ H−m(t) + h.c.),

l m>0
iH0 t
e
H−m e−iH0 t ,

(40)
(41)

where the summation in the second line is taken over l, m such that l, l +m ∈ Dlmax . Using

Accepted in Quantum 2023-03-16, click title to verify. Published under CC-BY 4.0.

10

the Dyson series expansion, we obtainlmax

hl| e−iHeff

t

|l0 i

hl|UI (t)|l0 i

=

∞ Z t
X

≤

dtn . . .

n=0 0
∞ Z t
X

=

dtn . . .

n=0 0

Z t2
0

Z t2
0

dt1 hl|HI (tn ) . . . HI (t1 )|l0 i
n
XY

dt1

hli |HI (ti )|li−1 i .

(42){li } i=1

In the last equality, we employ the identity li ∈Dlmax |li i hli | = I for n − 1 times. The
P
summation {li } is taken over li ∈ Dlmax for i = 1, 2, . . . , n − 1 under the fixed l0 = l0
Q
and ln = l. Each product ni=1 hli |HI (ti )|li−1 i represents a complex transition amplitude
from |l0 i to |li along the path |l0 i → |l1 i → . . . → |ln−1 i → |li under HI (ti ). Since the
Hamiltonian HI (ti ) shifts the Fourier index |li i by at-most mmax due to Eq. (18), the loworder terms labeled by n < |l − l0 |/mmax disappears. This results in
P

hl| e

lmax
−iHeff
t

0

|l i

∞
X

≤

Z t

dtn . . .

0n=d|l−l0 |/mmax e

∞
X
tn

≤

n=d|l−l0 |/mmax e

n!

Z t2
0

dt1 hl|nY

HI (ti )|0i

i=1

n

sup(kHI (t)k)

(43)

t

0

I (t) −
Here, HI (t) is a Toeplitz matrix that satisfies hl|HI (t)|l0 i = e−i(l−l )ωt Hl−l0
H0 δl−l0 ,0 . This leads to the upper bound of its operator as follows [51],kHI (t)k ≤ sup 
t0

X

0
I
(eimωt H−m(t) − H0 δm,0 )eimωt m∈ZH(t − t0 ) − H0

= sup
t0

= γ.

(44)

if n0 ≥ 2x ≥ 0.

(45)



We also use the following inequality;∞Xxn
n=n0 n!

≤2

xn0
,
n0 !

This relation can be easily confirmed by∞X
xn
n=n0 n!

=

∞
∞
xn0 X
xn−n0
xn0 X
1
≤
n0 ! n=n0 n(n − 1) . . . (n0 + 1)
n0 ! n=0 2

 n

.

(46)

For the indices l, l0 satisfying |l − l0 | ≥ 2mmax γt, we can substitute x = γt and n0 =d|l − l0 |/mmax e. This results in the bound,0

lmax

hl| e−iHefft|l0 i ≤ 2(γt)d|l−l |/mmax e.
(d|l − l0 |/mmax e)!

(47)

The inequality from the Stirling formula,n0
n0 ! ≥ 2e

ensures the inequality Eq. (36).

n0

,

(48)



Accepted in Quantum 2023-03-16, click title to verify. Published under CC-BY 4.0.

11

The Lieb-Robinson bound in the transition amplitude provides a guide for choosing aproper truncation order. By setting l0 = 0, Eqs. (35) and (36) say that the contributions to|ψ lmax (t)i from indices satisfying |l| ≥ 2mmax γt rapidly decay as l−l . Thus, O(γt) becomesa possible truncation order. However, we note that the change in lmax affects |ψ lmax (t)ilmax
. By taking it into
also via the change in the support of the effective Hamiltonian Heff
account with using the Lieb-Robinson bound, we obtain the exact upper bound on theerror between |ψ(t)i and |ψ lmax (t)i as follows.
Theorem 4. (Floquet-Hilbert space truncation)
We consider the approximate time-evolved state obtained from the truncated FloquetHilbert space, |ψ lmax (t)i, given by Eq. (19). Then, its deviation from the exact one |ψi isbounded by(γt)dlmax /mmax e|ψ(t)i − |ψ lmax (t)i ≤ 20mmax,(49)(dlmax /mmax e)!
if the truncation order lmax satisfies lmax ≥ 2mmax γt.
0
Proof.— We evaluate the convergence of |ψ lmax (t)i. For different orders lmax , lmax
0satisfying lmax
> lmax , we compute the difference,0|ψ lmax (t)i − |ψ lmax (t)i ≤ ε1 + ε2 ,

ε1 =

l0max

X

hl| eiHeff

t

(50)

(51)

|0i ,

0
l∈(Dlmax \Dlmax )

ε2 =

l0max
t

X

hl| eiHeff

lmax− eiHeff

t

(52)

|0i .

l∈Dlmax

The first error comes from truncating the order of the post-selected state |almax i. Using
Theorem 3 directly implies
X

ε1 ≤

0
l∈(Dlmax \Dlmax )∞X(γt)d|l|/mmax e(γt)dl/mmax e≤2.
(d|l|/mmax e)!
(dl/m
e)!
max
l=l

(53)

max

The summation ∞
l=lmax can be divided based on l (mod. mmax ), and each summation
corresponds to the left hand side of Eq. (45) with n0 ≥ dlmax /mmax e. In other words, forlmax ≥ 2mmax γt, we obtain the following inequality,P∞X(γt)dl/mmax el=lmax(dl/mmax e)!

≤ 2mmax(γt)dlmax /mmax e.
(dlmax /mmax e)!

(54)

As a result, the first error ε1 is bounded by 4mmax (γt)dlmax /mmax e /(dlmax /mmax e)!.
The second error ε2 comes from truncating the order of the effective Hamiltonian. In
a similar way to the proof of Theorem 3, each term is bounded byl0max
thl| eiHeff

lmax− eiHeff

t

|0i ≤

∞ Z t
X
n=0 0

dtn . . .

Z t2
0dt1 hl|nY

l0

HI max (ti ) −

i=1

n
Y

HIlmax (ti )|0i ,

i=1

(55)
0
where we explicitly write the size lmax , lmax
for HI (t) [See Eq. (40)]. In a similar mannerQQl0
to Eq. (42), hl| ni=1 HI max (ti )|0i and hl| ni=1 HIlmax (ti )|0i respectively represent transition amplitudes summed over all possible paths |0i → |l1 i → . . . → |ln−1 i → |li underAccepted in Quantum 2023-03-16, click title to verify. Published under CC-BY 4.0.

12

l0

lmax
(t). Their differences survive only when the path goes across the doHeffmax (t) and Heff
0
l
l
max
main D max \D, where these two Hamiltonians have different actions. In other words,
low order terms with n < {(lmax − 0) + (lmax − |l|)/mmax should vanish. Considering thatl0lmax
the norm of Heff(t) and Heffmax (t) is both bounded by γ, we obtain
∞
X

X

ε2 ≤

l∈Dlmax n=d(2lmax −|l|)/mmax

≤ 4

tn n
2γ
n!
e

(γt)d(2lmax −|l|)/mmax e(d(2lmax − |l|)/mmax e)!
lmax

X
l∈D

≤ 16mmax(γt)dlmax /mmax e.
(dlmax /mmax e)!

(56)

0
In the last inequality, we again use the relation Eq. (54). Taking the limit lmax
→ ∞ forε1 + ε2 reproduces the result of Theorem 4. 

4.2 Truncation order of Floquet-Hilbert space
We hereby determine the truncation order lmax so that |ψ lmax (t)i can reproduce the exact
time-evolved state with a desirable error up to O(ε). By using the inequality from theStirling formula, Eq. (48) to Theorem 4, we obtain the error bounded by|ψ(t)i − |ψ lmax (t)i ≤ 10mmax

emmax γtlmax

lmax /mmax

(57)

,

for lmax ≥ 2mmax γt. The truncation order lmax is chosen so that the right hand side canbe smaller than ε. To this aim, we should evaluate a function f (x) = (κ/x)x [κ > 0],
which is known to be dealt with the Lambert W function W (x) satisfying W (x)eW (x) = x[52]. Here, we rely on the resulting proposition [47].
Proposition 5.
The function f (x) = (κ/x)x [κ > 0] is monotonically decreasing in x ≥ κ/e, andsatisfies the following inequality for 0 < η < 1;
f (x) ≤ η,∀x ≥ eκ +4 log(1/η).
log(e + κ−1 log(1/η))(58)

See Lemma 59 in the full version of Ref. [47] for the proof. Based on the above
proposition, we choose lmax by
4mmax log(10mmax /ε)lmax = mmax + e mmax γt +,log(e + (eγt)−1 log(10mmax /ε))




2

(59)

so that the error can be bounded from above as
10mmax

emmax γtlmax − mmax

(lmax −mmax )/mmax

≤ ε.

(60)

We note that this choice does not violate the condition lmax ≥ 2mmax γt, which is required
for Theorem 4. The monotonicity of f (x) ensures |ψ(t)i − |ψ lmax (t)i ≤ ε under theabove choice (The additional term mmax is attached for later calculation, especially for
Appendix B.1) and B.2).
Accepted in Quantum 2023-03-16, click title to verify. Published under CC-BY 4.0.

13

Let us discuss how lmax increases in the time t and the inverse error 1/ε. The formof the error in lmax , given by Eq. (49), is the same as that for qubitization in the querycomplexity q, given by εq ∈ O((αt/q)q ). Therefore, its scaling can be determined in asimilar way, which results inlog(1/ε).
lmax ∈ Θ γt +log(e + (γt)−1 log(1/ε))




(61)

To reproduce the truncated Floquet-Hilbert space, we should prepare an ancillary systemlabeled by {|li}l∈Dlmax . The number of qubits for such ancillary system is at-mostO(log lmax ) ⊂ O(log(γt) + log log(1/ε)).

(62)

5 Amplitude amplification of the dynamics
In this section, we show two kinds of amplification to achieve sufficiently high successprobability of extracting the time-evolved state |ψ(t)i.
As we stated in Section 3.1, the approximate state |ψ lmax (t)i cannot be directly realizedon quantum circuits. After simulating the dynamics in the truncated Floquet-Hilbert
space to get |Ψlmax (t)i defined by Eq. (20), we make a projection to |almax i. Although
the resulting renormalized state is sufficiently close to |ψ lmax (t)i and also the time-evolved
state |ψ(t)i as
halmax |Ψlmax (t)i
|ψ lmax (t)i=
= |ψ(t)i + O(ε)(63)
khalmax |Ψlmax (t)ikk|ψ lmax (t)ik
for lmax given by Eq. (59), we should be careful of the low success probability; it becomes2

halmax |Ψlmax (t)i ∈ O (lmax )−1 .
We need approximately O(γt)-times trials of the post selection, and every trial islmax
expected to require at-least O(t) complexity for implementing |Ψlmax (t)i = e−iHeff t |0i ⊗|ψ(0)i. Therefore, the naive implementation based on Eq. (19) is not efficientfor the
time-evolved state |ψ(t)i in that the expected computational time reaches O t2 .
We resolve this problem by the amplification of the success probability up to 1−O(ε) inthe next section. The first one,which exploits the symmetry of the effective Hamiltonian−1
to O(1). The latter one following this, which is reminiscent
Heff , amplifies from O lmax
of the Grover’s search algorithm, allows the success probability 1 − O(ε). As discussed
later, using only either one fails to efficiently compute |ψ(t)i.


5.1 Amplification by symmetry
We introduce the amplification exploiting the symmetry of the effective Hamiltonian Heff .
First of all, we specify the symmetry here; it is about the translation of the photon number
from |li to |l + mi. When we define the translation operator on the Floquet-Hilbert spacePby Tm = l∈Z |li hl + m| ⊗ I, Heff satisfies the translation symmetry,Tm† Heff Tm = Heff + mω,hl|e−iHeff t 0|l i = eil0 ωt0(64)−iHeff t

hl − l |e

|0i .

(65)

Since the symmetry is present as long as the time-periodicity H(t + T ) = H(t) holds, theamplification discussed here is always available in our algorithm.
Here, we provide two tasks for the amplification by symmetry. First, we slightly extend
the truncated Floquet-Hilbert space to C8lmax ⊗ H, where the Fourier index |li is chosen
Accepted in Quantum 2023-03-16, click title to verify. Published under CC-BY 4.0.

14

from l ∈ D4lmax = {−4lmax + 1, −4lmax + 2, . . . , 4lmax }. The number of additional qubits
required for this extension from Clmax ⊗ H is three, and hence this does not affect thescaling of the computational resources. The second task is to preprocess the initial state
by
lmax
Uini
|0i |ψ(0)i = |almax i |ψ(0)i ,
(66)
lmax
where the unitary operator Uiniis expressed by





lmax
Uini
= |almax i h0| + . . . ⊗ I.

(67)

lmax
The unitary operator Uininontrivially acts only on the ancillary space, and it can belmax
implemented with O(log lmax ) gates (e.g. for lmax such that log2 lmax ∈ N, we have Uini=
⊗ log2 (2lmax )Had⊗ I with the Hadamard gate Had). This process plays a role of making theinitial state approximately translation-invariant with the width of Fourier indices, 2lmax .
The remaining procedures after the above two tasks is the same as those for the original
protocol in Section 3.1. Reflecting that the ancillary Hilbert space is extended to C8lmax ,P
4lmax
and l∈D4lmax lω |li hl| ⊗ I, andwe evolve the above uniform state |almax i |ψ(0)i by Heff
make a projection to |a4lmax i. The resulting state |ψ̄ lmax (t)i is calculated as follows.

|ψ̄ lmax (t)i ≡ ha4lmax | e−it

P
l

4lmax
lω|lihl| −iHeffte|almax i |ψ(0)i .

(68)

The point of this process is that the above state |ψ̄ lmax (t)i can reproduce |ψ(t)i with
O(1) amplitude owing to the translation symmetry of the effective Hamiltonian, Eqs.
(64) and (65). We give its brief explanation in this section, while the rigorous derivation
is provided in Appendix B.1. Although we do not have the exact translation symmetry
represented by Eqs. (65) due to the finite-size effect of lmax , it is expected to approximately
hold as
4lmax
4lmax
0
(69)
hl| e−iHeff t |l0 i ' eil ωt hl − l0 | e−iHeff t |0i .
if the truncation order lmax is sufficiently large. As a matter of fact, we can derive theexact upper bound on the difference between the left- and right-hand sides based on the
Lieb-Robinson bound (See Lemma 8 in Appendix B.1). We proceed the discussion withassuming the approximate relation Eq. (69) here. The resulting state of the process|ψ̄ lmax (t)i can be roughly computed as follows;|ψ̄ lmax (t)i =

1

X

X

4lmax 04lmax

e−ilωt hl| e−iHefft|l0 i |ψ(0)il ∈Dlmax l∈D4lmax

'

1

X

X

4lmax 004lmax

e−i(l−l )ωt hl − l0 | e−iHefft|0i |ψ(0)il ∈Dlmax l∈D4lmax

=

1

X

4lmax 0(|ψ(t)i + O(ε))

l ∈Dlmax=1|ψ(t)i + O(ε).
2(70)

The second equality comes from Theorem 4, considering that the summation of l − l0 overl ∈ D4lmax is sufficient to suppress the error up to O(ε). As a result, we obtain1|ψ̄ lmax (t)i ' (|ψ(t)i + O(ε)).
2(71)

Let us focus on the success probability of the projection onto |a4lmax i. It is provided
by hψ̄ lmax (t)|ψ̄ lmax (t)i = 1/4 + O(ε), which is much larger than the original one 1/(2lmax ).
Accepted in Quantum 2023-03-16, click title to verify. Published under CC-BY 4.0.

15

Therefore, the amplification protocol by symmetry enables us to avoid O(t)-times repetition of the time evolution e−iHeff t in contrast to the original protocol. Intuitively, thisdrastic improvement can be understood as a result of the interference under the translation
symmetry, as Fig. 2. As we can see from Eq. (20), the target state |ψ(t)i is extracted from|Ψlmax (t)i via the accompanying ancillary state |almax i, which is uniform in the Fourier index |li. When we begin with the non-uniform initial state |0i |ψ(0)i,the resulting state√Pafter the time evolution is also non-uniform. Since it involves (1/ 2lmax ) l∈Dlmax eikl |li
for k ∈ (2π/2lmax )Z (different eigenstates of the translation operator) with approximate
lmax i is comparably small as
equal√ weight, the amplitude of the desirable component |a
1/ 2lmax . In contrast, when we employ the uniform initial state |almax i |ψ(0)i, only the
uniform components that includes |almax i are amplified while the other components cancels
one another. It can be viewed also as the interference of different initial states |l0 i |ψ(0)iin |almax i |ψ(0)i as shown in Fig. 2, and this is why the amplification protocol achieves
O(1) amplitude of |almax i |ψ(0)i.
We summarize the amplification by symmetry with adding the exact results obtainedin Appendix B.1. We prepare the truncated Floquet-Hilbert space C8lmax ⊗ H, and makethe initial state uniform as |almax i |ψ(0)i. As a result, this process is described by a unitary
gate
P
4lmax
lmax
4lmax † −it l lω|lihl| −iHeff
lmax
t
Uamp1
(t) = (Uini
)e
e
Uini
,
(72)
4lmax †
where the last unitary gate (Uini) is added to replace the projection to |a4lmax i by
the one to |0i. This approximately realizes the time-evolved state |ψ(t)i with 1/4 + O(ε)
success probability as

1lmax
h0| Uamp1(t) |0i |ψ(0)i ' (|ψ(t)i + O(ε)).
2(73)

The probability 1/4 comes from the ratio of the width of the initial state |almax i to thatof the projected ancillary state |a4lmax i. The exact version, derived based on the LiebRobinson bound in Appendix B.1, is stated as follows, and ensures the validity of the
discussion here.
Theorem 6. (Amplification by symmetry)
We designate the truncation order lmax ∈ Θ(γt + log(1/ε)/ log log(1/ε)) by Eq. (59).
Then,
1
ε
lmax
h0| Uamp1(t) |0i |ψ(0)i − |ψ(t)i ≤ .
(74)
2
3
is satisfied for an arbitrary initial state |ψ(0)i ∈ H.
We finally note that the amplification solely by this process cannot achieve the successprobability 1−O(ε) with keeping the efficiency. Analogous to the above formulation, when
we set the widths of the supports of the initial ancillary state, the effective Hamiltonian,and the projected ancillary state to plmax , qlmax , qlmax (p, q ∈ N, q ≥ p + 2) respectively,
the success probability becomes p/q. However, the probability larger than 1 − ε demands
the relation,
p
p
.
(75)
1−ε< ≤q
p+2
This implies p ≥ 2(1 − ε)/ε, and hence the number of ancillary qubits reaches at-leastlog(plmax ) ∈ O(log(γt/ε)). In addition, the query complexity includes a term proportional to the dimension of the ancillary system, qlmax as discussed later in Section 6. Thisimplies that the query complexity increases linearly in 1/ε, destroying the original logarithmic scaling. Thus, relying solely on this amplification protocol is not suitable for
Accepted in Quantum 2023-03-16, click title to verify. Published under CC-BY 4.0.

16

(a)

(b)

Boundaryterms

Figure 2: Schematic picture of how the dynamics in the truncated Floquet-Hilbert space provides
the time-evolved state |ψ(t)i. (a) The dynamics when we begin with the initial state |0i |ψ(0)i [SeeEq. (20)]. Following the Lieb-Robinson bound, Theorem √3, the state spreads out within the region|l| . lmax . This spread results in the low amplitude O 1/ lmax of the target output |almax i |ψ(t)i.
(b) The dynamics under the amplification by symmetry, where we employ the uniform initial state|almax i ||ψ(0)ii as Eq. (68). Interference takes place within the region |l| . lmax , and thereby enhances
the amplitude of the uniform target state |a4lmax i |ψ(t)i up to O(1). The gray region represents the
4lmax
support of boundary terms in the refined effective Hamiltonian Heff,pbc[See Eq. (101)]. It is sufficiently
far from the light-blue region relevant to the dynamics, which is determined by the Lieb-Robinson bound.
This results in the validity of the refined effective Hamiltonian as Theorem 7 says.

efficient simulation of the dynamics. We resolve this problem by the oblivious amplitude
amplification below.

5.2 Oblivious amplitude amplification
The above amplification
based on the translation symmetry of Heff enhances the amplitude√of |ψ(t)i from 1/ 2lmax to 1/2−O(ε). Here, we introduce another amplification, called the
oblivious amplitude amplification [16], to achieve the amplitude (or the success probability)
1 − O(ε).
The starting point of this amplification is the result of the previous section 5.1. Equation (73), or equivalently Theorem 6, indicates that the consequence of the amplification
by symmetry is written as
lmax
Uamp1
(t) |0i |ψ(0)i =1|0i (|ψ(t)i + O(ε)) + |Ψ⊥ i ,2(76)with an additional term |Ψ⊥ i ∈ C8lmax ⊗ H satisfying(|0i h0| ⊗ I) |Ψ⊥ i = 0.

(77)

The state |Ψ⊥ i generally depends on |ψ(0)i and t. The oblivious amplitude amplification
takes a similar strategy to that of the Grover’s search algorithm, in which we compose thefollowing two unitary operators;
R = (2 |0i h0| − I) ⊗ I,
lmax
Uamp2
(t)

=

lmax
lmax
lmax
−Uamp1(t)R[Uamp1(t)]† RUamp1(t).

Accepted in Quantum 2023-03-16, click title to verify. Published under CC-BY 4.0.

(78)
(79)

17

The first one R reverses the sign of |li for l 6= 0, and it is implemented with O(log lmax )lmax
gates. The second one Uamp2(t) plays a role in enhancing the amplitude of |ψ(t)i up to
1 − O(ε). Its action on any initial state |0i |ψ(0)i is computed as follows;
1
lmax
lmax
lmax †
Uamp2|0i |ψ(0)i = Uamp1
R(Uamp1) − |0i |ψ(t)i + |Ψ⊥ i + O(ε)
2




n

o

lmax
lmax †= Uamp1
R |0i |ψ(0)i − (Uamp1) |0i |ψ(t)i + O(ε).

(80)

lmax †In the second equality, we use the relation obtained by applying (Uamp1) to Eq. (76).
Next, we evaluatelmax †
lmax †
lmax †
R(Uamp1) |0, ψ(t)i = 2 |0i h0| (Uamp1) |0i |ψ(t)i − (Uamp1) |0, ψ(t)i .

(81)

Theorem 6 indicates that the time evolution operator U (t) is approximated asε1lmax
h0|Uamp1(t)|0i − U (t) ≤ ,
2
3

(82)

and hence the relation

1
ε
lmax
h0|[Uamp1(t)]† |0i − U (t)† ≤2
3is also satisfied. This provides the relation,(83)

lmax †
lmax †
R(Uamp1) |0i |ψ(t)i = |0i |ψ(0)i − (Uamp1) |0i |ψ(t)i + O(ε),(84)

and substituting this into Eq. (80) results in
lmax(t) |0i |ψ(0)i = |0i |ψ(t)i + O(ε),
Uamp2

(85)

for an arbitrary initial state |ψ(0)i.
lmax
This result indicates that the operation Uamp2(t) generates the time-evolved state|ψ(t)i with the amplitude 1 − O(ε). Or equivalently, as an exact bound for the error, we
can derive the inequality,
lmax
Uamp2
(t) |0i |ψ(0)i − |0i |ψ(t)i ≤ ε.

(86)

The coefficient comes from the fact that an error bounded by ε/3 appears due to Theoremlmax
6 every time we call Eq. (76). The amplification protocol Uamp2(t) employs 3 times queries
lmax
lmax
to Uamp1 (t). Reflecting that the operations Uini and R require relatively a little resource
lmax(at-most O(log lmax ) elementary gates and complexity), the resource for Uamp2(t) has thePsame scaling as the one for implementing the time evolution operators e−it l lω|lihl| and4lmax
e−iHeff t .
As well as the amplification by symmetry, relying only on the oblivious amplitude
amplification fails to efficiently enhances the success probability to obtain |0i ⊗ |ψ(t)i−1
from O lmax
to 1 − O(ε). When we do not use the first amplification, the protocol of the
oblivious amplitude amplification is given by
lmax
lmax
lmax †U˜amp2,p= − Uorig
R(Uorig
)R

n

op

lmax
Uorig
,

(87)

lmax
where the operation Uorig
represents the time evolution without the amplification bysymmetry, defined bylmax
lmax † −it
Uorig(t) = (Uini
)e

P
l

lmax
lω|lihl| −iHeffte.

Accepted in Quantum 2023-03-16, click title to verify. Published under CC-BY 4.0.

(88)
18

lmax
Applying U˜amp2,pto the initial state |0i ⊗ |ψ(t)i returns |0i ⊗ |ψ(t)i whose amplitude√
increasesfrom 1/ 2lmax approximately in proportion to p. The integer p should be√O lmax to achieve the amplitude 1 − O(ε), P
reminiscent of the Grover’s search algorithm.
√
lmax−it l lω|lihl|
In other words, O lmax times call of e
and e−iHeff t is required when welmax
use U˜amp2,p. This is why we suggest the combination of the two amplification protocols,

P

lmax

−it l lω|lihl|
with which O(1) times
and e−iHeff√ call of e|0i ⊗ |ψ(t)i from 1/ 2lmax to 1 − O(ε).

tcan amplify the amplitude of

6 Block-encoding of Effective Floquet HamiltonianIn the previous section, we show that the combination of the two kinds of amplification
lmax
protocols, implemented by Uamp2(t) [See Eq. (86)], provides the target time-evolved statelmax|ψ(t)i with an arbitrarily small error O(ε). The computational resource for Uamp2(t) isP

4lmax

mostly determined by that for e−it l lω|lihl| and e−iHeff t . Our strategy is to implement
these two time evolution operators by the qubitization technique [18]. As discussed in
Section 5.2, the computational resources are determined by how to introduce the blockP
lmax= l∈Dlmax lω |li hl| ⊗ I (linear potential
encoding of the two static Hamiltonians, HLP
lmax(effective Hamiltonian). The aim of this section is to obtain
Hamiltonian) and Heff
P an ef−it l lω|lihl|
ficient block-encoding of them and to evaluate the costs required to implement e
4lmax
and e−iHeff t .

6.1 Block-encoding of linear potential Hamiltonian
We compose block-encoding of the linear potential Hamiltonian,
4lmax
HLP
=

X

(89)

lω |li hl| ⊗ I.

l∈D4lmax

By simple calculation, it can be written in the form of an LCU,
4lmax
HLP

=

l∈D

Vl4lmax=

ω 4lmax
V
⊗ I,
2 l
4lmax

(90)

X
4lX
max

l−1
X

|l0 i hl0 | −l0 =l|l0 i hl0 | .

(91)

l0 =−4lmax +1

Using the block-encoding formalism of LCUs by Eq. (10), we immediately obtain theoracle unitary gate and the oracle state by
4lmax
OLP
=

X

(|li hl|)b ⊗ Vl4lmax ⊗ I,

(92)l∈D4lmax
4lmax 4lmax
ha4lmax |OLP
|a
ib =

4lmax
HLP
.
4lmax ω(93)

The subscript b for the states |lib and |a4lmax ib represents a new ancillary system introduced
for the block-encoding, requiring the number of qubits nb ∈ O(log lmax ). Combining Eqs.
(91) and (92), the oracle unitary gate is rewritten by

4lmax
OLP
=


X
l,l0 ;l0 −l≥0|l, l0 i hl, l0 | −X|l, l0 i hl, l0 | ⊗ I.

(94)

l,l0 ;l0 −l<0

Accepted in Quantum 2023-03-16, click title to verify. Published under CC-BY 4.0.

19

To implement this oracle, we use a comparator defined by(Comp |l, l0 i |0ib0 =|l, l0 i |0ib0|l, l0 i |1ib0

if l0 ≥ l,
if l0 < l,(95)

with a single-qubit ancillary system b0 [53]. We can immediately confirm the relation,
4lmax|l, l0 i |0ib0 |ψi ,Comp† Zb0 Comp |l, l0 i |0ib0 |ψi = OLP(96)

for arbitrary inputs l, l0 ∈ D4lmax and |ψi ∈ H, where Zb0 denotes a Pauli Z operator on
the system b0 . Since the comparator on n qubits can be composed of O(n) elementary
4lmax
requires at-most O(log lmax ) gates. The oracle state
gates, the oracle unitary gate OLP
4l
max
|a
ib , which has equal weights in |li, can be prepared by O(log lmax ) elementary gates.
4lmax
We implement the time-evolution e−iHLP t with accuracy 1 − O(ε) by qubitization.
The number of additional gates, other than the block-encoding, is proportional to the
4lmax
number of ancillary qubits for the oracle OLP, that is, we need O(log lmax ) additionalgates per query. Based on Section 2.2, it can be executed with the following resources;• Number of ancillary qubits; O(log lmax )• Scaling of overall complexity;log(1/ε)lmax ωt +log lmax .
log(e + (lmax ωt)−1 log(1/ε))(97)

6.2 Block-encoding of effective Hamiltonian
4lmax

We hereby provide a way to efficiently implement e−iHeff t by the qubitization technique.
4lmax, given by Eq.
We find out that naive block-encoding of the effective Hamiltonian Heff(18), faces at the severe increase of the complexity in lmax . This originates from the fact
4lmax
that the number of nontrivial terms in Heff
increases linearly in lmax . To avoid thisproblem, we take an alternative approach which yields only O(log lmax ) complexity as
follows;
1. Find a refined effective Hamiltonian H on the truncated Floquet-Hilbert space
4lmax
C8lmax ⊗ H, such that e−iH t accurately reproduces e−iHeff t with an arbitrarily
small error O(ε).
2. Construct block-encoding of H with composing of efficient implementation of its
oracle unitary gate and oracle state.
First, we specify the assumption for the time-periodic Hamiltonian H(t). Here, we
suppose two kinds of the oracles that can be accessed. The first one is about the blockencoding of each Fourier component Hm (|m| ≤ mmax ), given by
hGm |a Om |Gm ia =

Hm
,
αm

αm > 0,(98)
nawith an na -qubit ancillary system. Each oracle unitary gate on C2 ⊗ H is representednaby Om , and each oracle state |Gm ia ∈ C2 is generated by a quantum circuit Gm as|Gm ia = Gm |0i⊗na . The second oracle is a query to coefficients of Fourier components,

Accepted in Quantum 2023-03-16, click title to verify. Published under CC-BY 4.0.

20

represented by αm . We define the oracle unitary circuit Gcoef , acting on C2mmax +1 [i.e.
O(1)-qubit system], by
Gcoef |0i =X|m|≤mmax

r

αm
|mi ,
α

α=

X

αm .

(99)

|m|≤mmax

The parameter α, defined by α = |m|≤mmax αm , characterizes the energy scale of the
Hamiltonian H(t). In fact, considering kHm k ≤ αm for each m, it is bounded by kH(t)k ≤P|m|≤mmax kHm k ≤ α. The parameter α also provides the bound on the energy scale oftime-dependent terms γ, as
P

X

γ≤

(100)

kHm k ≤ α.

m6=0;|m|≤mmax

Implementation of the oracle Gcoef is equivalent to embedding probability distributions inquantum states [54], which also appears as the oracle state for the qubitization techniquein the case of a LCU [See Eq. (10)]. Stimulated by its various applications covering linear
equation solver [55] and quantum singular value transformation [47, 48], there have beenvarious efficient implementations [56, 54, 46, 57, 58]. Here, we suppose that the number of
elementary gates required for the oracles Om , Gm , and Gcoef is at-most C. As discussed
later in Section 8, the form of time-periodic Hamiltonians H(t) designated by Eq. (98) isreasonable in that they involve various important classes such as LCUs.
6.2.1

Construction of a refined effective Hamiltonian

Let us compose a refined effective Hamiltonian, which accurately reproduces the dynamics
4lmax
of e−iHeff t . To come to the point, such a Hamiltonian equipped with suitability for
efficient block-encoding is given by
4lmax
4lmax
+
= HeffHeff,pbc

X(|li hl ⊕ m| ⊗ H−m + h.c.),(101)(l,m)∈∂F 4lmax

with ∂F 4lmax = {(l, m) | 4lmax − mmax + 1 ≤ l ≤ 4lmax , 4lmax − l + 1 ≤ m ≤ mmax }. The
4lmax
integer l ⊕ m ∈ D4lmax is defined modulo 8lmax . The additional terms in Heff,pbc
connects
4l
max
the boundary of D
so that the hopping terms, |li hl + m|, induced by H−m can betranslation-invariant under the shift of |li. The subscript “pbc” comes from the periodic
boundary conditions (PBC) for the hopping terms.
4lmax

In our algorithm, we employ the time evolution e−iHeff,pbc t instead of the original
4lmax
one e−iHeff t , which is derived by Floquet theory. In other words, we organize the
4lmax
amplification protocols with the refined effective Hamiltonian Heff,pbc
as
4lmax

lmax
4lmax † −iHLP
Uamp1,pbc(t) = (Uini
)e

4l

max t
t −iHeff,pbc

lmax
Uini
,

(102)

lmax
lmax
lmax
lmax
Uamp2,pbc(t) = −Uamp1,pbc
R(Uamp1,pbc)† RUamp1,pbc
,

(103)

e

and apply them to the initial state |0i |ψ(0)i. The refined effective Hamiltonian is valid
in a sense that it can also provide the exact time-evolved state as
lmax
Uamp2,pbc(t) |0i |ψ(0)i = |0i |ψ(t)i + O(ε).

(104)

The rigorous upper bound on the error is provided by the following theorem.

Accepted in Quantum 2023-03-16, click title to verify. Published under CC-BY 4.0.

21

Theorem 7. (Refined effective Hamiltonian)
We choose lmax ∈ O(γt + log(1/ε)/ log log(1/ε)) by Eq. (59). Then, the refined ef4lmax
4lmax
fective Hamiltonian Heff,pbc
well reproduces the dynamics by the original one Heff
as
1
ε
lmax
(t) |0i |ψ(0)i − |ψ(t)i ≤h0| Uamp1,pbc
,
(105)
2
3
lmax
(t) |0i |ψ(0)i − |0i |ψ(t)i
Uamp2,pbc

≤ ε,

(106)

for an arbitrary initial state |ψ(0)i ∈ H.
We deliver its detailed proof in Appendix B.2, and instead we briefly explain why therefined effective Hamiltonian is valid. The proof relies mainly on the Lieb-Robinson bound,
lmax
lmax
(t),
stated by Theorem 3. Focusing on the amplification protocols Uamp1(t) and Uamp1,pbc
the resulting difference caused by them originates from their actions on the uniform initialstate |almax i |ψ(0)i, as Eq. (72). As we can see from Fig. 2 (b), the Lieb-Robinson boundlmax
dictates that the dynamics during Uamp1(t) is almost closed within −2lmax . l . 2lmax
in the Fourier indices. In contrast, the support of the additional terms of the refined
effective Hamiltonian in Eq. (101) is located at {|li}l for l ' ±4lmax , which is sufficiently
far from the Fourier indices relevant for the dynamics. As a result, the refined protocol
lmax(t) transforms the state |almax i |ψ(0)i in almost the same way as the original
Uamp1,pbc
lmax(t), and hence it also outputs the target state |ψ(t)i as Eq. (105). Discussion
one Uamp1
similar to Section 5.2 ensures that the oblivious amplitude amplification under the refinedlmax
effective Hamiltonian, represented by Uamp2,pbc(t), provides |ψ(t)i with a sufficiently small
error O(ε) as Eq. (106).
6.2.2

Block-encoding of a refined effective Hamiltonian

4lmax
The benefit of employing the refined effective Hamiltonian Heff,pbc
instead of the original
4lmax
is reduction of resources for implementing block-encoding. We construct
one Heff
4lmax
max
i for it. With the help of theand an oracle state |G4l
an oracle unitary gate Oeff
eff
additional terms in Eq. (101), the Hamiltonian is written in the following form;
4lmax=
Heff,pbc

4lmax
max
Add4l
⊗ Hm − HLP
,
m

(107)

|l ⊕ mi hl| .

(108)

X
|m|≤mmax

max
Add4l
m

=

X
l∈D4lmaxmax
The unitary gate Add4l
is a full quantum adder that translates an index l by m modulom8lmax . It can be implemented by O(log lmax ) complexity and elementary gates.
Then, to construct the block-encoding, we prepare four kinds of auxiliary systems
labeled by a, b, c, d. The system a is an na -qubit system prepared for the block-encoding
4lmax
of Hm [See Eq. (98)], and the system b is the one for HLP
having 8lmax dimension [SeeSection 6.1]. The other systems c and d have the dimension 2mmax + 1 and 2 respectively,
yielding the number of qubits ncd ∈ O(log mmax ) (it is a constant here). We define a
4lmax
unitary gate Oeff
by

4lmax
Oeff
= |0i h0|d ⊗

X

4lmax
4lmax
|mi hm|c ⊗ Addm
⊗ Om − |1i h1|d ⊗ OLP
,

(109)

|m|≤mmax
max
where we omit the identity operators. We also provide a state |G4l
i by
eff
√
P √
αm |0id |mic |a4lmax ib |Gm ia + 4lmax ω |1id |0ic |a4lmax ib |0ia
max
√
|G4l
i= m
.
eff
α + 4lmax ω

Accepted in Quantum 2023-03-16, click title to verify. Published under CC-BY 4.0.

(110)
22

In the above formula, we implicitly include ancillary qubits to efficiently implement the
4lmax
oracles [e.g. the comparator for OLP
in Eq. (95)]. The number of such ancillary qubits
is smaller than O(log lmax ), and hence we neglect it below. We can confirm that they4lmax
as follows;
provide the block-encoding of Heff
4lmax
max
max
i=
|G4l
|Oeff
hG4l
eff
eff

4lmax
Heff,pbc

α + 4lmax ω

,

(111)

4lmax
represented by Eqs. (98) and (93).
where we use the block-encoding for Hm and HLP
4lmax
We evaluate the resources for implementing the oracles. The oracle unitary gate Oeff
given by Eq. (109) requires elementary gates at-most O(mmax (log lmax + C) + log lmax ),
4lmax
max
respectively yield O(log lmax ), C,, Om , and OLP
since the unitary operators Add4lmmax
O(log lmax ) gates. On the other hand, the oracle state |G4l
i can be prepared by
eff
max
max
|wid (Gcoef |0ic ) |a4lmax ib |0ia⊗na ,i = G4l
|G4l
eff
eff

(112)

max
and the state |wid defined bywith using the quantum circuit G4l
eff
max
G4l
eff

= Id ⊗

X|mi hm|c ⊗ Ib ⊗ (Gm )a ,(113)|m|≤mmaxr|wid =α|0i +α + 4lmax ω ds

4lmax ω|1i .
α + 4lmax ω d

(114)

max
can be composed of O(mmax C) elementary gates. The cost of
The quantum circuit G4leff
the state preparation |wid is negligible compared to others since it requires only a singlequbit rotation. The cost for preparing |a4lmax ib is at-most O(log lmax ) complexity and ele4lmax
maxi, which embody
and |G4l
mentary gates. To summarize, the cost for the oracles Oeff
eff
4lmax
the refined effective Hamiltonian Heff,pbc via Eq. (111), is at-most O(mmax (C + log lmax ))elementary gates. The additional gates per query, other than the block-encoding, amountsto O(na + log(lmax )).
4lmax
As a result, preparing a unitary circuit corresponding to e−iHeff t with an allowableerror up to O(ε) is summarized as follows;• Number of ancillary qubits;
na (lmax ) = na + O(log lmax + log mmax ).

(115)• Scaling of query complexity q(lmax );(α + lmax ω)t +log(1/ε).
log(e + {(α + lmax ω)t}−1 log(1/ε))(116)• Number of overall gates;
O(q(lmax )mmax {C + na + log lmax })(117)

Since mmax is supposed to be a O(1) constant, we will omit it from the cost in the rest of
the paper.

Accepted in Quantum 2023-03-16, click title to verify. Published under CC-BY 4.0.

23

Trotterization(for H and H(t))
Time-independent H(Qubitization [18])Time-periodic H(t)(Adiabatic, Theorem 1)
Time-periodic H(t)(Generic, Theorem 2)
Time-dependent H(t)(Dyson series [19, 20])

Ancillary qubits

Query complexity

0αt · (αt/ε)1/p , (p-th order)

na + O(1) (Independent of t, ε)log(1/ε)αt + log(e+(αt)−1 log(1/ε))

na + O(log(γt) + log log(1/ε))log(1/ε)αt + log(e+{αt+o(log(1/ε))}−1 log(1/ε))

na + O(log(γ/ω) + log log(ωt/ε))ωt log(ωt/ε)
αt + log(e+{α/ω+o(log(ωt/ε))}−1 log(ωt/ε))na + O(log{(γωt/α + αt)/ε})αt loglog(αt/ε)log(αt/ε)Table 1: Scaling of the computational resources for Hamiltonian simulation. The definitions for several
parameters are a little different between the above studies. In that case, for simplicity, we replace themdby those which have the similar scales (e.g. kH(t)k ↔ α, dt
H(t) ↔ ωγ). We also note that thequery complexity is measured by the oracles, but their definitions depends on the algorithms. See Eq.
(133) [for Trotterization], Eq. (8) [for qubitization], Eq. (24) and Eq. (26) [for our algorithm], andEq. (135) [for the truncated Dyson-series algorithm].

7 Algorithm and Its computational costThis section provides the main result of this paper; we compose the efficient quantum
algorithm for simulating the time-evolved state |ψ(t)i = U (t) |ψ(0)i under the HamiltonianH(t + T ) = H(t). We note that we take two different approaches depending on the timescale of interest. The first case is the adiabatic-like case, in which we are interested in O(1)period dynamics with t/T ∈ O(1). We call it “adiabatic-like” since long-time dynamics
during 0 ≤ t ≤ T under the sufficiently large period T , exemplified by Thouless pumping
and adiabatic state preparation, is a typical target. The second case is the generic longtime case, where we are interested in multiple-period dynamics at t/T ∈ Ω(1). In thatcase, the the period T is not so large and we often consider long-time dynamics at t  T ,exemplified by laser-irradiated materials.

7.1 Adiabatic-like cases
We consider the adiabatic-like cases, where long-time dynamics over O(1)-periods is ofinterest. For simplicity, we first consider the dynamics within one period at t ∈ [0, T ].
The algorithm in this case is composed of the following steps;
1. Determine the truncation order of the Fourier index lmax by Eq. (59), for the given
time t.
4lmax
4lmax
2. Compose two unitary gates ULP(t) and Ueff(t), corresponding to the time
4lmax

evolution operators e−iHLP t and e−iHeff,pbc t by the qubitization technique; We canfind such operators satisfying
4lmax
ULP
(t) |0i⊗na (4lmax ) |li |ψi = |0i⊗na (4lmax ) e−ilωt |li |ψi + O(ε),
4lmax
Ueff
(t) |0ina (4lmax ) |li |ψi = |0i⊗na (4lmax ) e

4lmax−iHeff,pbc
t

(118)

|li |ψi + O(ε), (119)for arbitrary inputs l ∈ D4lmax and |ψi ∈ H.
4lmax
4lmax
3. Execute the amplification protocols with substituting ULP(t) and Ueff(t) re-Accepted in Quantum 2023-03-16, click title to verify. Published under CC-BY 4.0.

24

4lmaxspectively for e−iHLP t and e−iHeff,pbc t [See Section 5];
lmax
4lmax
lmax
4lmax †
4lmax
Uini
,
Ueff
U¯amp1= (Uini
) ULP
lmax
lmax
lmax †lmax
U¯amp2= −U¯amp1
R(U¯amp1) RUamp1
.

(120)
(121)

lmax
4. Apply the unitary operation U¯amp2(t) to the initial state |0i⊗na (4lmax ) |0i |ψ(0)i;lmax
U¯amp2(t) |0i⊗na (4lmax ) |0i |ψ(0)i = |0i⊗na (4lmax ) |0i |ψ(t)i + O(ε).

(122)

Step 1. determines the dimension of the truncated Floquet-Hilbert space as C8lmax ⊗ H
so that it can accurately reproduce |ψ(t)i with the precision 1 − O(ε), as discussed in
Section 4. In Step 2., we employ the qubitization technique to realize the time evolution
4lmax
(t)
in the truncated Floquet-Hilbert space, as discussed in Section 6. We note that ULP
4lmax(t) can exploit a common O(log(lmax ))-qubit ancillary system. It is sufficientand Ueff
to prepare an {na (4lmax )}-qubit ancillary state |0i⊗na (4lmax ) [See Eq. (115)]. Steps 3. and4. execute the amplification of the time-evolved state |ψ(t)i as discussed in Section 5. By
projecting the ancillary state to |0i⊗na (4lmax ) |0i, we succeed in preparing the target state|ψ(t)i with precision 1 − O(ε), where the success probability is 1 − O(ε). These steps do
not change the scaling of the required resource from Eqs. (120) and (121), reflecting thatthe unitary operations Uini and R yield at-most O(log lmax ) gates.
We finally determine the resource for simulating time-periodic Hamiltonians in adiabaticlike regimes. The ancillary system should involve the degree of freedom for Fourier indices{|li}l∈D4lmax and na (4lmax ) qubits for qubitization. The total number of ancillary qubitsis(123)
na (4lmax ) + dlog2 (8lmax )e ∈ na + O(log(γt) + log log(1/ε)).
The numbers of the oracles Om and Gcoef are dominated by the qubitization under the refined effective Hamiltonian. Considering the assumption ωt ∈ O(1), the query complexity
given by Eq. (116), scales aslog(1/ε)log(1/ε)+,log(e + (γt)−1 log(1/ε)) log(e + {(α + γ)t + o(log(1/ε))}−1 log(1/ε))
log(1/ε)∈ O αt +,(124)
log(e + {αt + o(log(1/ε))}−1 log(1/ε))(α + γ)t +where the o(log(1/ε)) term scales as Eq. (29) originating from the scaling of lmax . The
scaling dependent on γ can be absorbed into that of α, due to the relation γ ≤ α. We
summarize the results in Theorem 1 and Table 1. The number of overall elementary gates
there is determined by the query complexity and the set of additional quantum gates per
query, each of which at most yields O(log lmax ) gates. By the O(1)-times repetition of thisprocedure for t ∈ [0, T ], the same result for O(1)-period dynamics is obtained.

7.2 Generic long-time cases
We next consider the other generic cases where we are interested in long-time dynamics
over multiple periods as ωt ∈ Ω(1). In this regime, we take a different strategy from that for
the adiabatic cases; We split the time t by t = (n+δ)T with n ∈ N and δ ∈ [0, 1). Following
this separation, we implement the time-evolution operator U (t) by n-times operation of
U (T ) and single operation of U (δT ). The algorithm is composed of the following steps;

Accepted in Quantum 2023-03-16, click title to verify. Published under CC-BY 4.0.

25

1. Split the time by t = (n + δ)T with n ∈ N and δ ∈ [0, 1). Determine the truncationTorder lmax
by substituting T and ε/n into t and ε of Eq. (59);T
lmax

log(n/ε)log(e + (γT )−1 log(n/ε))
log(ωt/ε)= Θ γ/ω +.
log(e + (γ/ω)−1 log(ωt/ε))∈ Θ γT +

(125)

lmax
2. Compose a unitary gate U¯amp2(T ), which satisfies
T

T
T
lmax
U¯amp2(T ) |0i⊗na (4lmax ) |0i |ψ(0)i = |0i⊗na (4lmax ) |0i |ψ(T )i + O(ε/n),
T

(126)

by the qubitization and the amplitude amplification protocols (Follow Steps 2.-4. ofT
Section 7.1 wih substituting T , ε/n, and lmax
into t, ε, and lmax respectively).
lmax
3. Apply the unitary gate U¯amp2(T ) to the initial state n times, which results inTmax[U¯amp2(T )]n |0i⊗na (4lmax ) |0i |ψ(0)i = |0i⊗na (4lmax ) |0i |ψ(nT )i + O(ε).

lT

T

T

(127)

lmax
T
4. Prepare the unitary gate U¯amp2(δT ) by substituting δT , ε, and lmax
into t, ε, andlmax respectively in Steps 2.-4. of Section 7.1. Applying it once to the above state
results in
T

T
T
lmax
lmax
U¯amp2(δT )[U¯amp2(T )]n |0i⊗na (4lmax ) |0i |ψ(0)i = |0i⊗na (4lmax ) |0i |ψ(t)i + O(ε).(128)TT

for arbitrary initial states |ψ(0)i.
We remark several points in each step. Steps 1. and 2. are executed to apply U (T ),giving the time-evolution over one period T . Here, we set the acceptable error to O(ε/n)
so that we can obtain the time-evolved state |ψ(nT )i with an error up to O(ε) after thelT
n-times repetition in Step 3. The cost for the unitary gate U¯ max (T ) is dominated only by
4lT



amp2
4lT



max
that of Ueff max , which reproduces the time evolution exp −iHeff,pbc
T by the qubitization.
4lT
max

This comes from the fact that the time evolution e−iHLP
4lTmax

e−iHLP

T=

X

T

is trivial due toeilωT |li hl| ⊗ I = 1,

(129)

T
l∈D4lmax

4lT

and we do not need ULPmax in Eq. (120). As a result, the query complexity for n-times
lT
implementation of U¯ max (T ) during Step 3. is proportional to
amp2

log(n/ε)T ω)T }−1 log(n/ε))log(e + {(α + lmaxωt log(ωt/ε)∈ O αt +.
log(e + {α/ω + o(log(ωt/ε))}−1 log(ωt/ε))Tn (α + lmax
ω)T +



(130)

We use the fact n ∈ O(ωt) in the above relation. Step 4. realizes the remaining micromotion U (δT ) for the duration δT with an error up to O(ε). We remark that the choice ofTthe truncation order lmax
so far is sufficient to achieve the precision 1 − O(ε) for U (δT )TT
δT
due to lmax > lmax . This implies that we can reuse the ancillary state |0i⊗na (4lmax ) for the
Accepted in Quantum 2023-03-16, click title to verify. Published under CC-BY 4.0.

26

qubitization approach to the time-evolution U (δT ) in Step 4. The cost for implementinglT
U¯ max (δT ) once is obtained by setting t = δT in adiabatic-like cases. It is smaller thanamp2

the cost for implementing the time-evolution U (T ), and does not affect the scaling of thecomputational resource.
Finally, we provide the computational resource for time-periodic Hamiltonian dynamics
in generic long-time regimes. The number of ancillary qubits is given by na (4lmax ) +T )e, which is bounded by
dlog2 (8lmax
na + O(log(γ/ω) + log log(ωt/ε)).

(131)

The query complexity is dominated by Eq. (130). We summarize them in Theorem 2 andTable 1.

7.3 Comparison with other algorithms
Let us compare our algorithm on time-periodic Hamiltonians with other quantum algorithms for Hamiltonian simulation, based on Table 1. We pick up the qubitization
technique [18] for time-independent Hamiltonian H and the truncated Dyson-series algorithm [19, 20] for generic time-dependent Hamiltonian H(t), whose resources have thebest scaling in t and 1/ε as far as we know. We also consider the standard way, Trotterization, which covers generic time-independent, time-periodic, and generic time-dependent
Hamiltonians.
Let us first compare our algorithm with Trotterization. When the Hamiltonian can bePdivided into H(t) = Γr=1 Hr (t), where every term in Hr (t) commutes with one anotherat every time, the first-order Trotterization roughly approximates the time-evolution U (t)
by
!#
" Γ
!
Rt
M
−1 Y
Y
(αt)2
−i l+1 dt0 Hr (t0 )tl+OU (t) =
Te
,
(132)
M
l=0 r=1
with M partitions of the time t as tl = lt/M [2]. Reflecting that the Trotterization error
polynomially decreases in the partition number M , we need queries to a layered quantum
circuit
!
R
Γ
Y

−i

Te

tl+1
dt0 Hr (t0 )
tl

(133)

r=1

for O (αt)2 /ε times to achieve the allowable error ε. While the coefficient α and thepowers in t and 1/ε can be improved by higher-order Trotterization [12], its resource
increases polynomially in 1/ε. Our algorithm has better scaling of elementary gates thanTrotterization in that the resource increases logarithmically in 1/ε.
The comparison with the qubitization and the truncated Dyson-series algorithm isinstructive for evaluating the efficiency of our algorithm due to the inclusion relation,
[Set of time-independent H] ⊂ [Set of time-periodic H(t) = H(t + T )]⊂ [Set of generic time-dependent H(t)].

(134)

The qubitization technique achieves the least number of ancillary qubits in a sense that
it is independent of t and 1/ε. It also has the best query complexity, in which the optimal scaling both in t and 1/ε appears in an additive way as Eq. (14). The number ofancillary qubits and the query complexity for the qubitization provide the best bound fortime-periodic Hamiltonians, while it is nontrivial whether or not it is actually achievable.
Comparison with the truncated Dyson-series algorithm tells us how efficiently we deal
Accepted in Quantum 2023-03-16, click title to verify. Published under CC-BY 4.0.

27

with the time-dependent Schrödinger equation. The cost reduction compared to them canbe interpreted as improvement of efficiency due to the time-periodicity.
Before going to the comparison respectively for the adiabatic-like regime and thegeneric long-time regime, we remark several points in common. First, we replace someparameters by those which have similar scales to simply compare these algorithms. For
instance, the parameters α and γ respectively give the energy scales of the whole Hamiltonian H(t) and the time-dependent terms H(t) − H0 according to Eqs. (25) and (23). Wedsubstitute α and ωγ respectively for supt (kH(t)k) and supt ( dt
H(t) ), which are characteristic values in the truncated Dyson-series algorithm.
The second point is the complexity of the oracles themselves. In the qubitization
technique, the query complexity is measured by the oracles to a static Hamiltonian H as
Eq. (8). In contrast, the truncated Dyson-series algorithm for generic cases employs theoracleHam−O =

MX

|li hl| ⊗ O(tl ),

(135)

l=1

where l labels the discretized time tl with the partition number M ∈ O(1/ε) [19]. The
oracle Ham−O includes multiple implementation of O(tl ), which is an oracle for a static
instantaneous Hamiltonian H(tl ) (some specific cases such as LCUs and sparse-accessmatrices with integrable time-dependency can be simplified [20]). Our algorithm for timeperiodic Hamiltonians uses the oracles {Om , |Gm i}, which gives block-encoding of each
Fourier component Hm as Eq. (98), and the oracle for the coefficients, Gcoeff given by Eq.
(99). The query complexity in our algorithm is roughly measured by the oracle for a static
operator Hm since the latter one Gcoeff , a quantum gate on at-most O(1) qubits, is usually
negligible. Therefore, note that our algorithm adopts essentially the same measure for thequery complexity as the query complexity, while the truncated Dyson-series algorithm
counts it by rather complicated oracles involving discretized time.
7.3.1

Adiabatic-like cases

We assess the cost for simulating time-periodic Hamiltonians in the adiabatic-like cases
based on Table 1. First, we compare the number of ancillary qubits with those for otheralgorithms. What should be noted in our algorithm is its scaling in the inverse error 1/ε.
Our algorithm requires O(log log(1/ε)) additional qubits, whose number lies just in the
middle of that for the qubitization (ε-independent) and that for the truncated Dyson-series
algorithm [O(log(1/ε))]. Importantly, the reduction compared to the latter one can be
attributed to the faster convergence of our formalism based on the Fourier indices. The
error arising from time discretization polynomially decays in the partition number M ;as a result, the truncated Dyson-series algorithm requires at-least O(log(1/ε)) ancillary
qubits via the oracle Eq. (135). In contrast, we introduce thetruncation order lmax for thelmaxFloquet-Hilbert space. The error by this cutoff scales as O (γt/lmax ), whose decay
is faster than the exponential function e−O(lmax ) . This leads to the reduction of ancillary
qubits from O(log(1/ε)) to O(log log(1/ε)).
Next, we compare the query complexity. When we fix the inverse error 1/ε and increasetime t, it scales as O(αt). The linear increase in time t implies the optimality of ouralgorithm in time. When we consider 1/ε → ∞ for given time, the query complexity
scales aslog(1/ε),(136)
log log log(1/ε)

Accepted in Quantum 2023-03-16, click title to verify. Published under CC-BY 4.0.

28

where we substitute the form of the o(log(1/ε)) term, Eq. (29), into Eq. (28). This scalingis nearly-optimal in 1/ε. Importantly, their contributions affect the query complexity inan additive way as αt + o(log(1/ε)). This means that our algorithm deals with timeperiodic systems with the cost sufficiently close to qubitization [See Eq. (14)], and saves
resource compared to the truncated Dyson-series algorithm requiring the multiplicative
query complexity αt × o(log(αt/ε)). As long as we suppose polynomial accuracy ε ∈1/poly (N ), implying log(1/ε) . αt, the query complexity can achieve the best scaling fortime-independent systems.
While we hereby provide comparison with the truncated Dyson-series algorithm, ouralgorithm outperforms other time-dependent Hamiltonian simulation algorithms [21, 22,
23, 24] in the presence of time-periodicity. For instance, Ref. [24] can achieve the additive
query complexity when the second derivative of H(t) in t vanishes (i.e., linear dependence),
but time-periodicity prohibits vanishing derivatives. For time-periodic Hamiltonian simulation, the additive query complexity close to that of the qubitization is unique to ouralgorithm. Any of the other algorithms using the discretized time requires at-least O(1/ε)
degrees of freedom, while we need O(log(1/ε)) for the Fourier indices. Thus, our algorithm
always saves the number of ancillary qubits.
7.3.2

Generic long-time cases

Here, we evaluate the computational resources for generic long-time regimes. We supposethat the frequency ω is constant and hence small compared to the whole-system energy
scales α and γ, which usually increase polynomially in the system size N . This assumption
comes from the fact that ω typically represents the frequency of external drives (e.g. light),
which is size-independent. More importantly, when the frequency ω is poly (N ) so thatit becomes comparable to α and γ, the dynamics can be efficiently simulated by themethods for time-independent Hamiltonians with the help of high-frequency expansions[59, 60, 61, 62, 63], as we prove in Appendix
D. Therefore, it is sufficient to consider thecomputational resources under ω ∈ O N 0 and α, γ ∈ poly (N ).
The number of ancillary qubits is similar to that for the adiabatic cases. We notethat its scaling in the time [O(log log t)] overwhelms that of the truncated Dyson-series
algorithm [O(log t)] in addition to the inverse error 1/ε. The number of ancillary qubits
for time-periodic Hamiltonians lie between those for time-independent and time-dependent
Hamiltonians in terms of both t and 1/ε.
We discuss the query complexity given by Theorem 2. Under the fixed allowable errorε, the scaling of the query complexity in sufficiently large time t is given by
αt + ωtlog(ωt/ε)≤ {α + ω log(1/ε)}t + ωt log(ωt).
log log log(ωt/ε)(137)

Due to the second term in the right hand side, the query complexity has at least nearlyoptimal dependence in t, which involves logarithmiccorrection log(ωt). However, reflecting
the assumptions α ∈ poly (N ) and ω ∈ O N 0 , the O(ωt log(ωt)) term can be nonα/ωnegligible compared to O(αt) only when the time reaches t ∼ e /ω ∈ O epoly(N ) T .
Therefore, unless we consider unpractical time scales t ∈ eΘ(N ) , this nearly-optimal scaling
is formal, and the query complexity actually increases linearly in time t, showing the
optimal scaling. On the other hand, when we focus on the scaling in the inverse error 1/ε,
the query complexity scales aslog(1/ε)ωt.
(138)
log log log(1/ε)
Accepted in Quantum 2023-03-16, click title to verify. Published under CC-BY 4.0.

29

Thus, our algorithm achieves nearly-optimal scaling in 1/ε with slight deviation from theoptimal dependence by log log(1/ε)/ log log log(1/ε).
When it comes to the combined scaling of the query complexity in t and 1/ε, we notethat it has an additive form as
αt + ωt × o(log(1/ε)),(139)

where we neglect O(ωt log(ωt)). Although linear increase in t couples with logarithmic
increase in 1/ε, the coefficient ωt ∈ O N 0 t is much smaller than the whole energy-timescale αt ∈ poly (N ) t. As a consequence, also in generic long-time cases, our algorithm
achieves the query complexity close to that of qubitization [See Eq. (14)] and saves muchcost compared to the truncated Dyson-series algorithm, where the query complexity scales
as poly (N ) t × o(log(αt/ε)).

8 Illustrative examples
In this section, we briefly discuss some potential applications of the algorithm. We expectthat it can be applied to nonequilibrium quantum many-body phenomena, which are often
of interest in condensed matter physics and quantum chemistry. In terms of quantum
computation, it will offer an efficient protocol for adiabatic state preparation, which canbe applied to quantum phase estimation for instance. We suggest the simplest examples
for both applications below.

8.1 Nonequilibrium quantum many-body phenomena
The first application is to simulate nonequilibrium dynamics of periodically-driven quantum materials. Optical responses are typical but of great interest both in condensed matterphysics and quantum chemistry. We pick up an N -site Fermi-Hubbard model under laser
light as the simplest case;H(t) = HHub + Hext (t),HHub =

X X(140)

k n̂kσ + U

Xn̂x↑ n̂x↓ .

(141)

x

k σ=↑,↓

Here, n̂kσ = ĉ†kσ ĉkσ and n̂xσ = ĉ†xσ ĉxσ are number operators of electrons in the momentum
and real spaces respectively, generated by fermionic annihilation operators ĉkσ and ĉxσ .
The time-periodic term Hext (t) represents the coupling with light. When we shine
linearly-polarized light with the frequency Ω, it is given by
Hext (t) = sin(Ωt)

X

Vx n̂xσ ,

(142)

x,σ

which results in T = 2π/Ω and mmax = 1. The Fourier components Hm are given byP
H0 = HHub and H±1 = (±i) x,σ (Vx /2)n̂xσ . To evaluate the cost of simulating |ψ(t)i,
we compose the oracles for them. We employ a unitary operation, called fermionic fast
Fourier transform (FFFT) [64, 65, 66], which transforms the basis in the momentum space
ki to that in the real space xi as n̂ki σ = FFFT† n̂xi σ FFFT. We map the fermionic system
to a spin system by Jordan-Wigner transformation as n̂xσ = (1 + Zxσ )/2 (“xσ” denotes anPindex for qubits). By neglecting constant terms and a conserved particle number xσ n̂xσ ,Accepted in Quantum 2023-03-16, click title to verify. Published under CC-BY 4.0.

30

the block-encoding for {Hm } can be constructed by the technique for LCUs, as Eq. (10).
Assuming k , U, Vx ≥ 0 without loss of generality, this leads to
O0 = |0i h0| ⊗X|x, σi hx, σ| ⊗ FFFT† Zxσ FFFTx,σ

+ |1i h1| ⊗X|x, σi hx, σ| ⊗ Zx↑ Zx↓ ,

(143)

x,σ

|G0 i =

X r εxα0x,σO±1 = I ⊗

X

s

|0i +

!

U
|1i |x, σi ,2α0|x, σi hx, σ| ⊗ (±iZxσ ),

(144)
(145)

x,σ

s

|G±1 i = |0i

X
x,σ

Vx
|x, σi ,
2α±1

(146)

with α0 = 2 x x + U N and α±1 = x Vx . In the above oracles, the unitary gate O0
requires much cost due to O(1)-times usage of FFFT, which can be implemented with
at-most O(N log N )-depth quantum circuits composed of adjacent two-qubit gates [66].
The depth needed for each oracle is at-most O(N log N ).
The number of ancillary qubits for them to express {|0i |x, σi , |1i |x, σi} amounts tona ∈ O(log N ). The energy scales of the whole Hamiltonian, α, determined by Eq. (25)asXX
x + U N + 2
Vx .
(147)
α=2
P

P

x

x

Although it is difficult to obtain γ from its definition Eq. (23), we can easily obtain its
upper bound, which results in
γ≤

X
m=±1

kHm k =

X

Vx .

(148)

x

It is adequate for determining the computational resource, since α is rather dominant
in the query complexity and γ appears as O(log γ) in the number of ancillary qubits.
When we define the characteristic local energy scale by αloc = max(k , U, Vx ), they areapproximately described by α, γ ∈ O(αloc N ).
Finally, if we are interested in the time-evolved state |ψ(t)i over multiple periods, thefollowing resources are required to achieve the precision 1 − O(ε);• Number of ancillary qubits;
O(log(αloc N/Ω) + log log(Ωt/ε)).

(149)

O({αloc N t + Ωt log(Ωt/ε)} N log N ).

(150)

• Overall gate complexity;

The above results are based on Eqs. (130) and (131), where we neglect size-independent
additional gates for each query to the oracles. According to classical numerical calculations[67, 68, 69], the above model is expected to host high-harmonic generation, where intense
oscillation with the frequency nΩ (n = 2, 3, . . .) arises in response to laser light with thefrequency Ω. Our algorithm allows to efficiently identify such a nontrivial nonequilibrium phenomenon with guaranteed accuracy by simulating the Fourier spectrum of someobservables hψ(t)|O|ψ(t)i (e.g. electric current) [70].
Accepted in Quantum 2023-03-16, click title to verify. Published under CC-BY 4.0.

31

Another interesting example is a discrete time crystal as a phase of matter inherent
in nonequilibrium, where time-translation symmetry is spontaneously broken [28, 29, 30].
When we choose Hext (t) by uniform circularly-polarized ac field represented by
Hext (t) = VXeiΩt ĉ†x↑ ĉx↓ + e−iΩt ĉ†x↓ ĉx↑ ,

(151)

x

the Fermi-Hubbard model becomes a potential platform for a time-crystalline phenomenon
protected by Floquet dynamical symmetry [71]. Its signature can be detected by subharmonic oscillations of spatio-temporal correlation functions and local observables. These
values are both efficiently computed via the time-evolution operator U (t) by our algorithm,
requiring the computational resources similar to Eqs. (149) and (150).
As nonequilibrium systems dominated by time-periodic Hamiltonians have been vigorously explored as Floquet systems, our algorithm will cover various phenomena. We also
note that our algorithm is extended to time-periodic Hamiltonians with exponentiallydecaying Fourier components kHm k . e−O(|m|) . Since we often face at situations where
high-frequency components of H(t) rapidly diminish, our result will be useful for a variety
of nonequilibrium phenomena in condensed matter physics and quantum chemistry, otherthan the above examples (See also Appendix C.5).

8.2 Adiabatic state preparation
Adiabatic state preparation is a protocol to obtain a preferable quantum state by adiabatically evolving quantum systems with time-dependent Hamiltonians H(t). While it has
been originally developed in the context of adiabatic quantum computation relying on the
adiabatic theorems of quantum dynamics [42], it can be exploited also on circuit-based
quantum computers by mimicking the adiabatic dynamics under H(t). One of typical aims
is to prepare initial states required for quantum simulation [41, 6, 58]. Here, we provide
the simplest application of our algorithm in this field.
We prepare a time-periodic Hamiltonian H(t), which continuously connects two different time-independent Hamiltonians H̄0 and H̄1 . A certain eigenstate of H̄0 , denotedby |ψ̄0 i, is supposed to be easily prepared on quantum circuits, and we assume that the
eigenstate of H̄1 , which is continuously connected to |ψ̄0 i, corresponds to the target state|ψ̄1 i. As the simplest case, we organize such a time-periodic Hamiltonian by
H(t) = H̄0 (1 − sin ωt) + H̄1 sin ωt,(152)with satisfies H(0) = H̄0 and H(T /4) = H̄1 . Our algorithm for the adiabatic-like cases
in Section 7.1 enables us to efficiently execute the adiabatic state preparation. Since theFourier components of H(t) is simply given by H0 = H̄0 and H±1 = ±(H̄0 − H̄1 )/2i withmmax = 1, the following oracles are necessary;
O0 = I ⊗ Ō0 ,|G0 i = |0i |Ḡ0 i ,O±1 = |0i h0| ⊗ (∓iŌ0 ) + |1i h1| ⊗ (±iŌ1 ),√√ᾱ0 |0i |Ḡ0 i + ᾱ1 |1i |Ḡ1 i√|G±1 i =.
ᾱ0 + ᾱ1

(153)
(154)
(155)

Here, (Ō0 , |Ḡ0 i , ᾱ0 ) and (Ō1 , |Ḡ1 i , ᾱ1 ) respectively provides the block-encoding of thetime-independent Hamiltonians H̄0 and H̄1 as Eq. (8). The oracles O±1 and |G±1 i gives
the block-encoding of H±1 as hG±1 |O±1 |G±1 i = 2H±1 /(ᾱ0 + ᾱ1 ). According to the standard adiabatic theorem [42], the duration required for approximating the target state |ψ̄1 i
Accepted in Quantum 2023-03-16, click title to verify. Published under CC-BY 4.0.

32

by the adiabatically-evolved state U (T /4) |ψ̄0 i with accuracy 1 − O(ε) roughly amounts
to T & H̄0 − H̄1 /ε∆2 , where ∆ denotes the minimal gap upon the instantaneous eigenstate of H(t) continuously connecting |ψ̄0 i and |ψ̄1 i. Therefore, in the adiabatic state
preparation based on our algorithm, the query complexity counted by the oracles of thetwo static Hamiltonians H̄0 and H̄1 amounts toᾱ0 + ᾱ1H̄0 − H̄1 + o(log(1/ε)).
ε∆2

(156)

Although the second term is buried by the first one polynomially increasing in 1/ε, thecost has better scaling in 1/ε compared to the cases where similar schedules for H(t) aretackled with the Trotterization [∼ (1/ε)1+2/p ] or the truncated Dyson-series algorithm[∼ (1/ε) log(1/ε)].
In adiabatic quantum computation, the required time in 1/ε can be improved via moresophisticated scheduling, such as local adiabatic interpolation [72], boundary cancellation[73, 74] and quasi-adiabatic processes [75]. In particular, the last one has achieved thequery complexity poly-logarithmic in 1/ε with the help of the truncated Dyson-series algorithm. If we can find time-periodic Hamiltonians which are consistent with such sophisticated schedules, it will offer a more efficient protocol for the adiabatic state preparation.

9 Discussion and Conclusion
We conclude our paper with summarizing the results. In this paper, we focus on timedependent systems with time-periodicity, and organize an efficient implementation of theirtime-evolution operators. Once we prepare the oracles which embeds each Fourier component and each coefficient of the Hamiltonian, a series of unitary operations on the truncatedFloquet-Hilbert space extract the time-evolved state |ψ(t)i with an allowable error 1/ε,
whose query complexity are both optimal or nearly-optimal in t and 1/ε. In addition, it
has an additive scaling αt + o(log(1/ε)) [adiabatic-like regime] or αt + ωt × o(log(1/ε))[generic long-time regime], and hence achieves the cost sufficiently close to that of the bestalgorithm for time-independent systems [18], despite the existence of time-dependency.
As exemplified by nonequilibrium quantum many-body phenomena and adiabatic state
preparation, our algorithm will contribute to pioneering applications of quantum computers for various aims, in condensed matter physics, quantum chemistry, and quantum
computation.
We finally discuss some potential future directions of our results. The first one is to
explore efficient implementation of meaningful functions
in time-periodic
systems, other
 R

t
0
0
than the time-evolution operator U (t) = T exp −i 0 H(t )dt . In our algorithm, we
lmax

exploit the qubitization technique to implement the exponential function, e−iHeff t , inlmax
the effective Hamiltonian Heff
. As the qubitization technique allows to efficiently apply
various polynomial functions, we expect that our algorithm can be extended for various
aims other than the unitary time evolution discussed here. For instance, our discussion on
the Lieb-Robinson bound in Section 4 is valid also when a time-periodic Hamiltonian H(t)is non-hermitian. This suggests that our algorithm may be available also for solving timeperiodic linear differential equations [76], exemplified by dissipative quantum may-bodysystems [77, 78]. While we do not expect exponential speedup due to the non-unitarity
in general, it will offer an efficient way both in time and desirable accuracy. It would
be important to clarify what kind of function is useful in time-periodic systems, how the
Lieb-Robinson bound and the amplitude amplification should be modified, and then how
the computational cost is affected.
Accepted in Quantum 2023-03-16, click title to verify. Published under CC-BY 4.0.

33

Another significant direction is to seek for useful tasks that can be efficiently tackledwith time-periodic Hamiltonians. As we have shown throughout the paper, time-periodic
Hamiltonians can be simulated more efficiently generic time-dependent Hamiltonians withcomputational resources close to those for time-independent Hamiltonians. This meansthat quantum tasks which inevitably requires time-dependent operations, such as adiabatic
state preparation, can be optimized by tuning their schedules in a time-periodic way. It will
be important to clarify what kind of tasks can be addressed by time-periodic Hamiltonians
and how our algorithm provides better scaling in 1/ε for their costs, which we leave for
future work.

Acknowledgment
We thank T. Kuwahara for fruitful discussion, especially on the Lieb-Robinson boundpresented in Section 4.1. This work is supported by MEXT Quantum Leap Flagship
Program (MEXTQLEAP) Grant No. JPMXS0118067394, JPMXS0120319794, and JST
COI-NEXT program Grant No. JPMJPF2014.

References[1] Richard P. Feynman. “Simulating Physics with Computers”. Int. J. Theor. Physics
21, 467 (1982).
[2] S Lloyd. “Universal Quantum Simulators”. Science 273, 1073–1078 (1996).
[3] Adam Smith, M S Kim, Frank Pollmann, and Johannes Knolle. “Simulating quantum many-body dynamics on a current digital quantum computer”. npj Quantum
Information 5, 1–13 (2019).
[4] Frank Arute et al. “Observation of separated dynamics of charge and spin in the
Fermi-Hubbard model” (2020). arXiv:2010.07965.
[5] A. Yu. Kitaev. “Quantum measurements and the Abelian Stabilizer Problem” (1995).
arXiv:quant-ph/9511026.
[6] Jiangfeng Du, Nanyang Xu, Xinhua Peng, Pengfei Wang, Sanfeng Wu, and Dawei Lu.
“NMR implementation of a molecular hydrogen quantum simulation with adiabatic
state preparation”. Phys. Rev. Lett. 104, 030502 (2010).
[7] B P Lanyon et al. “Towards quantum chemistry on a quantum computer”. Nat.
Chem. 2, 106–111 (2010).
[8] P. J. J O’Malley et al. “Scalable Quantum Simulation of Molecular Energies”. Phys.
Rev. X 6, 031007 (2016).
[9] Daniel S Abrams and Seth Lloyd. “Simulation of Many-Body Fermi Systems on a
Universal Quantum Computer”. Phys. Rev. Lett. 79, 2586–2589 (1997).
[10] A T Sornborger and E D Stewart. “Higher-order methods for simulations on quantum
computers”. Phys. Rev. A 60, 1956–1965 (1999).
[11] Earl Campbell. “Random Compiler for Fast Hamiltonian Simulation”. Phys. Rev.
Lett. 123, 070503 (2019).
[12] Andrew M. Childs, Yuan Su, Minh C. Tran, Nathan Wiebe, and Shuchen Zhu. “Theory of Trotter Error with Commutator Scaling”. Phys. Rev. X 11, 011020 (2021).

Accepted in Quantum 2023-03-16, click title to verify. Published under CC-BY 4.0.

34

[13] Andrew M Childs and Nathan Wiebe. “Hamiltonian simulation using linear combinations of unitary operations”. Quantum Information & Computation 12, 901 (2012).
[14] Dominic W Berry, Andrew M Childs, and Robin Kothari. “Hamiltonian Simulation with Nearly Optimal Dependence on all Parameters”. In 2015 IEEE 56th
Annual Symposium on Foundations of Computer Science. Pages 792–809. ieeexplore.ieee.org (2015).
[15] Dominic W Berry, Andrew M Childs, Richard Cleve, Robin Kothari, and Rolando D
Somma. “Simulating Hamiltonian dynamics with a truncated Taylor series”. Phys.
Rev. Lett. 114, 090502 (2015).
[16] Dominic W Berry, Andrew M Childs, Richard Cleve, Robin Kothari, and Rolando D
Somma. “Exponential improvement in precision for simulating sparse Hamiltonians”.
Forum of Mathematics, Sigma 5, E8 (2017).
[17] Guang Hao Low and Isaac L Chuang. “Optimal Hamiltonian Simulation by Quantum
Signal Processing”. Phys. Rev. Lett. 118, 010501 (2017).
[18] Guang Hao Low and Isaac L Chuang. “Hamiltonian simulation by qubitization”.
Quantum 3, 163 (2019).
[19] Guang Hao Low and Nathan Wiebe. “Hamiltonian simulation in the interaction
picture” (2018). arXiv:1805.00675.
[20] Mária Kieferová, Artur Scherer, and Dominic W Berry. “Simulating the dynamics
of time-dependent Hamiltonians with a truncated Dyson series”. Phys. Rev. A 99,
042314 (2019).
[21] Dominic W. Berry, Andrew M. Childs, Yuan Su, Xin Wang, and Nathan Wiebe.
“Time-dependent Hamiltonian simulation with L1 -norm scaling”. Quantum 4,
254 (2020).
[22] Jeongwan Haah, Matthew B Hastings, Robin Kothari, and Guang Hao Low. “Quantum algorithm for simulating real time evolution of lattice Hamiltonians”. SIAM J.
Comput.Pages FOCS18–250–FOCS18–284 (2021).
[23] Yi-Hsiang Chen, Amir Kalev, and Itay Hen. “Quantum Algorithm for TimeDependent Hamiltonian Simulation by Permutation Expansion”. PRX Quantum 2,030342 (2021).
[24] Jacob Watkins, Nathan Wiebe, Alessandro Roggero, and Dean Lee. “Timedependent Hamiltonian Simulation Using Discrete Clock Constructions” (2022).
arXiv:2203.11353.
[25] Takuya Kitagawa, Erez Berg, Mark Rudner, and Eugene Demler. “Topological characterization of periodically driven quantum systems”. Phys. Rev. B 82, 235114 (2010).
[26] Mark S Rudner, Netanel H Lindner, Erez Berg, and Michael Levin. “Anomalousedge states and the bulk-edge correspondence for periodically driven two-dimensional
systems”. Physical Review X 3, 031005 (2013).
[27] Fenner Harper, Rahul Roy, Mark S. Rudner, and S.L. Sondhi. “Topology and Broken
Symmetry in Floquet Systems”. Annual Review of Condensed Matter Physics 11,345–368 (2020).
[28] Vedika Khemani, Achilleas Lazarides, Roderich Moessner, and S L Sondhi. “Phase
Structure of Driven Quantum Systems”. Phys. Rev. Lett. 116, 250401 (2016).
[29] Dominic V Else, Bela Bauer, and Chetan Nayak. “Floquet Time Crystals”. Phys.
Rev. Lett. 117, 090402 (2016).
Accepted in Quantum 2023-03-16, click title to verify. Published under CC-BY 4.0.

35

[30] Vedika Khemani, Roderich Moessner, and S L Sondhi. “A Brief History of Time
Crystals” (2019). arXiv:1910.10745.
[31] Feng Mei, Qihao Guo, Ya-Fei Yu, Liantuan Xiao, Shi-Liang Zhu, and Suotang Jia.
“Digital Simulation of Topological Matter on Programmable Quantum Processors”.
Phys. Rev. Lett. 125, 160503 (2020).
[32] J Randall et al. “Many-body-localized discrete time crystal with a programmablespin-based quantum simulator”. Science 374, 1474 (2021).
[33] Xiao Mi et al. “Time-crystalline eigenstate order on a quantum processor”. Nature
601, 531–536 (2022).
[34] R M Potvliege and R Shakeshaft. “Multiphoton processes in an intense laser field:
Harmonic generation and total ionization rates for atomic hydrogen”. Phys. Rev. A
40, 3061–3079 (1989).
[35] F H M Faisal and J Z Kamiński. “Floquet-Bloch theory of high-harmonic generationin periodic structures”. Phys. Rev. A 56, 748–762 (1997).
[36] Takashi Oka and Hideo Aoki. “Photovoltaic Hall effect in graphene”. Phys. Rev. B
79, 081406 (2009).
[37] Takashi Oka and Sota Kitamura. “Floquet Engineering of Quantum Materials”.
Annual Review of Condensed Matter Physics 10, 387–408 (2019).
[38] D J Thouless. “Quantization of particle transport”. Phys. Rev. B 27, 6083 (1983).
[39] M Lohse, C Schweizer, O Zilberberg, M Aidelsburger, and I Bloch. “A Thoulessquantum pump with ultracold bosonic atoms in an optical superlattice”. Nat. Phys.
12, 350–354 (2015).
[40] Shuta Nakajima, Takafumi Tomita, Shintaro Taie, Tomohiro Ichinose, Hideki Ozawa,
Lei Wang, Matthias Troyer, and Yoshiro Takahashi. “Topological Thouless pumpingof ultracold fermions”. Nat. Phys. 12, 296–300 (2016).
[41] Alán Aspuru-Guzik, Anthony D Dutoi, Peter J Love, and Martin Head-Gordon. “Simulated quantum computation of molecular energies”. Science 309, 1704–1707 (2005).
[42] Tameem Albash and Daniel A. Lidar. “Adiabatic quantum computation”. Rev. Mod.
Phys. 90, 015002 (2018).
[43] Elliott H Lieb and Derek W Robinson. “The finite group velocity of quantum spin
systems”. Commun. Math. Phys. 28, 251–257 (1972).
[44] Hideo Sambe. “Steady States and Quasienergies of a Quantum-Mechanical System
in an Oscillating Field”. Phys. Rev. A 7, 2203–2213 (1973).
[45] T O Levante, M Baldus, B H Meier, and R R Ernst. “Formalized quantum mechanical
Floquet theory and its application to sample spinning in nuclear magnetic resonance”.
Mol. Phys. 86, 1195–1212 (1995).
[46] Ryan Babbush, Craig Gidney, Dominic W Berry, Nathan Wiebe, Jarrod McClean,
Alexandru Paler, Austin Fowler, and Hartmut Neven. “Encoding Electronic Spectrain Quantum Circuits with Linear T Complexity”. Phys. Rev. X 8, 041015 (2018).
[47] András Gilyén, Yuan Su, Guang Hao Low, and Nathan Wiebe. “Quantum singularvalue transformation and beyond: exponential improvements for quantum matrix
arithmetics”. Proceedings of the 51st Annual ACM SIGACT Symposium on Theoryof Computing, 193–204 (2019).

Accepted in Quantum 2023-03-16, click title to verify. Published under CC-BY 4.0.

36

[48] John M. Martyn, Zane M. Rossi, Andrew K. Tan, and Isaac L. Chuang. “Grand
Unification of Quantum Algorithms”. PRX Quantum 2, 040203 (2021).
[49] Lov K. Grover. “Quantum Mechanics Helps in Searching for a Needle in a Haystack”.
Phys. Rev. Lett. 79, 325–328 (1997).
[50] Zongping Gong and Ryusuke Hamazaki. “Bounds in Nonequilibrium Quantum Dynamics”. International Journal of Modern Physics B 36, 223007 (2022).
[51] Robert M Gray. “Toeplitz and Circulant Matrices: A Review”. Foundations andTrends® in Communications and Information Theory 2, 155–239 (2006).
[52] R M Corless, G H Gonnet, D E G Hare, D J Jeffrey, and D E Knuth. “On theLambertW function”. Adv. Comput. Math. 5, 329–359 (1996).
[53] Steven A Cuccaro, Thomas G Draper, Samuel A Kutin, and David Petrie Moulton.
“A new quantum ripple-carry addition circuit” (2004). arXiv:quant-ph/0410184.
[54] Lov Grover and Terry Rudolph. “Creating superpositions that correspond to efficiently integrable probability distributions” (2002). arXiv:quant-ph/0208112.
[55] Aram W. Harrow, Avinatan Hassidim, and Seth Lloyd. “Quantum algorithm forlinear systems of equations”. Phys. Rev. Lett. 103, 150502 (2009).
[56] L K Grover. “Synthesis of quantum superpositions by quantum computation”. Phys.
Rev. Lett. 85, 1334–1337 (2000).
[57] Yuval R Sanders, Guang Hao Low, Artur Scherer, and Dominic W Berry. “BlackBox Quantum State Preparation without Arithmetic”. Phys. Rev. Lett. 122,
020502 (2019).
[58] Arthur G Rattew and Bálint Koczor. “Preparing Arbitrary Continuous Functions in
Quantum Registers With Logarithmic Complexity” (2022). arXiv:2205.00519.
[59] Dmitry A Abanin, Wojciech De Roeck, and François Huveneers. “Exponentially
Slow Heating in Periodically Driven Many-Body Systems”. Phys. Rev. Lett. 115,
256803 (2015).
[60] Tomotaka Kuwahara, Takashi Mori, and Keiji Saito. “Floquet–Magnus theory andgeneric transient dynamics in periodically driven many-body quantum systems”. Ann.
Phys. 367, 96–124 (2016).
[61] Takashi Mori, Tomotaka Kuwahara, and Keiji Saito. “Rigorous Bound on Energy
Absorption and Generic Relaxation in Periodically Driven Quantum Systems”. Phys.
Rev. Lett. 116, 120401 (2016).
[62] Dmitry A Abanin, Wojciech De Roeck, Wen Wei Ho, and François Huveneers. “Effective Hamiltonians, prethermalization, and slow energy absorption in periodicallydriven many-body systems”. Phys. Rev. B 95, 014112 (2017).
[63] Dmitry Abanin, Wojciech De Roeck, Wen Wei Ho, and François Huveneers. “A
Rigorous Theory of Many-Body Prethermalization for Periodically Driven and Closed
Quantum Systems”. Commun. Math. Phys. 354, 809–827 (2017).
[64] Frank Verstraete, J Ignacio Cirac, and José I Latorre. “Quantum circuits for strongly
correlated quantum systems”. Phys. Rev. A 79, 032316 (2009).
[65] Andrew J Ferris. “Fourier transform for fermionic systems and the spectral tensor
network”. Phys. Rev. Lett. 113, 010401 (2014).

Accepted in Quantum 2023-03-16, click title to verify. Published under CC-BY 4.0.

37

[66] Ryan Babbush, Nathan Wiebe, Jarrod McClean, James McClain, Hartmut Neven,and Garnet Kin-Lic Chan. “Low-Depth Quantum Simulation of Materials”. Phys.
Rev. X 8, 011044 (2018).
[67] R E F Silva, Igor V Blinov, Alexey N Rubtsov, O Smirnova, and M Ivanov. “Highharmonic spectroscopy of ultrafast many-body dynamics in strongly correlated systems”. Nat. Photonics 12, 266–270 (2018).
[68] Yuta Murakami, Martin Eckstein, and Philipp Werner. “High-Harmonic Generationin Mott Insulators”. Phys. Rev. Lett. 121, 057405 (2018).
[69] Yuta Murakami, Shintaro Takayoshi, Akihisa Koga, and Philipp Werner. “Highharmonic generation in one-dimensional Mott insulators”. Phys. Rev. B 103,
035110 (2021).
[70] In order to capture nontrivial phenomena in the steady states, we should attach acoupling with an external bath to the time-independent term H0 . While it alters theblock-encoding of H0 , the scaling of the computational cost does not severely increase.
[71] Koki Chinzei and Tatsuhiko N. Ikeda. “Time Crystals Protected by Floquet Dynamical Symmetry in Hubbard Models”. Phys. Rev. Lett. 125, 060601 (2020).
[72] Jérémie Roland and Nicolas J. Cerf. “Quantum search by local adiabatic evolution”.
Phys. Rev. A 65, 042308 (2002).
[73] Daniel A Lidar, Ali T Rezakhani, and Alioscia Hamma. “Adiabatic approximation
with exponential accuracy for many-body systems and quantum computation”. J.
Math. Phys. 50, 102106 (2009).
[74] Nathan Wiebe and Nathan S Babcock. “Improved error-scaling for adiabatic quantum
evolutions”. New J. Phys. 14, 013024 (2012).
[75] Kianna Wan and Isaac H Kim. “Fast digital methods for adiabatic state preparation” (2020). arXiv:2004.04164.
[76] Di Fang, Lin Lin, and Yu Tong. “Time-marching based quantum solvers for timedependent linear differential equations”. Quantum 7, 955 (2023).
[77] Heinz-Peter Breuer, Francesco Petruccione, and School of Pure and Applied Physics
Francesco Petruccione. “The Theory of Open Quantum Systems”. Oxford University
Press. (2002).
[78] Angel Rivas and Susana F Huelga. “Open Quantum Systems”. Springer Berlin
Heidelberg. (2012).
[79] Derek W Robinson. “Properties of propagation of quantum spin systems”. ANZIAM
J. 19, 387–399 (1976).
[80] Bruno Nachtergaele and Robert Sims. “Lieb-Robinson bounds and the exponential
clustering theorem”. Commun. Math. Phys. 265, 119–130 (2006).
[81] Bruno Nachtergaele, Yoshiko Ogata, and Robert Sims. “Propagation of correlationsin quantum lattice systems”. J. Stat. Phys. 124, 1 (2006).

Accepted in Quantum 2023-03-16, click title to verify. Published under CC-BY 4.0.

38

Appendix
A Short explanation of Eq. (7)Our algorithm is mainly based on the fact that the dynamics in the Floquet-Hilbert space,|ψ ∞ (t)i =

lim

lmax

e−ilωt hl| e−iHeff

X

lmax →∞

t(157)|0i |ψ(0)i ,l∈Dlmax

exactly corresponds to the solution |ψ(t)i [45]. Here, we roughly explain this relation forthis paper to be self-contained. We differentiate the above formula in t;iXd ∞|ψ (t)i =e−ilωt hl| (Heff + lω)e−iHeff t |0i |ψ(0)i
dt
l∈Z

=

XX

e−ilωt H−m hl + m| e−iHeff t |0i |ψ(0)il∈Z m∈Z



= 


X

H−m eimωt  

m∈Z


X

e−ilωt hl| e−iHeff t |0i |ψ(0)il∈Z= H(t) |ψ ∞ (t)i .

(158)

Noting that |ψ ∞ (0)i = |ψ(0)i, the uniqueness of the solution of the Schrödinger equationsuggests the relation |ψ ∞ (t)i = |ψ(t)i.

B Proof of the theorems for formulation
B.1 Amplitude amplification by symmetry
In Section 5.1, we discuss the amplification of the time-evolved state |ψ(t)i with exploiting
the translation symmetry of generic time-periodic Hamiltonians. With the usage of the
protocol,
P
4lmax
lmax
lmax
4lmax † −it l lω|lihl| −iHeff
t
e
Uini
,
(159)
Uamp1
(t) = (Uini)ewe see that it generates the time-evolved state |ψ(t)i with O(1) amplitude as1lmax
h0| Uamp1(t) |0i |ψ(0)i ' (|ψ(t)i + O(ε)).
2

(160)

We rigorously prove this statement, which is summarized as Theorem 6.
In order to show Theorem 6, we begin with discussing the approximate translation
symmetry. To be precise, we evaluate how much error appears in the approximation,
4lmax
thl| e−iHeff

4lmax

0

|l0 i ' eil ωt hl − l0 | e−iHeff

t

(161)

|0i ,

which we refer to as Eq. (69) in the main text. The exact upper bound on this error isgiven by the following theorem.
Lemma 8. (Approximate translation symmetry)
We choose the truncation order lmax ∈ Θ(γt + log(1/ε)/ log log(1/ε)) by Eq. (59). For
indices l0 ∈ Dlmax and l ∈ D4lmax , the inequality
0

hl|e

4lmax
−iHefft 0

il0 ωt|l i − e

hl

0

l |e

4lmax
−iHefft(γt)n(l ,l)|0i ≤ 8n(l0 , l)!

Accepted in Quantum 2023-03-16, click title to verify. Published under CC-BY 4.0.

(162)
39

is satisfied with

(8lmax − 2mmax − |l| − |l0 |).
n(l , l) =mmax

Here, l





0

(163)

l0 ∈ D4lmax denotes the difference l − l0 defined modulo 8lmax .

Proof.— We consider dynamics under a perturbed Hamiltonian Hpert (t), defined by
4lmax
Hpert (t) = Heff
+ Hb (t),Hb (t) =(164)(|li hl ⊕ m| ⊗ e

Xi8lmax ωt

H−m + h.c.),(165)(l,m)∈∂F 4lmax

with ∂F 4lmax = {(l, m) | 4lmax − mmax + 1 ≤ l ≤ 4lmax , 4lmax − l + 1 ≤ m ≤ mmax }. In thePinteraction picture based on H0 = l∈D4lmax |li h0|⊗(H0 −lω), the perturbed Hamiltonian
Hpert,I (t) is exactly translation symmetric under arbitrary shift |li → |l ⊕ mi as
Hpert,I (t) =

4lmax
I
Addm
⊗ e−imωt Hm
(t).

X

(166)

m6=0;|m|≤mmax

In the original frame, this exact symmetry implies the satisfaction of0hl|Upert (t)|l0 i = eil ωt hl

l0 |Upert (t)|0i ,(167)

where Upert (t) indicates the time evolution operator under the perturbed Hamiltonian
Hpert (t). Upon this relation, we evaluate the upper bound on the error by a triangle
inequality,[l.h.s. of Eq. (162)] ≤4lmaxhl|e−iHefft− Upert (t)|l0 i + hl

4lmax

l0 |e−iHefft− Upert (t)|0i .
(168)

We begin with computing the bound of the first term. Using the Dyson series expansion
in the interaction picture, similarly to Eqs. (37)-(42), it can be bounded by
4lmaxthl|e−iHeff
≤

∞ Z t
X
n=0 0

− Upert (t)|l0 i

dtn . . .

Z t2
0dt1 hl|nY

Hpert,I (ti ) −

i=1

n
Y

HI4lmax (ti ) |l0 i .

(169)

i=1

As discussed in the proof of Theorem 3, each integrand is decomposed into contributions from the transition amplitudes via paths |l0 i → |l1 i → . . . → |li. The difference
4lmax
between the Hamiltonians Hpert (t) and Heffarises only when the path goes across
4l
4l
−m
max
max
max
D
\D
, which is the support of Hb (t). For low order terms with n < n(l0 , l),defined by Eq. (163), such a nontrivial path is absent, and we obtain
hl|e

4lmax
−iHefft∞
X

tn− Upert (t)|l i ≤n!
n=n(l0 ,l)0

sup(kHpert,I (t)k)

n

t

+γ

n



.

(170)

Since the Hamiltonian Hpert,I (t) given by Eq. (166) is a circulant matrix (i.e. invariant
under the shift of |li) [51], its operator norm can be bounded bykHpert,I (t)k ≤ sup 
t0

= sup
t0


X

0
I
e−imωt Hm(t)eimωt 

m6=0

H(t − t0 ) − H0



= γ.

Accepted in Quantum 2023-03-16, click title to verify. Published under CC-BY 4.0.

(171)
40

As a result, Eq. (170) is further bounded by
4lmaxthl|e−iHeff− Upert (t)|l0 i ≤ 2

∞
X

(γt)n
.
n!
n=n(l0 ,l)(172)

In a similar way, we can obtain the bound on the second term of the right hand sideof Eq. (168), which results in the above formula under the replacement of n(l0 , l) byn(0, l l0 ). Since n(0, l l0 ) is always larger than n(l0 , l) from their definitions, we obtain
the accuracy of the approximate translation symmetry as[l.h.s of Eq. (162)] ≤ 4

∞
X

0

(γt)n
(γt)n(l,l )≤8.
n!
n(l0 , l)!
n=n(l0 ,l)(173)

When we choose the truncation order lmax by Eq. (59), the integer n(l0 , l) is always larger
than 6γt for indices l ∈ D4lmax and l0 ∈ Dlmax . We use Eq. (45) for the last inequality.
This lemma ensures that the approximate translation symmetry in the truncatedFloquet-Hilbert
space is extremely accurate; the right hand side of Eq. (162) is approx
imately O ε3 under the choice of lmax , Eq. (59). As we intuitively discuss in Section
5.1, the approximate translation symmetry ensures the amplification by symmetry as Eq.
(73). We next prove the consequence of Theorem 6, which gives the exact description of
Eq. (73);1εlmax
h0| Uamp1(t) |0i |ψ(0)i − |0i |ψ(t)i ≤ .
(174)
2
3
Proof of Theorem 6.— We track Eq. (70) in a rigorous way. We begin with thedefinition,lmax
h0| Uamp1(t) |0i |ψ(0)i =

1

X

X

4lmax 04lmax

e−ilωt hl| e−iHefft|l0 i |ψ(0)i .

(175)

l ∈Dlmax l∈D4lmax

We separate the summation over l ∈ D4lmax in the above formula by
X

=

+

X
l;l−l0 ∈D3lmaxl∈D4lmax

X

(176)

.

l∈D4lmax ;l−l0 ∈D/ 3lmax

Let us focus on the first summation. For each l0 ∈ Dlmax , the summation can be approximated as
4lmax

e−ilωt hl| e−iHeff

Xt|l0 i |ψ(0)i − |ψ(t)il;l−l0 ∈D3lmax
0(γt)n(l ,l)+8n(l0 , l)!
3lmax

X

≤

l;l−l0 ∈D

8X≤l∈D4lmax

4lmax

e−ilωt hl| e−iHeff

Xt|0i |ψ(0)i − |ψ(t)il∈D3lmax(γt)n(lmax ,l)+ 10mmaxn(lmax , l)!



emmax γt3lmax

3lmax /mmax

(177)

.

In the first inequality, we employ Lemma 8 with taking l − l0 = l l0 for l − l0 ∈ D3lmax into
consideration. The last inequality comes from Theorem 4. The first term in Eq. (177) isfurther bounded by
16∞
X
l=3lmax −2mmax(γt)dl/mmax e≤ 16mmax(dl/mmax e)!

eγtd3lmax /mmax − 2ed3lmax /mmax −2e

Accepted in Quantum 2023-03-16, click title to verify. Published under CC-BY 4.0.

,

(178)

41

where we use the relation Eq. (54) and the Stirling’s formula Eq. (48). This accomplishes
the evaluation of Eq. (177) with the usage of Eq. (60) as[Eq. (177)] ≤ 26mmax"≤ 26mmax≤eγt
3lmax /mmax − 3eγtlmax /mmax − 1

3lmax /mmax −3lmax /mmax −1 #313
ε3 .
500(mmax )2(179)

We next compute the second summation in Eq. (176), which is taken over l ∈ D4lmax
satisfying l − l0 ∈/ D3lmax . The Lieb-Robinson bound, dictated by Theorem 3, immediately
concludes its upper bound by
0

X

4lmax

e−ilωt hl| e−iHefft|l0 i |ψ(0)i

2

X

≤

(γt)d|l−l |/mmax e(d|l − l0 |/mmax e)!

l;l−l0 ∈D/ 3lmax∞
X (γt)dl/mmax el∈D4lmax ;l−l0 ∈D/ 3lmax≤ 4

l=3lmax≤ 4mmax≤(dl/mmax e)!

emmax γt
3lmax

ε3
.
250(mmax )2

3lmax /mmax(180)

Finally, we show that the amplification protocol relying on the symmetry reproduces
the time-evolved state |ψ(t)i. Combining the results of Eqs. (179) and (180), we arrive at1lmax
h0| Uamp1(t) |0i |ψ(0)i − |0i |ψ(t)i
2

≤

1

13ε3
ε3
+
500(mmax )2 250(mmax )2X

4lmax 0l ∈Dlmax≤3
ε3 ,200(mmax )2

which is smaller than ε/3 for mmax ∈ N and ε ∈ [0, 1].

!

(181)



B.2 Validity of refined effective Hamiltonianlmax

In Section 6.2, we raise an efficient block-encoding protocol for implementing e−iHeff t by
the qubitization technique. The strategy is to substitute the effective Hamiltonian under
4lmax
4lmax
the periodic boundary conditions, Heff,pbc, for the original one Heff
; The complexityof the oracles reduces from O(lmax ) to O(log lmax ). In order to justify this replacement of
4lmax
4lmax
Heff
with Heff,pbc, we provide Theorem 7, which says that its error is bounded fromabove bylmax
h0| Uamp1,pbc(t) |0i |ψ(0)i −1|ψ(t)i2

lmax
Uamp2,pbc(t) |0i |ψ(0)i − |0i |ψ(t)i

ε
,
3

(182)

≤ ε.

(183)

≤

lmax
The amplification protocols under the refined effective Hamiltonian, Uamp1,pbc(t) andlmax
Uamp2,pbc (t), are designated by Eqs. (102) and (103).

Accepted in Quantum 2023-03-16, click title to verify. Published under CC-BY 4.0.

42

lmax
(t), i.e. Eq. (183), is derived by Eq.
Proof of Theorem 7.— The validity of Uamp2,pbc(182) following the discussion in Section 5.2. It is sufficient to prove the relation, Eq.
(182). We first evaluate the upper bound of
4lmax
4lmax
(t) |0i |ψ(0)i − h0| Uamp1(t) |0i |ψ(0)ih0| Uamp1,pbc

1

≤

X

4lmax

4lmax

4lmax
t

hl|e−iHeff,pbc t − e−iHeffX|l0 i .

(184)l∈D4lmax l0 ∈Dlmax

Each transition amplitude in the summation can be computed in a similar way to the proof
4lmax
4lmax
+ Hb
= Heffof Lemma 8. We decompose the refined effective Hamiltonian by Heff,pbc
with the perturbation
Hb =(|li hl ⊕ m| ⊗ H−m + h.c.).

X

(185)

(l,m)∈∂F 4lmax

Following the discussion in the proof of Lemma 8, the perturbation Hb which locally actson the domain ∂F 4lmax hardly changes the transition amplitude from the original effective
4lmax
Hamiltonian Heff
. As a result, the deviation can be bounded by
4lmax

4lmax

hl|e−iHeff,pbc t − e−iHeff

t 0|l i ≤∞
Xtn
n!
n=n(l0 ,l)

sup(kHpbc,I (t)k)

n

t



+ γn ,(186)

where the Hamiltonian in the interaction picture is defined by

Hpbc,I (t) = [U0 (t)]† 


X

max
Add4l
⊗ Hm  U0 (t).
m

(187)

m6=0

We again use the property of circulant matrices, and obtain kHpbc,I (t)k ≤ γ, which leads
to
hl|e

4lmax
−iHeff,pbc
t

−e

4lmax
−iHeff
t 0

|l i

∞
X

γt
(γt)n
≤4
2
≤
0 , l)n!
n(ln=n(l0 ,l)n(l0 ,l)

.

(188)

As a result, we get the difference between the original and refined effective Hamiltonians
as
4lmax
4lmax
h0| Uamp1,pbc(t) − Uamp1(t) |0i |ψ(0)i

≤

1
2

≤ 4

X
l∈D4lmax∞X4(γt)n(lmax ,l)n(lmax , l)!

l=3lmax −2mmax

3γt
3lmax /mmax − 2
ε3.
250(mmax )2≤ 4≤(γt)dl/mmax e(dl/mmax e)!

3lmax /mmax −2

(189)

4lmax
As stated by Theorem 6, the original amplification protocol Uamp1(t) accurately generates|ψ(t)i as Eq. (181). Therefore, the one under the refined effective Hamiltonian is also

Accepted in Quantum 2023-03-16, click title to verify. Published under CC-BY 4.0.

43

justified as
lmax(t) |0i |ψ(0)i −h0| Uamp1,pbc

1|ψ(t)i
2

4lmax
4lmax
4lmax
(t) |0i |ψ(0)i −≤ h0| Uamp1,pbc(t) − Uamp1(t) |0i |ψ(0)i + h0| Uamp1

≤

19
ε3 .
1000(mmax )2

1|ψ(t)i2(190)This bound is actually smaller than ε/3, and we complete the proof of Theorem 7.

C Extension to exponentially-decaying Fourier components
In the main text, we focus on time-periodic Hamiltonians which have vanishing Fouriercomponents Hm = 0 for |m| > mmax . Here, we generalize our results to time-periodic
Hamiltonians H(t) with exponentially-decaying Fourier components as
kHm k ≤ he−|m|/ζ ,h, ζ > 0,(191)

for |m| > 0. The norm of H0 is arbitrary as long as it is bounded. The effective Hamillmax, is the same as Eq. (18), but it
tonian in the truncated Floquet-Hilbert space, Heff
is dense in the basis {|li}l∈Dlmax compared to the cases where Hm = 0 is satisfied for|m| > mmax .
We formulate the protocol in a similar manner to the main text, and show that thecomputational resources for the time-evolved state have nearly-optimal dependence bothin t and 1/ε. The difference mainly comes from the form of the Lieb-Robinson boundand the infinite series of {Hm } needed for designating H(t). The former one affects thetruncation order lmax for the Floquet-Hilbert space. The latter one yields the change
in block-encoding so that almost all the information about H(t) can be embedded withkeeping the efficiency.

C.1 Truncation order of Floquet-Hilbert space
We first determine the proper truncation order lmax for the Floquet-Hilbert space, as we
did in Section 4. To this aim, we begin with deriving the Lieb-Robinson bound on thetransition rate, corresponding to Theorem 3.
Theorem 9. (Bound on transition rate)
We assume kHm k ≤ he−|m|/ζ with certain positive constants h and ζ. Then, for l, l0
such that |l|, |l0 | ≤ lmax , the transition rate is bounded from above bylmax

hl| e−iHefft|l0 i ≤ exp −|l − l0 | − 2βζ 0 ht
+ 2/β .
ζ0


(192)

Here, β and ζ 0 are positive constants defined by
β = (1 − e−1/ζ )−1 ,ζ 0 = (1/ζ − 1 + e−1/ζ )−1 .

(193)

Proof.— We prove the theorem in a similar way to Theorem 3. We start from theDyson series expansion in the interaction picture [See Eq. (42)],
hl| e

lmax
−iHeff
t

0

|l i

≤

∞ Z t
X
n=0 0

≤

∞ Z t
X
n=0 0

dtn . . .
dtn . . .

Z t2
0

Z t2
0

dt1
dt1

n
XY

hli |HI (ti )|li−1 i{li } i=1n
XYkhli |HI (ti )|li−1 ik,(194){li } i=1

Accepted in Quantum 2023-03-16, click title to verify. Published under CC-BY 4.0.

44where we insert the identity li ∈Dlmax |li i hli | = I for n − 1 times. The summation {li } istaken over li ∈ Dlmax for i = 1, 2, . . . , n − 1 under fixed l0 = l0 and ln = l. We introduce a
new variable mi = li − li−1 instead of using {li }, and then the above integrand is bounded
by
P

Fn ≡

P

n
XY

khli |HI (ti )|li−1 ik{li } i=1X≤δl0 +m1 +...+mn ,ln{mi }ni=1 ∈ZkH−mi ki=1

δl0 +m1 +...+mn ,l hn e−

X

≤

n
Y

Pn
i=1

|mi |/ζ

(195)

.

n
{mi }ni=1 ∈Z

Since l0 + m1 + . . . + mn = l implies |m1 | + . . . + |mn | ≥ |l − l0 |, it can be further bounded
as follows,
X

Fn ≤

θ

n
X

n
{mi }ni=1 ∈Z≤ (2h)≡|mi | − |l − l0 | hn e−

Pni=1|mi |/ζ

i=1

∞
X

n

!

θ

m1 ,...,mn =0n(2h) Sn (|l − l0 |).

nX!

mi − |l − l | e−
0

Pn
i=1

mi /ζ

.

i=1

(196)

Here, θ(x) is a step function, defined by θ(x) = 1 for x ≥ 0 and θ(x) = 0 otherwise. The
summation Sn (M ) is defined by
Sn (M ) =∞Xθ

m1 ,...,mn =0nX!

mi − M e−

P
i

mi /ζ

,

(197)

i=1

−m/ζ = βe−M/ζ .
which satisfies S1 (M ) = ∞m=M e
Let us evaluate the upper bound on Sn (M ). We split the summation over mn into theone over mn ≤ M − 1 and the one over mn ≥ M , which results inP

Sn (M ) =M
−1
X

e−mn /ζ Sn−1 (M − mn ) + β n e−M/ζ .

(198)

mn =0

We use this relation recursively until S1 appears. After the single use of this equality, we
obtainSn (M ) ≤M
−1 M −mX
Xn −1e−(mn +mn−1 )/ζ Sn−2 (M − mn − mn−1 )mn =0 mn−1 =0+β n−1 e−M/ζM−1X

1 + β n e−M/ζ ,

(199)

mn =0

where we use e−mn /ζ ≤ 1 for the second term. By repeating this calculation n − 1 times,

Accepted in Quantum 2023-03-16, click title to verify. Published under CC-BY 4.0.

45

we arrive at

Sn (M ) ≤ 

M −1−Mn−2

M
−1
X

...

m1 =0


 e−Mn−1 /ζ S1 (M − Mn−1 )X
mn−1 =0+e−M/ζ

n−1X



β n−k+1 = e−M/ζ

M −1−Mk−2



X

1

...

m1 =0

k=1
n
X

M
−1
X



β n−k+1 

M
−1
X

...

m1 =0

k=1

mk−1 =0

M −1−Mk−2



X

 1,

(200)

mk−1 =0k−1where ki=1 mk is denoted by Mk . In the above formula, the summation of 1 over {mi }i=1
M
represents the number of lattice points included in a (k − 1)-dimensional pyramid Pk−1 =P{(x1 , . . . , xk−1 ) ∈ [0, M − 1]k | 0 ≤ i xi ≤ M − 1}, and hence it is bounded by the volume
M +2
of Pk−1. This leads to the relation,P

Sn (M ) ≤ β n e−M/ζnX{β −1 (M + 2)}k−1(k − 1)!

k=1≤ β n e−(1/ζ−1/β)M +2/β ,(201)

We define a positive constant ζ 0 by 1/ζ 0 = 1/ζ −1/β as Eq. (193). The above inequality
enables to evaluate the integrand Fn as
00

Fn ≤ (2h)n Sn (|l − l0 |) ≤ (2βh)n e−|l−l |/ζ +2/β .

(202)

Finally, going back to the inequality Eq. (194), we arrive at the bound on the transition
rate,
lmax
t

hl| e−iHeff|l0 i ≤

∞ n
X
t
n=0

n!

0

0

0

0

(2βh)n e−|l−l |/ζ +2/β ≤ e2βht−|l−l |/ζ +2/β .

(203)This completes the proof of the theorem. 
Theorem 9 says that the bound on the transition amplitude from |l0 i to |li exponentiallydecays in the distance |l − l0 | for Hamiltonians with kHm k . e−O(|m|) .  The decay is0relatively slow compared to the cases Hm = 0 (|m| > mmax ) showing O |l − l0 |−|l−l | .
In fact, this result is reminiscent of the Lieb-Robinson bound for short-ranged interacting
systems [79, 80, 81], which provides exponentially-decaying correlation functions underinteractions Uij ∼ e−O(|i−j|) . Based on Theorem 9, we assess how the exact time-evolved
state can be approximated by the truncated Floquet-Hilbert space, and determine theproper truncation order lmax . We summarize the result by the following theorem, which
is a counterpart of Theorem 4.
Theorem 10. (Floquet-Hilbert space truncation)
We assume kHm k ≤ Ce−|m|/ζ . The exact time-evolved state |ψ(t)i is approximated by
the truncated state |ψ lmax (t)i [See Eq. (19)] as
0|ψ(t)i − |ψ lmax (t)i ≤ 4ζ 0 e2βht−(lmax −1)/ζ +2/β .

(204)

Proof.— The derivation is similar to the one for Theorem 4. According to Eqs. (50)(52), we should evaluate the following two values,
ε1 =

X
0

l∈(Dlmax \Dlmax )

l0max
t

hl| eiHeff

|0i ,

ε2 =

X

l0max
t

hl| eiHeff

lmax− eiHefft|0i , (205)l∈Dlmax

Accepted in Quantum 2023-03-16, click title to verify. Published under CC-BY 4.0.

46

0where the left hand side of Eq. (204) is bounded by ε1 + ε2 under lmax
→ ∞. Theorem 9soon concludes the upper bound of ε1 ;
0

X

ε1 ≤

e2βht−|l|/ζ +2/β0l∈(Dlmax \Dlmax )Z ∞2βht−x/ζ 0 +2/β

≤ 2

dxe

lmax −10 2βht−(lmax −1)/ζ 0 +2/β

≤ 2ζ e(206).

The second error ε2 can be evaluated by the Dyson series expansion like Eq. (42); eachterm of ε2 in Eq. (205) is bounded from above by
∞ Z t
X
n=0 0

dtn . . .

Z t2
0

n
Y
0
lmax
(ti ) −
HIlmax (ti )|0i .
HI
dt1 hl|
i=1
i=1
n
Y

l0

(207)

As stated in the proof of Theorem 3, each of hl| ni=1 HI max (ti )|0i and hl| ni=1 HIlmax (ti )|0i
is decomposed into a product of transition amplitudes via a path |0i → |l1 i → . . . →0|ln−1 i → |li. Their difference appears only when the path goes across Dlmax \Dlmax . In
other words, denoting the summation over the set of {li } taking such nontrivial paths by
P0{li } , we reachQ[Eq. (207)] ≤

∞ Z t
X
n=0 0

dtn . . .

Z t2
0

dt1

n
X Y
0

Q

l0

hli |HI max (ti )|li−1 i .

(208)

i=1

{li }

When we define the hopping distance mi by mi = li − li−1 , it should satisfy ni=1 |mi | ≥(lmax − |l|) + (lmax − 0) for {li } such that a nontrivial path is organized. This results in
the relation described by
P

ε2 ≤
≤

l∈Dlmax n=0
∞
X X
l∈D

≤

∞ n
X
t

X

n
Y

X

n! {m }ni i=1 ∈Zn!

kH−mi k × θnX!

|mi | − 2lmax + |l|

i=1

i=1

(2ht)nSn (2lmax − |l|)n!
lmax n=0∞X(2βht)nX

n!

l∈Dlmax n=00 2βht−(lmax −1)/ζ 0 +2/β

≤ 2ζ e0e−(2lmax −|l|)/ζ +2/β
(209)

.

Combining this inequality with |ψ(t)i − |ψ lmax (t)i ≤ ε1 + ε2 and Eq. (206), we obtain
Eq. (204), which completes the proof. 
Theorem 10 determines the proper truncation order of the Floquet-Hilbert space forthe algorithm. When we aim at the desirable accuracy as |ψ(t)i − |ψ lmax (t)i ≤ ε, it issufficient to choose the truncation order lmax bylmax =2ζ 0+ 1 ∈ Θ (ht + log(1/ε)) . (210)2βζ ht + ζ log(1/ε) + ζ log 4ζ +
β
0

0

0

0



The result has similar scaling in t and 1/ε to the one for Hamiltonians with a finite numberof Fourier components, described by Eq. (59), other than the factor of log log(1/ε). While
the actual mmax is infinite for exponentially-decaying Fourier components as kHm k ≤he−|m|/ζ , the parameters h and ζ 0 play a role of γ and mmax respectively.
Accepted in Quantum 2023-03-16, click title to verify. Published under CC-BY 4.0.

47

C.2 Amplitude amplification
We next verify the validity of the amplitude amplification protocols, discussed in Section 5,
for time-periodic Hamiltonians with exponentially-decaying Fourier components kHm k ≤he−|m|/ζ .
The first amplification relies on the translation symmetry of the effective Hamiltonian,as discussed in 5.1. To show its extension, we begin with discussing the approximatetranslation symmetry in the truncated Floquet-Hilbert space, which is a counterpart ofLemma 8.
Lemma 11. (Approximate translation symmetry)
We assume kHm k ≤ he−|m|/ζ , and consider the truncated Floquet-Hilbert space C8lmax ⊗H. The transition rate has an approximate translation symmetry in that it satisfies
4lmaxthl| e−iHeff

4lmax

0

|l0 i − eil ωt hl

l0 | e−iHeff

t

0

0

|0i ≤ 2e2βht−(8lmax −|l|−|l |)/ζ +2/β .

(211)

Proof.— The proof is done in a similar manner to that of Lemma 11. We consider aperturbation H˜b (t) designated byH˜b (t) =|li hl ⊕ m| ⊗ e8ilmax ωt H−m + h.c.,

X

(212)

(l,m)∈∂ F̃ 4lmax

with ∂ F̃ 4lmax = {(l, m) | l ∈ D4max , 8lmax − l + 1 ≤ m ≤ 8lmax − 1}. This Hamiltonian
indicates hopping terms that go across the boundaries |4lmax i and |−4lmax + 1i. Let
4lmax
U˜pert (t) denote a time evolution operator under Heff
+ H˜b (t). Then, due to the exacttranslation symmetry in the interaction picture, the transition amplitude hl|U˜pert (t)|l0 isatisfies
0hl|U˜pert (t)|l0 i = eil ωt hl l0 |U˜pert (t)|0i .
(213)
4l
4lmax
+ H˜b (t) and H max is bounded
The difference of transition amplitudes between H
eff

eff

in a similar way to Eqs. (207) and (209). The difference survives only when the trajectory|l0 i → |l1 i → |ln−1 i → |li pass through the boundaries |4lmax i and |−4lmax + 1i viaPH˜b (t), and then its length ni=1 |mi | with mi = li − li−1 should be equal to or larger than
(4lmax − |l|) + (4lmax − |l0 |). We obtain its upper bound in a similar way to Eq. (209),hl|U˜pert (t) − e−iHeff

4lmax

t 0

|l i

≤

n
X
(2ht)nn!

Sn (8lmax − |l| − |l0 |)n=02βht−(8lmax −|l|−|l0 |)/ζ 0 +2/β

≤ e

.

(214)

Using this relation twice and the symmetry Eq. (213) leads to Eq. (211). 
With the usage of this approximate translation symmetry, we can organize the amplitude amplification by symmetry like Theorem 6. We prepare the truncated Floquet-Hilbert
space C8lmax ⊗H, and make the initial state uniform
in |liwith lmax ∈ Θ(ht+log(1/ε)). The
4lmax
state resulting from the time evolution exp −iHeff t outputs the target time-evolved
state |ψ(t)i with amplitude 1/2 as follows.
Theorem 12. (Amplification by symmetry)
We assume kHm k ≤ he−|m|/ζ , and choose the truncation order lmax ∈ O(ht + log(1/ε))
by Eq. (210). Let us summarize the amplification protocol relying on the symmetry by
lmax
Uamp1
(t), whose explicit formula is given by Eq. (72). Then, it generates the time-evolved
state |ψ(t)i with O(1) amplitude as
51|ψ(t)i ≤ε3 .
264(ζ 0 )2(215)

Accepted in Quantum 2023-03-16, click title to verify. Published under CC-BY 4.0.

48

lmax
h0|Uamp1(t)|0i |ψ(0)i −It is reasonable to√assume that the desirable error ε is sufficiently small, not greater thanthe constant 8ζ 0 / 15. Then, the right hand side of the above formula is bounded by ε/3,
which reproduces Theorem 6.
Proof.— The proof follows that of Theorem 6, described in Appendix B.1. The left
hand side of the above inequality can be bounded by two contributions determined by Eq.
(176). The first contribution, which corresponds to Eq. (177), is given by
4lmax

e−ilωt hl| e−iHeff

Xt|l0 i |ψ(0)i − |ψ(t)il;l−l0 ∈D3lmax

≤

X

2βht−

2e

7lmax −|l|+ β2ζ0+ |ψ 3lmax (t)i − |ψ(t)il;l−l0 ∈D3lmax
0≤ 8ζ 0 e2βht−(3lmax −1)/ζ +2/β ≤ε3.
8(ζ 0 )2

(216)

Here, we used Lemma 11 in the second line, Theorem 10 in the third line, and the choiceof lmax , Eq. (210), in the last line. The second contribution, which is a counterpart of Eq.
(180), is bounded by
4lmax

e−ilωt hl| e−iHeff

Xt|l0 i |ψ(0)i

X

≤

0

0

e2βht−|l−l |/ζ +2/βl;l−l0 ∈D/ 3lmaxl;l−l0 ∈D/ 3lmax0≤ 2ζ 0 e2βht−(3lmax −1)/ζ +2/β ≤ε3.
32(ζ 0 )2(217)

By combining these results like Eq. (181), we complete the proof of Theorem 12. 
The above theorem ensures the amplification by symmetry also for time-periodic
−1Hamiltonians with kHm k ≤ he−|m|/ζ ; it enhances the amplitude from O lmax
to 1/2
only with additional O(log lmax ) elementary gates. In order to bring the amplitude up to
1 − O(ε), we need the oblivious amplitude amplification. Its validity immediately follows
from Theorem 12, as we discussed in 5.2. In other words, when we apply the protocol
lmax
Uamp2
(t), defined by Eq. (79), the time-evolved state |ψ(t)i can be obtained with accuracy
1 − O(ε);
15
lmax
Uamp2
(t) |0i |ψ(0)i − |0i |ψ(t)i ≤
ε3 .
(218)
64(ζ 0 )2
The right hand side is bounded by ε for sufficiently small ε.

C.3 Qubitization technique for effective Hamiltonian
lmax
Simulating |ψ(t)i via the amplitude amplification Uamp2(t) requires the time evolution









4lmax
4lmax
operators exp −iHLPt and exp −iHeff
t . The former one is the same as Section

6.1. We hereby present how the latter one is implemented for the cases kHm k ≤ he−|m|/ζ .
We take a similar strategy to Section 6.2, that is, we first organize a refined effective
4lmax
4lmax
Hamiltonian Heff,pbc, which can accurately reproduce the dynamics under Heff
. After
that, we compose its block-encoding which can be efficiently achievable.

Accepted in Quantum 2023-03-16, click title to verify. Published under CC-BY 4.0.

49

C.3.1

Refined effective Hamiltonian

4lmax, which acts on the truncated FloquetWe introduce a refined effective Hamiltonian Heff,pbc
Hilbert space C8lmax ⊗ H, by
4lmax
4lmax
Heff,pbc= Heff+ H˜b ,H˜b =X|li hl ⊕ m| ⊗ H−m + h.c.,(219)(l,m)∈∂ F̃ 4lmax

where ∂ F̃ 4lmax shares the definition with the one in Eq. (212). Owing to the additionalterm Hb , the hopping terms induced by H−m become translation symmetric in |li as
4lmax
Heff,pbc
=

X

4lmax
max
Add4l
⊗ Hm − HLP
.
m

(220)

m∈D4lmaxmaxHere, a full quantum adder Add4l, defined by Eq. (108), appears and it allows efficientmimplementation of block-encoding as discussed later. Before going to its block-encoding,we prove the validity of the refined effective Hamiltonian as a counterpart of Theorem 7.

Theorem 13. (Refined effective Hamiltonian)lmax
We assume kHm k ≤ he−|m|/ζ , and organize the two amplification protocols Uamp1,pbc(t)lmax
and Uamp2,pbc (t) respectively based on Eqs. (102) and (103) with using the refined effective
4lmax
given by Eq. (219). When the truncation order lmax ∈ Θ(ht +Hamiltonian Heff,pbc
log(1/ε)) is chosen by Eq. (210), they also provide the exact time-evolved state |ψ(t)i as
1|ψ(t)i2≤

lmax
Uamp2,pbc(t) |0i |ψ(0)i − |0i |ψ(t)i≤

lmax(t) |0i |ψ(0)i −h0| Uamp1,pbc

11
ε3 ,128(ζ 0 )233
ε3 ,128(ζ 0 )2

(221)
(222)

for arbitrary initial states |ψ(0)i ∈ H. For a allowable error ε sufficiently small comparedto the constant ζ 0 , both of the left hand sides are smaller than ε.
Proof.— We prove the theorem as we do for Theorem 7 in Appendix B.2. First, we
4lmax
4lmax
evaluate the difference of the transition rates between Heff,pbc
and Heff. We replace
the perturbation H̃b (t) by the boundary term H̃b in the proof of Lemma 11 (See AppendixB.1). We obtain the same result as Eq. (214),
4lmax

4lmax

hl| e−iHeff,pbc t |l0 i − hl| e−iHeff

t

0

0

|l0 i ≤ e2βht−(8lmax −|l|−|l |)/ζ +2/β .

(223)

Once we obtain this bound, we can track the proof of Theorem 7, composed of Eqs.
lmax(184)-(190). First, we evaluate the deviation from Uamp1(t) as
0

4lmax
4lmax
h0| Uamp1,pbc(t) − Uamp1(t) |0i |ψ(0)i≤X
l∈D4lmax ,l0 ∈D≤0e2βht−(8lmax −|l|−|l |)/ζ +2/β
4lmax
lmax

1 0 2βht−(3lmax −1)/ζ 0 +2/βε3ζe≤
. (224)2
128(ζ 0 )2

4lmax
Since h0| Uamp1(t) |0i |ψ(0)i accurately provides the state |ψ(t)i /2 as Eq. (215) according
to Theorem 12, a triangle inequality concludes Eq. (221). Since the oblivious amplitude
amplification generates the error at-most three times larger than Eq. (221) according to
the discussion in Section 5.2, Eq. (222) is immediately derived. 

Accepted in Quantum 2023-03-16, click title to verify. Published under CC-BY 4.0.

50

C.3.2Block-encoding

4lmax
to implement its
We compose block-encoding of the refined effective Hamiltonian Heff,pbc
time evolution operator by qubitization.
We begin with describing the assumption on time-periodic Hamiltonians H(t). If we
adopt the one for the main text as Eq. (98), we should have the knowledge about blockencoding for the infinite series {Hm }m∈Z . To avoid this complexity, we instead consider
time-periodic Hamiltonians H(t) written by

H(t) =jXmax

αj (t)Mj ,αj (t + T ) = αj (t).

(225)

j=1
max
max
with coefficients {αj (t)}jj=1
and operators {Mj }jj=1. We assume the knowledge of blockencoding for each operator Mj as

hGj |Oj |Gj i = Mj ,|Gj i = Gj |0i⊗na ,

na

(226)
na

with an oracle state |Gj i ∈ C2 and an oracle unitary gate Oj on C2 ⊗ H. We canalways set the denominator in block-encoding by 1 since its absolute value and its phasecan be absorbed respectively into the coefficient αj (t) and the oracle Oj . The Fouriercomponents are given by
Hm =

jX
max

αjm =

αjm Mj ,

j=1

1
T

Z T
0

dtαj (t)eimωt ,(227)

and it is always possible to impose αjm ≥ 0 for every j, m. This is because, when a certain
coefficient αjm is complex, we can divide it into αjm = Re αjm +i Im αjm and redefine Mj withthe signs of Re αjm , Im αjm and the phase i. In other words, adding −Mj , iMj , and −iMj
to the set {Mj }j reproduces the form of Eq. (227) with αjm ≥ 0, where jmax is replaced by
at-most 4lmax . In the following discussion, the oracles Oj and Gj are respectively supposed
to be composed of at-most C elementary gates. In addition, we also assume the access tomax
the coefficients αjm for m ∈ Dlmax by the oracle unitary gate Glcoef
asmax
Glcoef|0i |0i =

X

jX
max

m∈Dlmax j=1

s

αjm|mi |ji ,αlmax

αlmax =

X

jX
max

αjm .

(228)

m∈Dlmax j=1

max
The oracle state Glcoef|0i |0i ∈ C2lmax +jmax is a O(log lmax + log jmax )-qubit ancillary state,max
and its preparation Glcoef
is supposed to require C lmax elementary gates. The parameterαlmax is bounded by the whole-system energy scale asαlmax ≤max
X jX

αjm ≡ α.

(229)m∈Z j=1

Although the definition of α is different from that of the main text, it again provides thebound on the scale of time-dependent terms as
γ≤

max
X jX

|αjm | · kMj k ≤ α,(230)m∈Z j=1

where we use kMj k ≤ 1 from Eq. (226).
Accepted in Quantum 2023-03-16, click title to verify. Published under CC-BY 4.0.

51

We move to how we organize an oracle state and oracle unitary gate for the refined
4lmax
Heff,pbc. Substituting the Fourier components Hm into Eq. (220) results in

4lmax
Heff,pbc
=

jX
max

X

4lmax
4lmax
αjm Addm⊗ Mj − HLP
.

(231)

m∈D4lmax j=1

We introduce five kinds of ancillary systems labeled by a, b, c, d, and e. The systems a, b,
and e are respectively described by 2na -, 8lmax -, and 2-dimensional Hilbert spaces, as we
consider in Section 6.2. States of the system c is spanned by {|mi}m∈D4lmax due to themaxabsence of mmax . The system d is characterized by {|ji}jj=1. The oracle unitary gate for4lmax
the refined effective Hamiltonian Heff,pbc is given by
4lmax
Oeff

OJ

4lmax
= |0i h0|e ⊗ OJ ⊗ Add4lmax ⊗ Ib − |1i h1|e ⊗ Id ⊗ Ic ⊗ OLP⊗ Ia ,=

jX
max

Add4lmax =|ji hj|d ⊗ Oj ,

j=1

X

max
|mi hm|c ⊗ Add4l
.
m

(232)
(233)

m∈D4lmax

4lmax
max
are respectively implemented by O(jmax C),
Here, the unitary gates OJ , Add4l
and OLPm
O(log lmax ), and O(log lmax ) elementary gates (the second one represents addition of a varimax
i on the auxiliary systems
able m ∈ D4lmax [53]). Next, we provide the oracle state |G4l
eff
a, b, c, d, and e by
4lmax
max
max
max
ib |0ia ,|w4lmax ie (G4l
i = G4l|G4l
eff
coef |0id |0ic ) |a
eff
max
G4l
eff

= Ie ⊗
√

|w

4lmax

ie =

jX
max

|ji hj|d ⊗ Ic ⊗ Ib ⊗ (Gj )a ,j=1

α4lmax |0ie +
p

(234)
(235)

√

4lmax ω |1ie.
4l
α max + 4lmax ω(236)

We can confirm the relation,
4lmax
max
max
hG4l
|Oeff
|G4l
i=
eff
eff

4lmax
Heff,pbc

α4lmax + 4lmax ω

(237)

,

by directly substituting the above equations. The cost of implementing the oraclestate
4lmax|Geff i from the trivial state |0ie |0id |0ic |0ib |0ia is O jmax C + C 4lmax + log lmax elementary gates.
Therefore, implementing the time evolution operator exp(−iHeff,pbc t) with the qubitization technique requires the following resources for time-periodic Hamiltonians satisfying
kHm k ≤ he−|m|/ζ ;• Number of ancillary qubits; O(log jmax + log lmax )• Number of overall gates;O {(α + lmax ω)t + o(log(1/ε))} jmax C + C 4lmax + na + log jmax + log lmax.

We replace α4lmax by α based on the inequality, Eq. (229).

Accepted in Quantum 2023-03-16, click title to verify. Published under CC-BY 4.0.

52

C.4 Algorithm and comparison
The algorithm for computing the time-evolved state for time-periodic Hamiltonian withkHm k ≤ he−|m|/ζ is almost the same as those in Section 7. We separately consider
adiabatic-like cases and generic long-time cases. The computational resources are obtainedby replacing γ in Table 1 by h. The query complexity is determined by the coefficients ofjmax C and C 4lmax in the number of overall elementary gates. Their scaling is summarized
as follows;
Adiabatic-like cases.—
The scaling of the query complexity is optimal in the time t and nearly-optimal in the
inverse error 1/ε.
• Number of ancillary qubits; na + O(log jmax + log(γt) + log log(1/ε)).
• Query complexity; O(αt + log(1/ε)).
• Additional gates per query; O(na + log jmax + log(γt) + log log(1/ε)).
Generic long-time cases.—
The scaling of the query complexity is nearly-optimal in the time t and the inverse error
1/ε. For practical problems up to poly (N )-time, it becomes optimal in the time t.
• Number of ancillary qubits; na + O(log jmax + log(γ/ω) + log log(ωt/ε)).
• Query complexity; O(αt + ωt log(ωt/ε)).
• Additional gates per query; O(na + log jmax + log(γ/ω) + log log(ωt/ε)).
We note that log log(1/ε) factor in the query complexity is buried by lmax . Thisdifference from the cases Hm = 0 (|m| > mmax ) essentially originates from the form of
the Lieb-Robinson bound dictated by 9. A slightly slow decay in the transition amplitude
results in slightly larger cost for the cases kHm k ≤ he−|m|/ζ , compared to the results in themain text. This leads to the nearly-optimal dependence of the query complexity in 1/ε.
Anyway, since the factor log log(1/ε) is not so large compared to others, the relation to
the quantum algorithms for time-independent and generic time-dependent Hamiltonians,
discussed in Section 7, is maintained. Importantly, in the query complexity, the wholeenergy scales αt and γt are separated from the O(log(1/ε)) term; the query complexity
is usually much smaller than that of the truncated Dyson-series algorithm, and close
to that of the qubitization technique. The number of ancillary qubits, which scales as
O(log log(1/ε)) for an allowable error ε, also lies between those for these algorithms.
We also mention about the complexity of the oracles for the cases kHm k ≤ he−|m|/ζ . In
contrast to the cases Hm = 0 (|m| > mmax ), we require a quantum circuit which embeds
the function αjm to the amplitude of a quantum state with O(jmax (γt + log(1/ε))) ormax
O(jmax (γ/ω + log(ωt/ε))) degrees of freedom, by the oracle G4l
coef defined by Eq. (228).
While integrable functions are efficiently implemented [54], the number of elementarygates for the worst cases amounts to the dimension of the quantum state, i.e. C 4lmax ∈O(jmax (γt + log(1/ε))) or O(jmax (γ/ω + log(ωt/ε))). While the preparation of oracles
is relatively difficult compared to the cases Hm = 0 (|m| > mmax ), we emphasize that
it is much easier than that of the truncated Dyson-series algorithm. In the truncatedDyson-series algorithm, we need more or less an oracle that embeds the information of
Accepted in Quantum 2023-03-16, click title to verify. Published under CC-BY 4.0.

53

time-dependent Hamiltonians at every discretized time into a quantum state. Such an
oracle is exemplified by Eq. (135) for generic Hamiltonians [19]. For time-dependent LCU
Hamiltonians or sparse-access Hamiltonians, it requires a quantum state holding all thecoefficients αj (t) at every discretized time t as an oracle [20]. The number of elementarygates per oracle can be proportional to that of the discretized time, O(jmax ωγt/αε), at
the worst case. Therefore, our algorithm also improves the complexity of the oracles aswell as the query complexity, compared to the truncated Dyson-series algorithm.

C.5 Example
For the cases kHm k ≤ he−|m|/ζ , we adopt the assumption of block-encoding as Eq. (227),in contrast to the cases Hm = 0 (|m| > mmax ). To show that it is reasonable and versatilein condensed matter physics and quantum chemistry, we provide a simple example.
We consider an N -site Fermi-Hubbard model HHub , given by Eq. (141). While weconsider an idealized laser light which has a constant amplitude in Section 8.1, we focus
on a rather realistic case where Gaussian wave packets of lights are shone to materials[69]. Then, if the time scale of the wave packets is sufficiently larger than the light withthe frequency Ω, the external drive Hext (t) is well described by

Hext (t) = 

X

2

2

e−(t−(n+1/2)T ) /2τ sin ΩtXVx n̂xσ ,

(238)

x,σ

n∈Z

where peaks of the wave packets are located at t = (n + 1/2)T . The period T is designatedby T = p(2π/Ω) with the number of the waves of the light for each wave packet p ∈ N.
Let us organize the block-encoding required for simulation. By carefully computing
Fourier components of H(t) = HHub + Hext (t), we arrive at H0 = HHub and
n
o
X
1 2
ωτ
2
2
Hm = i(−1)p+m sgn(m)AmVx n̂xσ , Am = √ e− 2 (p +m )(ωτ ) sinh p|m|(ωτ )2 ,
2π
x,σ
(239)
for m 6= 0. We define sgn(m) by sgn(m) = m/|m| for m 6= 0 and sgn(0) = 0. Since Am ≥2
0 rapidly decays as e−O(m ) , this Hamiltonian is suitable for our algorithm, satisfying
kHm k ≤ he−|m|/ζ . This satisfies the assumption of Eq. (227) with non-negative coefficients
αjm by choosing jmax = 3;
M1 =

HHub
P
,
2 x x + U N

M2 = +iM3 = −iPx,σ Vx n̂xσ

2

P

x Vx

,

2

P

x Vx

α1m =

2

x + U N δm0 ,

(240)

x

α2m = 2X

Vx Am δ(−1)p+m sgn(m),+1 ,(241)

Vx Am δ(−1)p+m sgn(m),−1 .

(242)

x

P

x,σ Vx n̂xσ!
X,α3m = 2Xx

The block-encoding for {Mj } can be composed by Eqs. (143)-(146) in Section 8.1. Therefore, by resorting to our algorithm for the adiabatic-like cases, we can efficiently simulate
the response to a realistic wave packet of laser light with the number of ancillary qubits
and the query complexity provided in Appendix C.4.

D Efficient algorithm for high-frequency regimes with ω ' α, γ
In Section
7.2, we concentrate on the generic long-time regime with the frequency ω ∈
0
O N for the number of sites, while the whole energy scales α, γ are given by poly (N ).
Accepted in Quantum 2023-03-16, click title to verify. Published under CC-BY 4.0.

54

Here, we show its reasonableness by showing that quantum simulation of broad timeperiodic Hamiltonians under the frequency comparable to the whole energy scale can bealmost attributed to that of time-independent Hamiltonians.
As a generic setup for time-periodic Hamiltonians, we assume the locality [= thesystem involves at-most O(1)-body interactions] and the extensiveness. The latter one
means that the energy scale per site, represented by λ, is bounded. In our definitions, thescale λ is approximated by O(α/N, γ/N ). When the ratio of the frequency ω to the local
energy scale λ is sufficiently large as ω/λ  1, we can apply the perturbative expansionin (ω/λ)−1 , called the high-frequency expansion [59, 60, 61, 62, 63]. Recent developments
in generic local and extensive Floquet systems have revealed the accurate upper bound on
the error as follows [60];l

U (nT ) − e−iHFM nT ≤ C1 λN nT e−C2 l0 + C3 nN (λT )l+2 (l + 1)!.

(243)

In the above inequality, C1 , C2 , and C3 represent some positive constants, and the integerl0 ∈ O(ω/λ) provides the optimal truncation order of the perturbation theory, while thelactual truncation order l (≤ l0 ) is optional. The time-independent Hamiltonian HFM
isthe l-th order Floquet-Magnus expansion, given by
l
HFM
=

l
X
l0 =0

(l0 )

HFM ,

(0)

HFM = H0 ,

(1)

HFM =

1
2iT

Z T
0

dt1

Z t1
0

dt2 [H(t1 ), H(t2 )],

(244)

(l0 )

where each l0 -th order term in ω/λ, represented by HFM , involves l0 -fold multi-commutatorsof H(t) (See Ref. [60] for the explicit formula).
Based on the above formalism, we see that the long-time dynamics at t = (n+δ)T  Twith n ∈ N and δ ∈ [0, 1) can be computed in a trivial way when ω is comparable toα, γ ∈ poly (N ). In this regime, the integer l0 is O(N ) due to the local energy scale λ ∈O(α/N, γ/N ). Then, we simulate the dynamics U (t) = U (δT )U (nT ) with approximatingthe latter one U (nT ) by Eq. (243). We choose the truncation order l by O N 0 . SincelH(t) generally
has poly (N ) local terms, the Floquet-Magnus expansion HFM
for such0ll ∈ O N can be identified with poly (N )-time classical computation, and HFM
is alsocomposed of poly (N ) localterms. Therefore,
we can apply the qubitization technique tolthe implementation of exp −iHFM nT with the cost designated by Eq. (13). Since theright hand side of Eq. (243) is bounded by O nN −l−1 , the replacement of U (nT ) byl nT
exp −iHFM

with an allowable error ε is valid as long as
n ≤ Const. × εN l+1 .

(245)

We usually employ a polynomial accuracy ε ∈ O(N −ν1 ) in quantum computation. The assumption for the frequency, represented by ω ∈ O(N ν2 ), implies the period T ∈ O(N −ν2 ).
Therefore, by properly choosing the O N 0 truncation order l ≥ ν1 + ν2 −1, we canefficiently simulate U (nT ) for arbitrary polynomial time nT ∈ O N l+1−ν1 −ν2 with thequbitization technique for time-independent Hamiltonians. The remaining micromotion
U (δT ) can be implemented with our protocol for the adiabatic-like regime (See Section7.1), but it does not affects the computational cost due to δT  nT .
In short, the dynamics under time-periodic Hamiltonians H(t) with ω ∼ α, γ can be
simulated with the time-independent Hamiltonian approaches. As a result, it is reasonable
to concentrate on the cases where ω is negligible compared to α and γ, as we do in Section
7.2.
Accepted in Quantum 2023-03-16, click title to verify. Published under CC-BY 4.0.

55