Optimal Hamiltonian simulation for time-periodic systems

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

Keisuke Fujii
Graduate School of Engineering Science, Osaka University, 1-3 Machikaneyama, Toyonaka, Osaka 560-8531, Japan.
Center for Quantum Information and Quantum Biology, Osaka University, Japan.
RIKEN Center for Quantum Computing (RQC), Hirosawa 2-1, Wako, Saitama 351-0198, Japan
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 $\varepsilon$.
In contrast, those for time-dependent systems require larger cost due to the difficulty of handling time-dependency.
In this paper, we establish optimal/nearly-optimal Hamiltonian simulation for generic time-dependent systems with time-periodicity, 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 query complexity, which measures the cost for implementing the time-evolution, has optimal and nearly-optimal dependency respectively in time $t$ and inverse error $\varepsilon$, and becomes sufficiently close to that of qubitization.
Thus, our protocol tells 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 and quantum chemistry, and quantum tasks yielding time-dependency in quantum computation.

§ 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 <cit.>.
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 application of Hamiltonian simulation nowadays ranges from condensed matter physics to quantum chemistry, such as quantum dynamics <cit.> and quantum phase estimation <cit.>.
For long time, Trotterization has been a standard way of Hamiltonian simulation for both time-independent and time-dependent systems, which can provide simple realization available in today's quantum computers <cit.>.
Instead, it requires a huge number of elementary gates up to $\poly{1/\varepsilon}$ to achieve an acceptable error $\varepsilon$.

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}$ <cit.>.
Among them, the qubitization technique <cit.> achieves the best cost in that its query complexity (a measure of cost) has optimal scaling both in time $t$ and inverse error $1/\varepsilon$ in an additive way. <cit.>.
In contrast, for time-dependent systems, while there exist several efficient Hamiltonian simulation for constructing time-evolution $U(t)$ <cit.>, they largely rely on the truncated Dyson-series expansion.
Although the query complexity of the truncated Dyson-series algorithm has optimal/nearly-optimal dependence in $t$ and $1/\varepsilon$ <cit.>, it scales in a productive way 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 a period.
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 phases of matter absent in time-independent systems, such as topological phases <cit.> and time crystals <cit.>.
Time-periodic Hamiltonians also cover optical responses in solids and molecules <cit.>, exemplified by high-harmonic generation and photo-chemical reactions.
In addition, adiabatic quantum dynamics such as Thouless pumping <cit.> and adiabatic state preparation for quantum computation <cit.> can be regarded as a part of time-periodic Hamiltonian dynamics under sufficiently large period $T$.
Despite fundamental significance and potential application of time-periodic Hamiltonians in broad fields, there is no Hamiltonian simulation protocol
which can efficiently handle their time-dependency, while their simulation itself is possible by the truncated Dyson-series algorithm.

In this paper, we develop an optimal and/or nearly-optimal quantum algorithm for time-periodic Hamiltonian simulation.
The key idea of our formalism is mapping time-dependent systems to effective time-independent systems in the infinite-dimensional Floquet-Hilbert space, obtained by the Fourier series expansion in time.
Although the infinite-dimensionality prohibits its simulation on quantum computers, we resolve this by formulating the Lieb-Robinson bound <cit.>, amplitude amplification, and the qubitization technique.
Consequently, we can simulate the target time-evolved state with arbitrary small 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/\varepsilon$, 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 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 matter physics and quantum chemistry, and the optimal control in quantum computation.

The rest of this paper is organized as follows.
In Section <ref>, we briefly review Floquet theory and the qubitization for this paper to be self-contained.
We provide our main results in Sections <ref>-<ref>, with firstly summarizing them in Section <ref>.
In Section <ref>, 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 <ref> and <ref> are devoted to deriving the subroutines of the algorithm. They respectively provide the amplification protocol that can realize the time-evolved state with sufficiently high success probability (Section <ref>) and efficient implementation of the time-evolution operator in the Floquet-Hilbert space (Section <ref>).
Section <ref> 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 <ref>, and concludes our paper in Section <ref>.

§ PRELIMINARIES

§.§ Review on Floquet theory

We briefly review Floquet theory to analyze Schrödinger equation under a time-periodic Hamiltonian,

\begin{equation}
i \dv{t} \ket{\psi(t)} = H(t) \ket{\psi(t)}, \quad H(t+T) = H(t),
\end{equation}

and its solution

\begin{equation}
\ket{\psi(t)} = U(t) \ket{\psi(0)}, \quad U(t) = \mcl{T} \exp \left( -i \int_0^t \dd{t^\prime} H(t^\prime) \right).
\end{equation}

Here, the time-evolved state of the system of interest is represented by $\ket{\psi(t)}$, and we assume that it is defined on a finite-dimensional Hilbert space $\mcl{H}$.
The set of states $\{ \ket{\psi_j} \}_{j=1}^{\mr{dim}(\mcl{H})}$ denotes a certain choice of the basis of $\mcl{H}$.
Floquet theorem says that the solution $\ket{\psi(t)}$ is always written in the form of

\begin{equation}
\ket{\psi(t)} = \sum_{\alpha=1}^{\mr{dim}(\mcl{H})} c_\alpha e^{-i \varepsilon_\alpha t} \ket{\phi_\alpha (t)}, \quad \ket{\phi_\alpha (t+T)}=\ket{\phi_\alpha (t)},
\end{equation}

similar to Bloch theorem for spatially-periodic systems. Here, $\varepsilon_\alpha \in [-\pi/T, \pi/T )$ and $\ket{\phi_\alpha(t)} \in \mcl{H}$ are respectively called quasienergy and Floquet state.
Time-periodicity of $H(t)$ and $\ket{\phi_\alpha(t)}$ allows the Fourier series expansion as

\begin{equation}
H(t) = \sum_{m \in \bbZ} H_m e^{-im\omega t}, \quad \ket{\phi_\alpha(t)} = \sum_{m \in \bbZ} \ket{\phi_\alpha^m} e^{-im\omega t},
\end{equation}

with the frequency $\omega = 2\pi /T$. The hermiticity of $H(t)$ implies $H_{-m}=H_m^\dagger$. We assume $\norm{H(t)} < \infty$ ($\norm{\cdot}$ ; operator norm), which results in $\sum_{m \in \bbZ} \norm{H_m}^2 < \infty$. Introducing an auxiliary degree of freedom $\{ \ket{l} \}_{l \in \bbZ}$ (sometimes called photon number) to relate $\ket{\phi_\alpha^l} \leftrightarrow \ket{l}\ket{\phi_\alpha^l}$, every pair of $\varepsilon_\alpha, \ket{\phi_\alpha(t)}$ can be obtained by the time-independent eigenvalue equation,

\begin{equation}
\ms{H}_\eff \left( \sum_{l \in \bbZ} \ket{l} \ket{\phi_\alpha^l} \right) = \varepsilon_\alpha \left( \sum_{l \in \bbZ} \ket{l} \ket{\phi_\alpha^l} \right),
\end{equation}

with the effective Hamiltonian defined by

\begin{equation}\label{Eq2A:def_effective_Hamiltonian}
\ms{H}_\eff = \sum_{l \in \bbZ} \ket{l}\bra{l} \otimes (H_0-l\omega) + \sum_{l,m \in \bbZ} \ket{l}\bra{l+m} \otimes H_{-m}.
\end{equation}

Intuitively, this time-independent Hamiltonian $\ms{H}_\eff$ describes a static one-dimensional system, where each $l$-th site has potential energy $H_0-l\omega$, as Fig. <ref>.
Off-diagonal terms $H_{-m}$ represent 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 vanish by this mapping; it is translated into  infinite dimensionality of the space spanned by $\{ \ket{l} \ket{\psi_j} \, | \, l \in \bbZ, \, j=1,2,\hdots, \mr{dim}(\mcl{H}) \}$, which is called Floquet-Hilbert space or Sambe space <cit.>.
While the effective Hamiltonian $\ms{H}_\eff$ is often used for identifying $\varepsilon_\alpha, \ket{\phi_\alpha(t)}$, it can be employed also for directly determining the dynamics as follows <cit.>;

\begin{equation}\label{Eq2A:psi_t_infinite}
\ket{\psi(t)} = \sum_{l \in \bbZ} e^{-il\omega t} \braket{l |e^{-i \ms{H}_\eff t}| 0}\ket{\psi(0)}.
\end{equation}

We can extract the dynamics under a time-dependent Hamiltonian $H(t)$ from that under the time-independent one $\ms{H}_\mr{eff}$ by preparing auxiliary systems labeled by $\ket{l}$.
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 infinite-dimensional.

Intuitive understanding of the effective Hamiltonian $\ms{H}_\eff$ [Eq. (<ref>)] in the Floquet-Hilbert space. It can be seen as a single-particle quantum system with the potential $H_0-l\omega$ and the hopping $H_{-m}$ on a infinite one-dimensional chain. We should introduce the cutoffs $\ket{-l_\imax + 1}$ and  $\ket{l_\imax}$ to embed it on quantum computers, which leads to $\ms{H}_\eff^{l_\imax}$ [Eq. (<ref>)].

§.§ 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 $\mcl{H}$ <cit.>.
It begins with assuming the existence of block-encoding, that is, we suppose an $n_a$-qubit auxiliary state $\ket{G}_a$ (called an oracle state) and a unitary gate $O$ on the space $\bbC^{2^{n_a}} \otimes \mcl{H}$ (called an oracle gate) satisfying

\begin{equation}\label{Eq2B:block_encode}
\braket{G|O|G}_a = \frac{H}{\alpha}, \quad \alpha \geq \norm{H}.
\end{equation}

The inequality $\alpha \geq \norm{H}$ comes from the norm of the unitary gate as $\norm{O}=1$, and the parameter $\alpha$ represents a typical energy scale of the whole system.
Here, the oracles $O$ and $\ket{G}_a$ are supposed to be efficiently implemented; the unitary gates $O$ and $G$, which realizes the oracle state from a reference state $\ket{0}_a$ as $G\ket{0}_a = \ket{G}_a$, require at most $C \in \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 of unitary (LCU) as

\begin{equation}\label{Eq3B:def_LCU}
H = \sum_{j=1}^{j_\imax} \alpha_j U_j, \quad \alpha_j \geq 0, \quad \text{$U_j$; unitary},
\end{equation}

a possible choice of the oracles is

\begin{equation}\label{Eq3B:block_encode_LCU}
O = \sum_{j=1}^{j_\imax} \ket{j}\bra{j}_a \otimes U_j, \quad \ket{G}_a = \sum_{j=1}^{j_\imax} \sqrt{\frac{\alpha_j}{\alpha}} \ket{j}_a,
\end{equation}

with $\alpha = \sum_{j=1}^{j_\imax} \alpha_j$.
The number of ancillary qubits $n_a$ scales as $\log j_\imax$.
Since $j_\imax$ is $\poly{N}$ for typical $N$-site systems, the number of elementary gates for the oracles amounts to $C \in \poly{N}$.
LCUs cover various types of Hamiltonians exemplified by quantum spin systems composed of local Pauli operators and fermionic systems in condensed matter physics and quantum chemistry <cit.>.
The block-encoding of other Hamiltonians, such as sparse-access matrices and purified density matrices, has been also revealed <cit.>.

Once we find out the oracles $O$ and $\ket{G}_a$, we can organize a unitary gate $W_q$, implemented by $\order{1}$ additional ancillary qubits and $\order{q}$-times usage of $O$, $G$, and $\order{n_a}$ elementary gates, which satisfies

\begin{equation}
W_q \ket{0}^{\otimes \left\{ n_a + \order{1} \right\}} \ket{\psi} = \ket{0}^{\otimes \left\{ n_a + \order{1} \right\} }  e^{-i H t} \ket{\psi} + \order{\varepsilon_q},
\end{equation}

for every state $\ket{\psi} \in \mcl{H}$.
Here, the error $\varepsilon_q$ comes from approximating $e^{-iHt}$ by a certain degree-$q$ polynomial of $H$, and decays as $\varepsilon_q \in \order{(\alpha t / q)^q}$.
The scaling of $q$ required for making the error $\varepsilon_q$ smaller than an acceptable error $\varepsilon$ is given by the following formula;

\begin{equation}\label{Eq2B:query_qubitization}
q \in \Theta \left( \alpha t + \frac{\log (1/\varepsilon)}{\log (e+(\alpha t)^{-1} \log (1/\varepsilon))}\right).
\end{equation}

In order to determine the explicit form of $W_q$, we resort to quantum signal processing <cit.>.
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 $\varepsilon$ (See also Table <ref>);

* Number of ancillary qubits; $n_a+\order{1}$

* Query complexity; $q$ [Eq. (<ref>)]

* Number of elementary gates;

\begin{equation}\label{Eq2B:qubitization_cost}
\order{\left\{ \alpha t + \frac{\log (1/\varepsilon)}{\log (e + (\alpha t)^{-1} \log (1/\varepsilon) )}\right\} (C+n_a)}.
\end{equation}

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

\begin{equation}\label{Eq2:qubitization_additive_scale}
\alpha t + o(\log (1/\varepsilon)) \in \poly{N} t + o(\log (1/\varepsilon)).
\end{equation}

The $o(\log (1/\varepsilon))$ term scales as $\log (1/\varepsilon) / \log \log (1/\varepsilon)$ under $1/\varepsilon \to \infty$ with fixed time $t$.
Its scaling is known to be optimal both in $t$ and $1/\varepsilon$ for simulating generic time-independent Hamiltonians <cit.>.

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. (<ref>), while the conventional truncated Dyson-series algorithm for generic time-dependent systems needs the query complexity in a productive form $\alpha t \times o(\log (\alpha t / \varepsilon))$.

§ SUMMARY OF THIS PAPER

In this section, we overview our main results on efficient quantum simulation of time-periodic Hamiltonians.
We will provide their detailed derivation in the following Sections <ref>-<ref>.

§.§ 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 $m_\imax \in \order{1}$;

\begin{equation}
H_m = 0, \quad \text{if} \quad |m| > m_\imax.
\end{equation}

(This discussion can be extended to cases where $H_m$ exponentially decays in $|m|$. See Appendix <ref>.).

The central attempt in the algorithm is to efficiently simulate Eq. (<ref>),

\begin{equation}\label{Eq3A:psi_t_infinite}
\ket{\psi(t)}  = \sum_{l = - \infty}^\infty e^{-il\omega t} \bra{l} e^{-i \ms{H}_\eff t} \ket{0} \ket{\psi(0)},
\end{equation}

on quantum circuits.
The time-dependency is erased by mapping the dynamics to the one on the Floquet-Hilbert space, which results in the following two benefits.
First, we do not discretize the time with infinitesimal intervals.
We instead rely on Fourier indices $l \in \bbZ$, which are originally discrete, and actually they lead to much better accuracy with the same 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/\varepsilon$, exploiting it for $\ms{H}_\eff$ accelerate the simulation of time-periodic systems.

On the other hand, in order to simulate the time-evolved state $\ket{\psi(t)}$ via Eq. (<ref>) 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 $l_\imax$, and consider a finite set of indices defined by

\begin{equation}
D^{l_\imax} = \{ -l_\imax +1, -l_\imax+2, \hdots, l_\imax \} \subset \bbZ.
\end{equation}

Then, we focus on the approximate dynamics by

\begin{eqnarray}
\ms{H}_\eff^{l_\imax} &=& \sum_{l \in D^{l_\imax}} \ket{l}\bra{l} \otimes (H_0-l\omega) \nonumber \\
&\quad& + \sum_{l \in D^{l_\imax}} \sum_{m \neq 0; l+m \in D^{l_\imax}} \ket{l}\bra{l+m} \otimes H_{-m}, \nonumber \\
&& \label{Eq3A:def_effective_H_l_max}\\
\ket{\psi^{l_\imax}(t)}  &=& \sum_{l \in D^{l_\imax}} e^{-il\omega t} \bra{l} e^{-i \ms{H}_\eff^{l_\imax} t} \ket{0}\ket{\psi(0)}. \label{Eq3A:def_psi_l_max}
\end{eqnarray}

The finite-dimensional space, spanned by $2l_\imax \times \mr{dim}(\mcl{H})$ states $\{ \ket{l} \ket{\psi_j} \, | \, l \in D^{l_\imax}, \, j=1,\hdots,\mr{dim}(\mcl{H})\}$, is called the truncated Floquet-Hilbert space.
While $l_\imax \to \infty$ reproduces $\ket{\psi^{l_\imax}(t)} \to \ket{\psi(t)}$, it is nontrivial how we should choose $l_\imax$ to achieve the accuracy $\norm{\ket{\psi(t)} - \ket{\psi^{l_\imax}(t)}} \leq \varepsilon$.
The second problem is the small success probability of post-selecting the ancillary state.
Even if we succeed in approximation with finite $l_\imax$, Eq. (<ref>) requires a projection to a unnormalized ancillary state $\sum_{l \in D^{l_\imax}} \ket{l}$ for the state in the truncated Floquet-Hilbert space,

\begin{equation}\label{Eq3A:def_large_psi_l_max}
\ket{\Psi^{l_\imax}(t)} = e^{-i \omega t \sum_l \ket{l} \bra{l}} e^{-i \ms{H}_\eff^{l_\imax} t} \ket{0} \ket{\psi(0)}.
\end{equation}

In the actual computation, we need post-selection to an ancillary state given by

\begin{equation}\label{Eq3A:def_a_l_max}
\ket{a^{l_\imax}} = \frac{1}{\sqrt{2l_\imax}}\sum_{l \in D^{l_\imax}} \ket{l}.
\end{equation}

Success probability of the projection is given by $\norm{\braket{a^{l_\imax}|\Psi^{l_\imax}(t)}}^2 \simeq \braket{\psi(t)|\psi(t)}/(2l_\imax) = 1/(2l_\imax)$.
As $l_\imax$ increases to ensure the accuracy, the expected time to successfully obtain $\ket{\psi(t)}$ becomes longer in proportion to it.
The final problem is about implementation of $e^{-i\ms{H}_\eff^{l_\imax}t}$.
While the effective Hamiltonian $\ms{H}_\eff$ is time-independent, its structure is complicated due to the additional degrees of freedom labeled by $\{\ket{l}\}$.
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 the time-evolved state $\ket{\psi(t)}$ with resolving the above problems.
The significant developments are composed of the following steps;

* Decision of the Fourier index for truncation, $l_\imax$, to achieve an allowable error $\varepsilon$.

* Amplification by symmetry and oblivious amplitude amplification to enhance the success probability up to $1-\order{\varepsilon}$.

* Efficient implementation of $\exp ( -i \ms{H}_\eff^{l_\imax} t)$ by the qubitization technique.

In Step (a), we explicitly derive the upper bound of $\norm{\ket{\psi(t)}-\ket{\psi^{l_\imax}(t)}}$ in a similar way to the Lieb-Robinson bound.
We show that the choice of $l_\imax$, based on

\begin{eqnarray}
l_\imax &\in& \Theta \left( \gamma t + \frac{\log (1/\varepsilon)}{\log ( e + (\gamma t)^{-1} \log (1/\varepsilon) )} \right), \\
\gamma &=& \sup_t (\norm{H(t)-H_0}) \label{Eq3A:def_gamma},
\end{eqnarray}

is sufficient to make the error smaller than $\varepsilon$.
Here, $\gamma$ gives the scale of the time-dependent terms in $H(t)$.
The second step (b) plays a role in amplifying the success probability of the post selection from $1/(2l_\imax)$ to $1-\order{\varepsilon}$.
The first protocol, which we call amplification by symmetry, exploits the symmetry of $\ms{H}_\eff$ which is always present and inherent in time-periodic systems.
It brings the success probability from $1/(2l_\imax)$ to $1/4-\order{\varepsilon}$ only with small cost of $\order{\log l_\imax}$ elementary gates.
Following this, we apply the oblivious amplitude amplification <cit.>, reminiscent of Grover's search algorithm <cit.>. The success probability is further amplified from $1/4-\order{\varepsilon}$ to $1-\order{\varepsilon}$.
Exploiting these two kinds of amplification, we can obtain the target state $\ket{\psi(t)}$ only with a little additional resource that does not change the scaling.
The final step (c) constructs conversion of the effective Hamiltonian so that $\exp ( -i \ms{H}_\eff^{l_\imax} t)$  can be efficiently implemented.
The original effective Hamiltonian $\ms{H}_\eff^{l_\imax}$ is not suitable for the qubitization technique, since the additional degrees of freedom $\ket{l}$ makes its block-encoding inefficient.
To resolve this problem, we derive an alternative effective Hamiltonian such that the desired evolution $\exp ( -i \ms{H}_\eff^{l_\imax} t)$ is accurately reproduced.
The block-encoding for the alternative effective Hamiltonian requires $\order{1}$-times queries to the oracles for each Fourier component $\{ H_m \}$, which is important for achieving optimal or nearly-optimal dependence in $t$ and $1/\varepsilon$.

§.§ Main results

Here, we summarize the computational resources for computing the dynamics $\ket{\psi(t)}$ 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 $\order{1}$ periods as $\omega t \in \order{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 $\omega t \in \Omega (1)$.
The algorithms for both cases follow Steps (a)-(c) of Section <ref>; the number of ancillary qubits is determined by the truncation order $l_\imax$ of Step (a), and the number of elementary gates comes mainly from the cost of implementing $\exp ( -i \ms{H}_\eff^{l_\imax} t)$  in Step (c).

We assume that the Fourier component $ H_m $ becomes zero for $|m|>m_\imax$ with an $\order{1}$ constant $m_\imax$. The block-encoding for each $H_m$ is supposed to be given by an oracle unitary gate $O_m$ and an $n_a$-qubit oracle state $\ket{G_m}$ as

\begin{equation}\label{Eq3B:def_block_encode_H_m}
\braket{G_m |O_m| G_m} = \frac{H_m}{\alpha_m}, \quad \alpha_m > 0.
\end{equation}

Each oracle state $\ket{G_m}$ is generated by trivial states as $\ket{G_m}=G_m\ket{0}$.
We define the energy scale of the Hamiltonian $H(t)$ by

\begin{equation}\label{Eq3B:def_alpha}
\alpha = \sum_{|m| \leq m_\imax} \alpha_m,
\end{equation}

and suppose that $\{ \alpha_m \}_m$ can be embedded into an $\order{1}$-qubit quantum system by

\begin{equation}\label{Eq3B:def_G_coef}
G_\mr{coef} \ket{0} = \sum_{|m| \leq m_\imax} \sqrt{\frac{\alpha_m}{\alpha}} \ket{m} \in \bbC^{2m_\imax+1}.
\end{equation}

The query complexity is defined by the number of queries to the oracles $O_m$, $G_m$, and $G_\mr{coef}$, where each of them is supposed to require at-most $C$ elementary gates.
We summarize the computational resources for simulating $\ket{\psi(t)}$ in the adiabatic-like regime and the generic long-time regime respectively by the following theorems.

(Resource for adiabatic-like regime)

Suppose that we are interested in the time-evolved state $\ket{\psi(t)}$ from an arbitrary initial state $\ket{\psi(0)}$ over $\order{1}$ periods.
The computational resources to obtain it with the acceptable error and the failure probability smaller than $\order{\varepsilon}$ are summarized as follows;

* Number of ancillary qubits;

\begin{equation}
n_a + \order{\log (\gamma t) + \log \log (1/\varepsilon)}.
\end{equation}

* Scaling of query complexity;

\begin{equation}\label{Eq3B:query_complexity_adiabatic}
\alpha t  + \frac{\log (1/\varepsilon)}{\log ( e + \{ \alpha t + o(\log (1/\varepsilon))\}^{-1} \log (1/\varepsilon) )}.
\end{equation}

The $o(\log (1/\varepsilon))$ term scales in the above as

\begin{equation}\label{Eq3B:o_log_term_adiabatic}
\frac{\log (1/\varepsilon)}{\log (e+ (\gamma t)^{-1} \log (1/\varepsilon))}.
\end{equation}

* Additional gates per query;

\begin{equation}
\order{n_a + \log (\gamma t) + \log \log (1/\varepsilon)}.
\end{equation}

The query complexity has optimal scaling in the time $t$ and nearly-optimal scaling in the inverse error $1/\varepsilon$.

(Resource for long-time regime)

Suppose that we are interested in the time-evolved state $\ket{\psi(t)}$ from arbitrary initial states $\ket{\psi(0)}$ over multiple periods $\omega t \in \Omega (1)$.
The computational resources to obtain it with the acceptable error and the failure probability smaller than $\order{\varepsilon}$ are summarized as follows;

* Number of ancillary qubits;

\begin{equation}
n_a + \order{\log (\gamma /\omega) + \log \log (\omega t/\varepsilon)}.
\end{equation}

* Scaling of query complexity;

\begin{equation}
\alpha t + \frac{ \omega t \log (\omega t/\varepsilon)}{\log (e+\{ \alpha /\omega + o(\log (\omega t /\varepsilon))\}^{-1} \log (\omega t/\varepsilon))}.
\end{equation}

The $\order{\log (\omega t/\varepsilon)}$ term scales as

\begin{equation}\label{Eq3B:o_log_term_long}
\frac{\log (\omega t /\varepsilon)}{\log (e + (\gamma / \omega)^{-1} \log (\omega t /\varepsilon))}.
\end{equation}

* Additional gates per query;

\begin{equation}
\order{n_a + \log (\gamma /\omega) + \log \log (\omega t /\varepsilon)}
\end{equation}

The scaling of the query complexity is optimal in time $t$ for practical problems up to $\poly{N}$-time, while it is formally nearly-optimal.
It is nearly-optimal in the inverse error $1/\varepsilon$.

The parameters $\alpha$ and $\gamma$, defined by Eqs. (<ref>) and (<ref>), respectively represent energy scales of the overall terms and the time-dependent terms in $H(t)$.
They typically scale as $\alpha, \gamma \in \mr{poly}(N)$ with the system size $N$.
In contrast, the frequency $\omega=2 \pi / T$ is typically an  $\order{N^0}$ value much smaller than $\alpha$ and $\gamma$ (The high-frequency cases where $\omega$ are comparable to or larger than $\alpha, \gamma$ are trivial, as discussed in Appendix <ref>).
While the above theorems are about time-periodic Hamiltonians with vanishing Fourier components at $m_\imax \in \order{1}$ as $H_m= 0$ ($|m| > m_\imax$), we obtain similar results for Hamiltonians with exponentially-decaying Fourier components $\norm{H_m} \lesssim e^{-\order{|m|}}$.
See Appendix <ref> for its derivation.

§ 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 <ref>.
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 <ref>, and determine a proper truncation order based on it in Section <ref>.

§.§ Lieb-Robinson bound in Floquet-Hilbert space

The proper choice of the truncation order $l_\imax$ is determined so that the approximate state $\ket{\psi^{l_\imax (t)}}$ [See Eq. (<ref>)] can reproduce the time-evolved state $\ket{\psi(t)}$ with $\norm{\ket{\psi(t)}-\ket{\psi^{l_\imax}(t)}} \leq \varepsilon$.
From the explicit formula Eq. (<ref>), we can see that it is important to observe how the transition amplitude $\braket{l|e^{-i \ms{H}_\eff^{l_\imax}t}|0}$ behaves for sufficiently large $l$.
We first show the upper bound on this transition amplitude in Theorem \ref{Thm4:Bound_transition_rate}.
Since it is reminiscent of the Lieb-Robinson bound in single-particle quantum systems \cite{Gong2019bound}, we call it the Lieb-Robinson bound in the Floquet-Hilbert space.

\begin{theorem}\label{Thm4:Bound_transition_rate}
\textbf{(Bound on transition amplitude)}

We assume $H_m = 0$ for $|m| > m_\imax$. Then, for indices $l, l^\prime \in \bbZ$ satisfying $|l-l^\prime| \geq 2 m_\imax \gamma t$, the transition amplitude is bounded by
\begin{eqnarray}
\norm{\bra{l} e^{-i \ms{H}_\eff^{l_\imax} t} \ket{l^\prime}} &\leq& 2 \frac{(\gamma t)^{\lceil |l-l^\prime| /m_\imax \rceil}}{(\lceil |l-l^\prime| /m_\imax \rceil)!} \label{Eq4A:thm_Lieb_Robinson_1} \\
&\leq& \left( \frac{\gamma t}{\lceil |l-l^\prime| /m_\imax \rceil}\right)^{\lceil |l-l^\prime| /m_\imax \rceil}, \label{Eq4A:them_Lieb_Robinson_2}\nonumber \\
\end{eqnarray}
where the parameter $\gamma$ is defined by Eq. (\ref{Eq3A:def_gamma}).
\end{theorem}

\textbf{Proof.---}
We will omit the superscripts $l_$ for some operators introduced here since they are not important. We first employ the interaction picture. With the unitary operation defined by
\begin{equation}\label{Eq4A:def_reference_hamiltonian}
\ms{U}_0(t) = e^{-i \ms{H}_0 t}, \quad \ms{H}_0 = \sum_{l \in D^{l_\imax}} \ket{l}\bra{l} \otimes (H_0-l \omega ),
\end{equation}
the time evolution operator $e^-iH_^l_t$ is represented as
\begin{eqnarray}
e^{-i\ms{H}_\eff^{l_\imax} t} &=& \ms{U}_0(t) \ms{U}_I (t), \\
\ms{U}_I(t) &=& \mcl{T} \exp \left( -i \int_0^t \dd t^\prime \ms{H}_I(t^\prime )\right).
\end{eqnarray}
Here, the Hamiltonian in the interaction picture is defined by
\begin{eqnarray}
\ms{H}_I(t) &=& \ms{U}_0(t)^\dagger (\ms{H}_\eff^{l_\imax}-\ms{H}_0) \ms{U}_0(t) \nonumber \\
&=& \sum_{l} \sum_{m>0} (e^{im\omega t}\ket{l}\bra{l+m} \otimes H_{-m}^I(t) + \text{h.c.} ), \nonumber \\
&& \label{Eq4B:effective_H_int_pic}\\
H_{-m}^I(t) &=& e^{i H_0 t} H_{-m} e^{-iH_0 t},
\end{eqnarray}
where the summation in the second line is taken over $l,m$ such that $l, l+m ∈D^l_$.
Using the Dyson series expansion, we obtain
\begin{eqnarray}
&& \norm{\bra{l} e^{-i \ms{H}_\eff^{l_\imax} t} \ket{l^\prime}} = \norm{\braket{l|\ms{U}_I(t)|l^\prime}} \nonumber \\
&\quad& \leq \sum_{n=0}^\infty \int_0^t \dd t_n \hdots \int_0^{t_2} \dd t_1 \norm{\braket{l|\ms{H}_I(t_n) \hdots \ms{H}_I(t_1)|l^\prime}}. \nonumber \\
&\quad& = \sum_{n=0}^\infty \int_0^t \dd t_n \hdots \int_0^{t_2} \dd t_1 \norm{\sum_{\{l_i\}}\prod_{i=1}^n \braket{l_i|\ms{H}_I(t_i)|l_{i-1}}}. \nonumber \\
&& \label{Eq4A:Dyson_transition_path}
\end{eqnarray}
In the last equality, we employ the identity $∑_l_i ∈D^l_ |l_i⟩⟨l_i|=I$ for $n-1$ times.
The summation $∑_{l_i }$ is taken over $l_i ∈D^l_$ for $i=1,2,,n-1$ under the fixed $l_0=l^'$ and $l_n=l$.
Each product $∏_i=1^n ⟨l_i|H_I(t_i)|l_i-1|$⟩ represents a complex transition amplitude from $\ket{l^\prime}$ to $\ket{l}$ along the path $\ket{l^\prime} \to \ket{l_1} \to \hdots \to \ket{l_{n-1}} \to \ket{l}$ under $\ms{H}_I(t_i)$.
Since the Hamiltonian $\ms{H}_I(t_i)$ shifts the Fourier index $\ket{l_i}$ by at-most $m_\imax$ due to Eq. (<ref>), the low order terms labeled by $n < |l-l^\prime|/m_\imax$ disappears.
This results in

\begin{eqnarray}
&& \norm{\bra{l} e^{-i \ms{H}_\eff^{l_\imax} t} \ket{l^\prime}} \nonumber \\
&\quad& \leq \sum_{n=\lceil |l-l^\prime|/m_\imax \rceil}^\infty \int_0^t \dd t_n \hdots \int_0^{t_2} \dd t_1 \norm{\braket{l|\prod_{i=1}^n \ms{H}_I(t_i) |0}} \nonumber \\
&\quad& \leq \sum_{n=\lceil |l-l^\prime|/m_\imax \rceil}^\infty \frac{t^n}{n!} \left( \sup_t (\norm{\ms{H}_I(t)})\right)^n
\end{eqnarray}

Here, $\ms{H}_I(t)$ is a Toeplitz matrix that satisfies $\braket{l|\ms{H}_I(t)|l^\prime} = e^{-i(l-l^\prime)\omega t} H_{l-l^\prime}^I(t) - H_0 \delta_{l-l^\prime,0}$.
This leads to the upper bound of its operator as follows <cit.>,

\begin{eqnarray}
\norm{\ms{H}_I(t)} &\leq& \sup_{t^\prime} \left( \norm{\sum_{m \in \bbZ} (e^{im\omega t} H_{-m}^I(t) - H_0 \delta_{m,0}) e^{im\omega t^\prime}}\right) \nonumber \\
&=& \sup_{t^\prime} \left( \norm{H(t-t^\prime) - H_0} \right) = \gamma. \label{Eq4A:H_eff_bound_gamma}
\end{eqnarray}

We also use the following inequality;

\begin{equation}\label{Eq4A:sum_Taylor}
\sum_{n=n_0}^\infty \frac{x^n}{n!} \leq 2 \frac{x^{n_0}}{n_0 !}, \quad \text{if} \quad  n_0 \geq 2x \geq 0.
\end{equation}

This relation can be easily confirmed by

\begin{eqnarray}
\sum_{n=n_0}^\infty \frac{x^n}{n!} &=& \frac{x^{n_0}}{n_0 !} \sum_{n=n_0}^\infty \frac{x^{n-n_0}}{n (n-1) \hdots (n_0+1)} \nonumber \\
&\leq&  \frac{x^{n_0}}{n_0 !} \sum_{n=0}^\infty \left(\frac12\right)^n.
\end{eqnarray}

For the indices $l, l^\prime$ satisfying $|l-l^\prime| \geq 2m_\imax \gamma t$, we can substitute $x=\gamma t$ and $n_0=\lceil |l-l^\prime| /m_\imax \rceil$. This results in the bound,

\begin{equation}
\norm{\bra{l} e^{-i \ms{H}_\eff^{l_\imax} t} \ket{l^\prime}} \leq 2 \frac{(\gamma t)^{\lceil |l-l^\prime| /m_\imax \rceil}}{(\lceil |l-l^\prime| /m_\imax \rceil)!}.
\end{equation}

The inequality from the Stirling formula,

\begin{equation}\label{Eq4A:Stirling_inequality}
n_0 ! \geq 2 \left( \frac{n_0}{e}\right)^{n_0},
\end{equation}

ensures the inequality Eq. (<ref>). $\quad \square$

The Lieb-Robinson bound in the transition amplitude provides a guide for choosing a proper truncation order.
By setting $l^\prime=0$, Eqs. (<ref>) and (<ref>) say that the contributions to $\ket{\psi^{l_\imax}(t)}$ from indices satisfying $|l| \geq 2 m_\imax \gamma t$ rapidly decay as $l^{-l}$.
Thus, $\order{\gamma t}$ becomes a possible truncation order.
However, we note that the change in $l_\imax$ affects $\ket{\psi^{l_\imax}(t)}$ also via the change in the support of the effective Hamiltonian $\ms{H}_\eff^{l_\imax}$.
By taking it into account with using the Lieb-Robinson bound, we obtain the exact upper bound on the error between $\ket{\psi(t)}$ and $\ket{\psi^{l_\imax}(t)}$ as follows.

(Floquet-Hilbert space truncation)

We consider the approximate time-evolved state obtained from the truncated Floquet-Hilbert space, $\ket{\psi^{l_\imax}(t)}$, given by Eq. (<ref>).
Then, its deviation from the exact one $\ket{\psi}$ is bounded by

\begin{equation}\label{Eq4B:deviation_truncation}
\norm{\ket{\psi(t)} - \ket{\psi^{l_\imax}(t)}} \leq 20 m_\imax \frac{(\gamma t)^{\lceil l_\imax/m_\imax \rceil}}{(\lceil l_\imax/m_\imax \rceil)!},
\end{equation}

if the truncation order $l_\imax$ satisfies $l_\imax \geq 2 m_\imax \gamma t$.

We evaluate the convergence of $\ket{\psi^{l_\imax}(t)}$.
For different orders $l_\imax, l_\imax^\prime$ satisfying $l_\imax^\prime > l_\imax$, we compute the difference,

\begin{equation}\label{Eq4A:thm_truncation_proof_1}
\norm{\ket{\psi^{l_\imax^\prime}(t)}-\ket{\psi^{l_\imax}(t)}} \leq \varepsilon_1 +\varepsilon_2,
\end{equation}

\begin{eqnarray}
\varepsilon_1 &=& \sum_{ l \in (D^{l_\imax^\prime} \backslash D^{l_\imax})} \norm{\bra{l} e^{i\ms{H}_\eff^{l_\imax^\prime}t} \ket{0}}, \label{Eq4A:thm_truncation_proof_2}\\
\varepsilon_2 &=& \sum_{l \in D^{l_\imax}} \norm{\bra{l} e^{i\ms{H}_\eff^{l_\imax^\prime}t}-e^{i\ms{H}_\eff^{l_\imax}t} \ket{0}}. \label{Eq4A:thm_truncation_proof_3}
\end{eqnarray}

The first error comes from truncating the order of the post-selected state $\ket{a^{l_\imax}}$.
Using Theorem <ref> directly implies

\begin{eqnarray}
\varepsilon_1 &\leq& \sum_{ l \in (D^{l_\imax^\prime} \backslash D^{l_\imax})} \frac{(\gamma t)^{\lceil |l|/m_\imax \rceil}}{(\lceil |l|/m_\imax \rceil)!} \nonumber \\
&\leq& 2 \sum_{l = l_\imax}^\infty  \frac{(\gamma t)^{\lceil l/m_\imax \rceil}}{(\lceil l/m_\imax \rceil)!} .
\end{eqnarray}

The summation $\sum_{l=l_\imax}^\infty$ can be divided based on $l$ $(\text{mod. } m_\imax)$, and each summation corresponds to the left hand side of Eq. (<ref>) with $n_0 \geq \lceil l_\imax / m_\imax \rceil$.
In other words, for $l_\imax \geq 2 m_\imax \gamma t$, we obtain the following inequality,

\begin{equation}\label{Eq4A:sum_Taylor_modulo}
\sum_{l = l_\imax}^\infty  \frac{(\gamma t)^{\lceil l/m_\imax \rceil}}{(\lceil l/m_\imax \rceil)!} \leq 2 m_\imax \frac{(\gamma t)^{\lceil l_\imax/m_\imax \rceil}}{(\lceil l_\imax/m_\imax \rceil)!}.
\end{equation}

As a result, the first error $\varepsilon_1$ is bounded by $4 m_\imax (\gamma t)^{\lceil l_\imax/m_\imax \rceil}/(\lceil l_\imax/m_\imax \rceil)!$.

The second error $\varepsilon_2$ comes from truncating the order of the effective Hamiltonian.
In a similar way to the proof of Theorem <ref>, each term is bounded by

\begin{eqnarray}
&& \norm{\bra{l} e^{i\ms{H}_\eff^{l_\imax^\prime}t}-e^{i\ms{H}_\eff^{l_\imax}t} \ket{0}} \nonumber \\
&\quad& \leq \sum_{n=0}^\infty \int_0^t \dd t_n \hdots \int_0^{t_2} \dd t_1 \nonumber \\
&& \qquad \quad \norm{\braket{l|\prod_{i=1}^n \ms{H}_I^{l_\imax^\prime}(t_i) - \prod_{i=1}^n \ms{H}_I^{l_\imax}(t_i)|0}},
\end{eqnarray}

where we explicitly write the size $l_\imax,l_\imax^\prime$ for $H_I(t)$ [See Eq. (<ref>)].
In a similar manner to Eq. (<ref>), $\braket{l|\prod_{i=1}^n \ms{H}_I^{l_\imax^\prime}(t_i)|0}$ and $⟨l|∏_i=1^n H_I^l_(t_i)|0|$⟩ respectively represent transition amplitudes summed over all possible paths $\ket{0} \to \ket{l_1} \to \hdots \to \ket{l_{n-1}} \to \ket{l}$ under $\ms{H}_\eff^{l_\imax^\prime}(t)$ and $\ms{H}_\eff^{l_\imax}(t)$.
Their differences survive only when the path goes across the domain $D^{l_\imax^\prime} \backslash D^{l_\imax}$, where these two Hamiltonians have different actions.
In other words, low order terms with $n < \{ (l_\imax-0) + (l_\imax - |l|)/m_\imax$ should vanish.
Considering that the norm of $\ms{H}_\eff^{l_\imax}(t)$ and $\ms{H}_\eff^{l_\imax^\prime}(t)$ is both bounded by $\gamma$, we obtain

\begin{eqnarray}
\varepsilon_2 &\leq& \sum_{l \in D^{l_\imax}} \sum_{n=\lceil (2l_\imax-|l|)/m_\imax \rceil}^\infty \frac{t^n}{n!} 2 \gamma^n \nonumber \\
&\leq& 4 \sum_{l \in D^{l_\imax}} \frac{(\gamma t)^{\lceil (2l_\imax-|l|)/m_\imax \rceil}}{(\lceil (2l_\imax-|l|)/m_\imax \rceil)!} \nonumber \\
&\leq& 16 m_\imax \frac{(\gamma t)^{\lceil l_\imax/m_\imax \rceil}}{(\lceil l_\imax/m_\imax \rceil)!}.
\end{eqnarray}

In the last inequality, we again use the relation Eq. (<ref>).
Taking the limit $l_\imax^\prime \to \infty$ for $\varepsilon_1 + \varepsilon_2$ reproduces the result of Theorem <ref>. $\quad \square$

§.§ Truncation order of Floquet-Hilbert space

We hereby determine the truncation order $l_\imax$ so that $\ket{\psi^{l_\imax}(t)}$ can reproduce the exact time-evolved state with a desirable error up to $\order{\varepsilon}$.
By using the inequality from the Stirling formula, Eq. (<ref>) to Theorem <ref>, we obtain the error bounded by

\begin{equation}
\norm{\ket{\psi(t)}-\ket{\psi^{l_\imax}(t)}} \leq 10 m_\imax \left( \frac{em_\imax \gamma t}{l_\imax}\right)^{l_\imax/m_\imax},
\end{equation}

for $l_\imax \geq 2 m_\imax \gamma t$.
The truncation order $l_\imax$ is chosen so that the right hand side can be smaller than $\varepsilon$.
To this aim, we should evaluate a function $f(x)=(\kappa/x)^x$ [$\kappa > 0$], which is known to be dealt with the Lambert W function $W(x)$ satisfying $W(x)e^{W(x)}=x$ <cit.>.
Here, we rely on the resulting proposition <cit.>.

The function $f(x)=(\kappa/x)^x$ [$\kappa>0$] is monotonically decreasing in $x \geq \kappa /e$, and satisfies the following inequality for $0 < \eta < 1$;

\begin{equation}
f(x) \leq \eta, \quad^\forall x \geq e \kappa + \frac{4 \log (1/\eta)}{\log (e + \kappa^{-1} \log (1/\eta))}.
\end{equation}

See Lemma 59 in the full version of Ref. <cit.> for the proof.
Based on the above proposition, we choose $l_\imax$ by

\begin{eqnarray}
&& l_\imax - m_\imax =  \nonumber \\
&& \quad \left\lceil e^2 m_\imax \gamma t + \frac{4 m_\imax \log (10m_\imax / \varepsilon )}{\log (e + (e\gamma t)^{-1} \log (10m_\imax / \varepsilon))}\right\rceil, \nonumber \\
&& \label{Eq4B:l_max_choice}
\end{eqnarray}

so that the error can be bounded from above as

\begin{equation}\label{Eq4B:error_condition}
10 m_\imax \left( \frac{em_\imax \gamma t}{l_\imax-m_\imax}\right)^{(l_\imax-m_\imax)/m_\imax} \leq \varepsilon.
\end{equation}

We note that this choice does not violate the condition $l_\imax \geq 2 m_\imax \gamma t$, which is required for Theorem <ref>.
The monotonicity of $f(x)$ ensures $\norm{\ket{\psi(t)}-\ket{\psi^{l_\imax}(t)}} \leq \varepsilon$ under the above choice (The additional term $m_\imax$ is attached for later calculation, especially for Appendix <ref>) and <ref>).

Let us discuss how $l_\imax$ increases in the time $t$ and the inverse error $1/\varepsilon$.
The form of the error in $l_\imax$, given by Eq. (<ref>), is the same as that for qubitization in the query complexity $q$, given by $\varepsilon_q \in \order{(\alpha t /q)^q}$.
Therefore, its scaling can be determined in a similar way, which results in

\begin{equation}\label{Eq4B:l_max_scaling}
l_\imax \in \Theta \left( \gamma t + \frac{\log (1/\varepsilon)}{\log ( e+ (\gamma t)^{-1} \log (1/\varepsilon) )} \right).
\end{equation}

To reproduce the truncated Floquet-Hilbert space, we should prepare an ancillary system labeled by $\{ \ket{l} \}_{l\in D^{l_\imax}}$.
The number of qubits for such ancillary system is at-most

\begin{equation}
\order{\log l_\imax} \subset \order{\log (\gamma t) + \log \log (1/\varepsilon)}.
\end{equation}

§ AMPLITUDE AMPLIFICATION OF THE DYNAMICS

In this section, we show two kinds of amplification to achieve sufficiently high success probability of extracting the time-evolved state $\ket{\psi(t)}$.

As we stated in Section <ref>, the approximate state $\ket{\psi^{l_\imax}(t)}$ cannot be directly realized on quantum circuits.
After simulating the dynamics in the truncated Floquet-Hilbert space to get $\ket{\Psi^{l_\imax}(t)}$ defined by Eq. (<ref>), we make a projection to $\ket{a^{l_\imax}}$.
Although the resulting renormalized state is sufficiently close to $\ket{\psi^{l_\imax}(t)}$ and also the time-evolved state $\ket{\psi(t)}$ as

\begin{equation}
\frac{\braket{a^{l_\imax}|\Psi^{l_\imax}(t)}}{ \norm{\braket{a^{l_\imax}|\Psi^{l_\imax}(t)}}} = \frac{\ket{\psi^{l_\imax}(t)}}{\norm{\ket{\psi^{l_\imax}(t)}}} = \ket{\psi(t)} + \order{\varepsilon}
\end{equation}

for $l_\imax$ given by Eq. (<ref>), we should be careful of the low success probability; it becomes $\norm{\braket{a^{l_\imax}|\Psi^{l_\imax}(t)}}^2 \in \order{(l_\imax)^{-1}}$.

We need approximately $\order{\gamma t}$-times trials of the post selection, and every trial is expected to require at-least $\order{t}$ complexity for implementing $\ket{\Psi^{l_\imax}(t)} = e^{-i \ms{H}_\eff^{l_\imax} t} \ket{0}\otimes \ket{\psi(0)}$.
Therefore, the naive implementation based on Eq. (<ref>) is not efficient for the time-evolved state $\ket{\psi(t)}$ in that the expected computational time reaches $\order{t^2}$.

We resolve this problem by the amplification of the success probability up to $1-\order{\varepsilon}$ in the next section.
The first one, which exploits the symmetry of the effective Hamiltonian $\ms{H}_\eff$, amplifies from $\order{l_\imax^{-1}}$ to $\order{1}$.
The latter one following this, which is reminiscent of the Grover's search algorithm, allows the success probability $1-\order{\varepsilon}$.
As discussed later, using only either one fails to efficiently compute $\ket{\psi(t)}$.

§.§ Amplification by symmetry

We introduce the amplification exploiting the symmetry of the effective Hamiltonian $\ms{H}_\eff$.
First of all, we specify the symmetry here; it is about the translation of the photon number from $\ket{l}$ to $\ket{l + m}$.
When we define the translation operator on the Floquet-Hilbert space by $\ms{T}_m = \sum_{l \in \bbZ} \ket{l}\bra{l+m} \otimes I$, $\ms{H}_\eff$ satisfies the translation symmetry,

\begin{eqnarray}
\ms{T}_m^\dagger \ms{H}_\eff \ms{T}_m &=& \ms{H}_\eff + m\omega, \label{Eq5A:symmetry_effective_H} \\
\braket{l|e^{-i \ms{H}_\eff t}|l^\prime} &=& e^{il^\prime \omega t} \braket{l-l^\prime|e^{-i \ms{H}_\eff t}|0}. \label{Eq5A:symmetry_transition_amplitude}
\end{eqnarray}

Since the symmetry is present as long as the time-periodicity $H(t+T)=H(t)$ holds, the amplification 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 $\bbC^{8l_\imax} \otimes \mcl{H}$, where the Fourier index $\ket{l}$ is chosen from $l \in D^{4l_\imax} = \{-4l_\imax+1, -4 l_\imax+2, \hdots, 4l_\imax\}$.
The number of additional qubits required for this extension from $\bbC^{l_\imax} \otimes \mcl{H}$ is three, and hence this does not affect the scaling of the computational resources.
The second task is to preprocess the initial state by

\begin{equation}
\ms{U}_\mr{ini}^{l_\imax} \ket{0}\ket{\psi(0)} = \ket{a^{l_\imax}} \ket{\psi(0)},
\end{equation}

where the unitary operator $\ms{U}_\mr{ini}^{l_\imax}$ is expressed by

\begin{equation}
\ms{U}_\mr{ini}^{l_\imax} = \left( \ket{a^{l_\imax}} \bra{0} + \hdots \right) \otimes I.
\end{equation}

The unitary operator $\ms{U}_\mr{ini}^{l_\imax}$ nontrivially acts only on the ancillary space, and it can be implemented with $\order{\log l_\imax}$ gates
(e.g. for $l_\imax$ such that $\log_2 l_\imax \in \bbN$, we have $\ms{U}_\mr{ini}^{l_\imax} = \mr{Had}^{\otimes \log_2 (2l_\imax)} \otimes I$ with the Hadamard gate $\mr{Had}$).
This process plays a role of making the initial state approximately translation-invariant with the width of Fourier indices, $2l_\imax$.
The remaining procedures after the above two tasks is the same as those for the original protocol in Section <ref>.
Reflecting that the ancillary Hilbert space is extended to $\bbC^{8l_\imax}$, we evolve the above uniform state $\ket{a^{l_\imax}} \ket{\psi(0)}$ by $\ms{H}_\eff^{4l_\imax}$ and $\sum_{l \in D^{4l_\imax}} l\omega \ket{l} \bra{l} \otimes I$, and make a projection to $\ket{a^{4l_\imax}}$.
The resulting state $\ket{\bar{\psi}^{l_\imax}(t)}$ is calculated as follows.

\begin{equation}\label{Eq5A:dynamics_uniform_initial_state}
\ket{\bar{\psi}^{l_\imax}(t)} \equiv \bra{a^{4l_\imax}} e^{-it \sum_l l \omega \ket{l}\bra{l} }e^{-i\ms{H}_\eff^{4l_\imax}t} \ket{a^{l_\imax}} \ket{\psi(0)}.
\end{equation}

The point of this process is that the above state $\ket{\bar{\psi}^{l_\imax}(t)}$ can reproduce $\ket{\psi(t)}$ with $\order{1}$ amplitude owing to the translation symmetry of the effective Hamiltonian, Eqs. (<ref>) and (<ref>).
We give its brief explanation in this section, while the rigorous derivation is provided in Appendix <ref>.
Although we do not have the exact translation symmetry represented by Eqs. (<ref>) due to the finite-size effect of $l_\imax$, it is expected to approximately hold as

\begin{equation}\label{Eq5A:approximate_symmetry}
\bra{l} e^{-i\ms{H}_\eff^{4l_\imax}t} \ket{l^\prime} \simeq e^{il^\prime \omega t}\bra{l-l^\prime} e^{-i\ms{H}_\eff^{4l_\imax}t} \ket{0}.
\end{equation}

if the truncation order $l_\imax$ is sufficiently large.
As a matter of fact, we can derive the exact upper bound on the difference between the left- and right-hand sides based on the Lieb-Robinson bound (See Lemma <ref> in Appendix <ref>).
We proceed the discussion with assuming the approximate relation Eq. (<ref>) here.
The resulting state of the process $\ket{\bar{\psi}^{l_\imax}(t)}$ can be roughly computed as follows;

\begin{eqnarray}
&& \ket{\bar{\psi}^{l_\imax}(t)} \nonumber \\
&&  = \frac{1}{4l_\imax}\sum_{l^\prime \in D^{l_\imax}}  \sum_{l \in D^{4l_\imax}} e^{-il\omega t} \bra{l} e^{-i\ms{H}_\eff^{4l_\imax}t} \ket{l^\prime} \ket{\psi(0)} \nonumber \\
&& \simeq \frac{1}{4l_\imax}\sum_{l^\prime \in D^{l_\imax}}   \sum_{l \in D^{4l_\imax}} e^{-i(l-l^\prime)\omega t} \nonumber \\
&& \qquad \qquad \qquad \qquad \qquad \quad \times \bra{l-l^\prime} e^{-i\ms{H}_\eff^{4l_\imax}t} \ket{0} \ket{\psi(0)} \nonumber \\
&& = \frac{1}{4l_\imax}\sum_{l^\prime \in D^{l_\imax}} (\ket{\psi(t)} + \order{\varepsilon}) \nonumber \\
&& = \frac12 \ket{\psi(t)} + \order{\varepsilon}. \label{Eq5A:amp1_calculation}
\end{eqnarray}

The second equality comes from Theorem <ref>, considering that the summation of $l-l^\prime$ over $l \in D^{4l_\imax}$ is sufficient to suppress the error up to $\order{\varepsilon}$ under $l_\imax \in \order{\gamma t + \log (1/\varepsilon)/\log \log (1/\varepsilon)}$.
As a result, we obtain

\begin{equation}
\ket{\bar{\psi}^{l_\imax}(t)} \simeq  \frac{1}{2}(\ket{\psi(t)}+\order{\varepsilon}).
\end{equation}

Let us focus on the success probability of the projection onto $\ket{a^{4l_\imax}}$.
It is provided by $\braket{\bar{\psi}^{l_\imax}(t)|\bar{\psi}^{l_\imax}(t)} = 1/4 + \order{\varepsilon}$, which is much larger than the original one $1/(2l_\imax)$.
Therefore, the amplification protocol by symmetry enables us to avoid $\order{t}$-times repetition of the time evolution $e^{-i \ms{H}_\eff t}$ in contrast to the original protocol.
Intuitively, this drastic improvement can be understood as a result of the interference under the translation symmetry, as Fig. <ref>.
As we can see from Eq. (<ref>), the target state $\ket{\psi(t)}$ is extracted from $\ket{\Psi^{l_\imax}(t)}$ via the accompanying ancillary state $\ket{a^{l_\imax}}$, which is uniform in the Fourier index $\ket{l}$.
When we begin with the non-uniform initial state $\ket{0}\ket{\psi(0)}$, the resulting state after the time evolution is also non-uniform.
Since it involves $(1/\sqrt{2l_\imax}) \sum_{l \in D^{l_\imax}} e^{ikl} \ket{l}$ for $k \in (2\pi / 2l_\imax) \bbZ$ (different eigenstates of the translation operator) with approximate equal weight, the amplitude of the desirable component $\ket{a^{l_\imax}}$ is comparably small as $1/\sqrt{2l_\imax}$.
In contrast, when we employ the uniform initial state $\ket{a^{l_\imax}}\ket{\psi(0)}$, only the uniform components that includes $\ket{a^{l_\imax}}$ are amplified while the other components cancels one another.
It can be viewed also as the interference of different initial states $\ket{l^\prime}\ket{\psi(0)}$ in $\ket{a^{l_\imax}}\ket{\psi(0)}$ as shown in Fig. <ref>, and this is why the amplification protocol achieves $\order{1}$ amplitude of $\ket{a^{l_\imax}}\ket{\psi(0)}$.

We summarize the amplification by symmetry with adding the exact results obtained in Appendix <ref>.
We prepare the truncated Floquet-Hilbert space $\bbC^{8 l_\imax} \otimes \mcl{H}$, and make the initial state uniform as $\ket{a^{l_\imax}}\ket{\psi(0)}$.
As a result, this process is described by a unitary gate

\begin{equation}\label{Eq5A:U_amp1_def}
\ms{U}_\mr{amp1}^{l_\imax} (t) = (\ms{U}_\mr{ini}^{4l_\imax})^\dagger e^{-i t \sum_l l \omega \ket{l}\bra{l}} e^{-i \ms{H}_\eff^{4l_\imax} t} \ms{U}_\mr{ini}^{l_\imax},
\end{equation}

where the last unitary gate $(\ms{U}_\mr{ini}^{4l_\imax})^\dagger$ is added to replace the projection to $\ket{a^{4l_\imax}}$ by the one to $\ket{0}$.
This approximately realizes the time-evolved state $\ket{\psi(t)}$ with $1/4+\order{\varepsilon}$ success probability as

\begin{equation}\label{Eq5A:approximate_amp1_result}
\bra{0} \ms{U}_\mr{amp1}^{l_\imax} (t) \ket{0} \ket{\psi(0)} \simeq \frac{1}{2} (\ket{\psi(t)}+\order{\varepsilon}).
\end{equation}

The probability $1/4$ comes from the ratio of the width of the initial state $\ket{a^{l_\imax}}$ to that of the projected ancillary state $\ket{a^{4l_\imax}}$.
The exact version, derived based on the Lieb-Robinson bound in Appendix <ref>, is stated as follows, and ensures the validity of the discussion here.

(Amplification by symmetry)

We designate the truncation order $l_\imax \in \Theta (\gamma t + \log (1/\varepsilon) / \log \log (1/\varepsilon))$ by Eq. (<ref>). Then,

\begin{equation}
\norm{\bra{0} \ms{U}_\mr{amp1}^{l_\imax} (t) \ket{0} \ket{\psi(0)} - \frac{1}{2}\ket{\psi(t)}} \leq \frac{\varepsilon}{3}.
\end{equation}

is satisfied for an arbitrary initial state $\ket{\psi(0)} \in \mcl{H}$.

We finally note that the amplification solely by this process cannot achieve the success probability $1-\order{\varepsilon}$ 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 $p l_\imax$, $q l_\imax$, $q l_\imax$ ($p, q \in \bbN$, $q \geq p +2$) respectively, the success probability becomes $p/q$.
However, the probability larger than $1 - \varepsilon$ demands the relation,

\begin{equation}
1- \varepsilon < \frac{p}{q} \leq \frac{p}{p+2}.
\end{equation}

This implies $p \geq 2(1-\varepsilon)/\varepsilon$, and hence the number of ancillary qubits reaches at-least $\log (p l_\imax) \in \order{\log (\gamma t/\varepsilon)}$.
In addition, the query complexity includes a term proportional to the dimension of the ancillary system, $q l_\imax$ as discussed later in Section <ref>.
This implies that the query complexity increases linearly in $1/\varepsilon$, destroying the original logarithmic scaling.
Thus, relying solely on this amplification protocol is not suitable for efficient simulation of the dynamics.
We resolve this problem by the oblivious amplitude amplification below.

Schematic picture of how the dynamics in the truncated Floquet-Hilbert space provides the time-evolved state $\ket{\psi(t)}$. (a) The dynamics when we begin with the initial state $\ket{0}\ket{\psi(0)}$ [See Eq. (<ref>)]. Following the Lieb-Robinson bound, Theorem <ref>, the state spreads out within the region $|l| \lesssim l_\imax$. This spread results in the low amplitude $\order{1/\sqrt{l_\imax}}$ of the target output $\ket{a^{l_\imax}}\ket{\psi(t)}$. (b) The dynamics under the amplification by symmetry, where we employ the uniform initial state $\ket{a^{l_\imax}}\ket{\ket{\psi(0)}}$ as Eq. (<ref>). Interference takes place within the region $|l| \lesssim l_\imax$, and thereby enhances the amplitude of the uniform target state $\ket{a^{4l_\imax}}\ket{\psi(t)}$ up to $\order{1}$. The gray region represents the support of boundary terms in the refined effective Hamiltonian $\ms{H}_{\eff,\mr{pbc}}^{4l_\imax}$ [See Eq. (<ref>)]. 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 <ref> says.

§.§ Oblivious amplitude amplification

The above amplification based on the translation symmetry of $\ms{H}_\eff$ enhances the amplitude of $\ket{\psi(t)}$ from $1/\sqrt{2l_\imax}$ to $1/2-\order{\varepsilon}$.
Here, we introduce another amplification, called the oblivious amplitude amplification <cit.>, to achieve the amplitude (or the success probability) $1-\order{\varepsilon}$.

The starting point of this amplification is the result of the previous section <ref>.
Equation (<ref>), or equivalently Theorem <ref>, indicates that the consequence of the amplification by symmetry is written as

\begin{equation}\label{Eq5B:action_amp1}
\ms{U}_\mr{amp1}^{l_\imax} (t) \ket{0}\ket{\psi(0)} = \frac{1}{2} \ket{0}  (\ket{\psi(t)}+\order{\varepsilon}) + \ket{\Psi^\perp},
\end{equation}

with an additional term $\ket{\Psi^\perp} \in \bbC^{8l_\imax} \otimes \mcl{H}$ satisfying

\begin{equation}
(\ket{0}\bra{0} \otimes I)\ket{\Psi^\perp} = 0.
\end{equation}

The state $\ket{\Psi^\perp}$ generally depends on $\ket{\psi(0)}$ and $t$.
The oblivious amplitude amplification takes a similar strategy to that of the Grover's search algorithm, in which we compose the following two unitary operators;

\begin{eqnarray}
\ms{R} &=& (2\ket{0}\bra{0}-I) \otimes I, \label{Eq5B:def_reflection} \\
\ms{U}_\mr{amp2}^{l_\imax} (t) &=& - \ms{U}_\mr{amp1}^{l_\imax} (t) \ms{R} [\ms{U}_\mr{amp1}^{l_\imax} (t)]^\dagger \ms{R} \ms{U}_\mr{amp1}^{l_\imax} (t). \label{Eq5B:def_amp2}
\end{eqnarray}

The first one $\ms{R}$ reverses the sign of $\ket{l}$ for $l \neq 0$, and it is implemented with $\order{\log l_\imax}$ gates.
The second one $\ms{U}_\mr{amp2}^{l_\imax} (t)$ plays a role in enhancing the amplitude of $\ket{\psi(t)}$ up to $1-\order{\varepsilon}$.
Its action on any initial state $\ket{0}\ket{\psi(0)} $ is computed as follows;

\begin{eqnarray}
&& \ms{U}_\mr{amp2}^{l_\imax} \ket{0}\ket{\psi(0)}  \nonumber \\
&& = \ms{U}_\mr{amp1}^{l_\imax}\ms{R} (\ms{U}_\mr{amp1}^{l_\imax})^\dagger \left( - \frac{1}{2} \ket{0}{\psi(t)} + \ket{\Psi^\perp} \right) + \order{\varepsilon} \nonumber \\
&& = \ms{U}_\mr{amp1}^{l_\imax}\ms{R} \left\{ \ket{0}\ket{\psi(0)} - (\ms{U}_\mr{amp1}^{l_\imax})^\dagger \ket{0}\ket{\psi(t)}  \right\} + \order{\varepsilon}. \nonumber \\
&& \label{Eq5B:amp2_middle_1}
\end{eqnarray}

In the second equality, we use the relation obtained by applying $(\ms{U}_\mr{amp1}^{l_\imax})^\dagger$ to Eq. (<ref>).
Next, we evaluate

\begin{eqnarray}
&& \ms{R} (\ms{U}_\mr{amp1}^{l_\imax})^\dagger \ket{0,\psi(t)} \nonumber \\
& \quad& = 2 \ket{0} \left( \bra{0} (\ms{U}_\mr{amp1}^{l_\imax})^\dagger \ket{0} \right) \ket{\psi(t)} - (\ms{U}_\mr{amp1}^{l_\imax})^\dagger \ket{0,\psi(t)}. \nonumber \\
\end{eqnarray}

Theorem <ref> indicates that the time evolution operator $U(t)$ is approximated as

\begin{equation}
\norm{\braket{0|\ms{U}_\mr{amp1}^{l_\imax}(t)|0} - \frac{1}{2} U(t)} \leq \frac{\varepsilon}{3},
\end{equation}

and hence the relation

\begin{equation}
\norm{\braket{0|[\ms{U}_\mr{amp1}^{l_\imax}(t)]^\dagger|0} - \frac{1}{2} U(t)^\dagger} \leq \frac{\varepsilon}{3}
\end{equation}

is also satisfied.
This provides the relation,

\begin{eqnarray}
&& \ms{R} (\ms{U}_\mr{amp1}^{l_\imax})^\dagger \ket{0}\ket{\psi(t)} \nonumber \\
&& \qquad = \ket{0}\ket{\psi(0)} - (\ms{U}_\mr{amp1}^{l_\imax})^\dagger \ket{0}\ket{\psi(t)} + \order{\varepsilon},
\end{eqnarray}

and substituting this into Eq. (<ref>) results in

\begin{equation}
\ms{U}_\mr{amp2}^{l_\imax}(t) \ket{0} \ket{\psi(0)} = \ket{0} \ket{\psi(t)} + \order{\varepsilon},
\end{equation}

for an arbitrary initial state $\ket{\psi(0)}$.

This result indicates that the operation $\ms{U}_\mr{amp2}^{l_\imax}(t)$ generates the time-evolved state $\ket{\psi(t)}$ with the amplitude $1-\order{\varepsilon}$.
Or equivalently, as an exact bound for the error, we can derive the inequality,

\begin{equation}\label{Eq5B:deviation_amp2}
\norm{\ms{U}_\mr{amp2}^{l_\imax}(t) \ket{0} \ket{\psi(0)} - \ket{0}  \ket{\psi(t)}} \leq \varepsilon.
\end{equation}

The coefficient comes from the fact that an error bounded by $\varepsilon /3$ appears due to Theorem <ref> every time we call Eq. (<ref>).
The amplification protocol $\ms{U}_\mr{amp2}^{l_\imax}(t)$ employs $3$ times queries to $\ms{U}_\mr{amp1}^{l_\imax}(t)$.
Reflecting that the operations $\ms{U}_\mr{ini}^{l_\imax}$ and $\ms{R}$ require relatively a little resource (at-most $\order{\log l_\imax}$ elementary gates and complexity), the resource for $\ms{U}_\mr{amp2}^{l_\imax}(t)$ has the same scaling as the one for implementing the time evolution operators $e^{-it \sum_l l\omega \ket{l} \bra{l}}$ and $e^{-i \ms{H}_\eff^{4l_\imax}t}$.

As well as the amplification by symmetry, relying only on the oblivious amplitude amplification fails to efficiently enhances the success probability to obtain $\ket{0} \otimes \ket{\psi(t)}$ from $\order{l_\imax^{-1}}$ to $1-\order{\varepsilon}$.
When we do not use the first amplification, the protocol of the oblivious amplitude amplification is given by

\begin{equation}
\tilde{\ms{U}}_{\mr{amp2}, p}^{l_\imax} = - \left\{ \ms{U}_\mr{orig}^{l_\imax} \ms{R} (\ms{U}_\mr{orig}^{l_\imax})^\dagger \ms{R} \right\}^p \ms{U}_\mr{orig}^{l_\imax},
\end{equation}

where the operation $\ms{U}_\mr{orig}^{l_\imax}$ represents the time evolution without the amplification by symmetry, defined by

\begin{equation}
\ms{U}_\mr{orig}^{l_\imax}(t) = (\ms{U}_\mr{ini}^{l_\imax})^\dagger e^{-it \sum_l l\omega \ket{l} \bra{l}} e^{-i \ms{H}_\eff^{l_\imax}t}.
\end{equation}

Applying $\tilde{\ms{U}}_{\mr{amp2}, p}^{l_\imax}$ to the initial state $\ket{0} \otimes \ket{\psi(t)}$ returns $\ket{0} \otimes \ket{\psi(t)}$ whose amplitude increases from $1/\sqrt{2l_\imax}$ approximately in proportion to $p$.
The integer $p$ should be $\order{\sqrt{l_\imax}}$ to achieve the amplitude $1-\order{\varepsilon}$, reminiscent of the Grover's search algorithm.
In other words, $\order{\sqrt{l_\imax}}$ times call of $e^{-it \sum_l l\omega \ket{l} \bra{l}}$ and $e^{-i \ms{H}_\eff^{l_\imax}t}$ is required when we use $\tilde{\ms{U}}_{\mr{amp2}, p}^{l_\imax}$.
This is why we suggest the combination of the two amplification protocols, with which $\order{1}$ times call of $e^{-it \sum_l l\omega \ket{l} \bra{l}}$ and $e^{-i \ms{H}_\eff^{l_\imax}t}$ can amplify the amplitude of $\ket{0} \otimes \ket{\psi(t)}$ from $1/\sqrt{2l_\imax}$ to $1-\order{\varepsilon}$.

§ BLOCK-ENCODING OF EFFECTIVE FLOQUET HAMILTONIAN

In the previous section, we show that the combination of the two kinds of amplification protocols, implemented by $\ms{U}_\mr{amp2}^{l_\imax}(t)$ [See Eq. (<ref>)], provides the target time-evolved state $\ket{\psi(t)}$ with an arbitrarily small error $\order{\varepsilon}$.
The computational resource for $\ms{U}_\mr{amp2}^{l_\imax}(t)$ is mostly determined by that for $e^{-it \sum_l l\omega \ket{l} \bra{l}}$ and $e^{-i \ms{H}_\eff^{4l_\imax}t}$.
Our strategy is to implement these two time evolution operators by the qubitization technique <cit.>.
As discussed in Section <ref>, the computational resources are determined by how to introduce the block-encoding of the two static Hamiltonians, $\ms{H}_\mr{LP}^{l_\imax}=\sum_{l \in D^{l_\imax}} l \omega \ket{l}\bra{l} \otimes I$ (linear potential Hamiltonian) and $\ms{H}_\eff^{l_\imax}$ (effective Hamiltonian).
The aim of this section is to obtain an efficient block-encoding of them and to evaluate the costs required to implement $e^{-it \sum_l l\omega \ket{l} \bra{l}}$ and $e^{-i \ms{H}_\eff^{4l_\imax}t}$.

§.§ Block-encoding of linear potential Hamiltonian

We compose block-encoding of the linear potential Hamiltonian,

\begin{equation}
\ms{H}_\mr{LP}^{4l_\imax} = \sum_{l \in D^{4l_\imax}} l\omega \ket{l} \bra{l} \otimes I.
\end{equation}

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

\begin{eqnarray}
\ms{H}_\mr{LP}^{4l_\imax} &=& \sum_{l \in D^{4l_\imax}} \frac{\omega}{2} V_l^{4l_\imax} \otimes I, \label{Eq6A:LCU_linear_potential} \\
V_l^{4l_\imax} &=& \sum_{l^\prime = l}^{4l_\imax} \ket{l^\prime} \bra{l^\prime} - \sum_{l^\prime = -4 l_\imax + 1}^{l-1} \ket{l^\prime} \bra{l^\prime}. \label{Eq6A:Unitaries_linear_potential}
\end{eqnarray}

Using the block-encoding formalism of LCUs by Eq. (<ref>), we immediately obtain the oracle unitary gate and the oracle state by

\begin{eqnarray}
&& \ms{O}_\mr{LP}^{4l_\imax} = \sum_{l \in D^{4l_\imax}} (\ket{l}\bra{l})_b \otimes V_l^{4l_\imax} \otimes I, \label{Eq6A:Oracle_gate_linear_potential} \\
&& \braket{a^{4l_\imax}|\ms{O}_\mr{LP}^{4l_\imax}|a^{4l_\imax}}_b = \frac{\ms{H}_\mr{LP}^{4l_\imax}}{8l_\imax \omega }. \label{Eq6A:Oracle_state_linear_potential}
\end{eqnarray}

The subscript $b$ for the states $\ket{l}_b$ and $\ket{a^{4l_\imax}}_b$ represents a new ancillary system introduced for the block-encoding, requiring the number of qubits $n_b \in \order{\log l_\imax}$.
Combining Eqs. (<ref>) and (<ref>), the oracle unitary gate is rewritten by

\begin{equation}
\ms{O}_\mr{LP}^{4l_\imax} = \left (\sum_{l,l^\prime ; l^\prime-l \geq 0} \ket{l,l^\prime}\bra{l,l^\prime} - \sum_{l,l^\prime ; l^\prime-l < 0} \ket{l,l^\prime}\bra{l,l^\prime} \right) \otimes I .
\end{equation}

To implement this oracle, we use a comparator defined by

\begin{equation}\label{Eq6A:comparator}
\mr{Comp} \ket{l,l^\prime} \ket{0}_{b^\prime} = \begin{cases}
\ket{l,l^\prime} \ket{0}_{b^\prime} & \text{if $l^\prime \geq l$}, \\
\ket{l,l^\prime} \ket{1}_{b^\prime} & \text{if $l^\prime < l$},
\end{cases}
\end{equation}

with a single-qubit ancillary system $b^\prime$ <cit.>.
We can immediately confirm the relation,

\begin{equation}
\mr{Comp}^\dagger Z_{b^\prime} \mr{Comp} \ket{l,l^\prime} \ket{0}_{b^\prime} \ket{\psi} = \ms{O}_\mr{LP}^{4l_\imax}  \ket{l,l^\prime} \ket{0}_{b^\prime} \ket{\psi},
\end{equation}

for arbitrary inputs $\ket{l,l^\prime} \in \bbC^{8l_\imax}$ and $\ket{\psi} \in \mcl{H}$, where $Z_{b^\prime}$ denotes a Pauli $Z$ operator on the system $b^\prime$.
Since the comparator on $n$ qubits can be composed of $\order{n}$ elementary gates, the oracle unitary gate $\ms{O}_\mr{LP}^{4l_\imax}$ requires at-most $\order{\log l_\imax}$ gates.
The oracle state $\ket{a^{4l_\imax}}_b$, which has equal weights in $\ket{l}$, can be prepared by $\order{\log l_\imax}$ elementary gates.

We implement the time-evolution $e^{-i \ms{H}_\mr{LP}^{4l_\imax} t}$ with accuracy $1-\order{\varepsilon}$ by qubitization.
The number of additional gates, other than the block-encoding, is proportional to the number of ancillary qubits for the oracle $\ms{O}_\mr{LP}^{4l_\imax}$, that is, we need $\order{\log l_\imax}$ additional gates per query.
Based on Section <ref>, it can be executed with the following resources;

* Number of ancillary qubits; $\order{\log l_\imax}$

* Scaling of overall complexity;

\begin{equation}\label{Eq6A:gate_count_linear_potential}
\left\{ l_\imax \omega t + \frac{\log (1/\varepsilon)}{\log ( e + (l_\imax \omega t)^{-1} \log (1/\varepsilon) )} \right\} \log l_\imax.
\end{equation}

§.§ Block-encoding of effective Hamiltonian

We hereby provide a way to efficiently implement $e^{-i \ms{H}_\eff^{4l_\imax}t}$ by the qubitization technique.
We find out that naive block-encoding of the effective Hamiltonian $\ms{H}_\eff^{4l_\imax}$, given by Eq. (<ref>), faces at the severe increase of the complexity in $l_\imax$.
This originates from the fact that the number of nontrivial terms in $\ms{H}_\eff^{4l_\imax}$ increases linearly in $l_\imax$.
To avoid this problem, we take an alternative approach which yields only $\order{\log l_\imax}$ complexity as follows;

* Find a refined effective Hamiltonian $\ms{H}$ on the truncated Floquet-Hilbert space $\bbC^{8l_\imax} \otimes \mcl{H}$, such that $e^{-i \ms{H} t}$ accurately reproduces $e^{-i \ms{H}_\eff^{4l_\imax}t}$ with an arbitrarily small error $\order{\varepsilon}$.

* Construct block-encoding of $\ms{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 block-encoding of each Fourier component $H_m$ ($|m| \leq m_\imax$), given by

\begin{equation}\label{Eq6B:block_encoding_H_m}
\bra{G_m}_a O_m \ket{G_m}_a  = \frac{H_m}{\alpha_m}, \quad \alpha_m > 0,
\end{equation}

with an $n_a$-qubit ancillary system.
Each oracle unitary gate on $\bbC^{2^{n_a}} \otimes \mcl{H}$ is represented by $O_m$, and each oracle state $\ket{G_m}_a \in \bbC^{2^{n_a}}$ is generated by a quantum circuit $G_m$ as $\ket{G_m}_a = G_m \ket{0}^{\otimes n_a}$.
The second oracle is a query to coefficients of Fourier components, represented by $\alpha_m$.
We define the oracle unitary circuit $G_\mr{coef}$, acting on $\bbC^{2m_\imax+1}$ [i.e. $\order{1}$-qubit system], by

\begin{equation}\label{Eq6B:oracl_coefficient}
G_\mr{coef} \ket{0} = \sum_{|m| \leq m_\imax} \sqrt{\frac{\alpha_m}{\alpha}} \ket{m}, \quad \alpha = \sum_{|m| \leq m_\imax} \alpha_m.
\end{equation}

The parameter $\alpha$, defined by $\alpha = \sum_{|m| \leq m_\imax} \alpha_m$, characterizes the energy scale of the Hamiltonian $H(t)$.
In fact, considering $\norm{H_m} \leq \alpha_m$ for each $m$, it is bounded by $\norm{H(t)} \leq \sum_{|m| \leq m_\imax} \norm{H_m} \leq \alpha$.
The parameter $\alpha$ also provides the bound on the energy scale of time-dependent terms $\gamma$, as

\begin{equation}\label{Eq6B:gamma_alpha}
\gamma \leq \sum_{m \neq 0; |m| \leq m_\imax} \norm{H_m} \leq \alpha.
\end{equation}

Implementation of the oracle $G_\mr{coef}$ is equivalent to embedding probability distributions in quantum states <cit.>, which also appears as the oracle state for the qubitization technique in the case of a LCU [See Eq. (<ref>)].
Stimulated by its various applications covering linear equation solver <cit.> and quantum singular value transformation <cit.>, there have been various efficient implementations <cit.>.
Here, we suppose that the number of elementary gates required for the oracles $O_m$, $G_m$, and $G_\mr{coef}$ is at-most $C$.
As discussed later in Section <ref>, the form of time-periodic Hamiltonians $H(t)$ designated by Eq. (<ref>) is reasonable in that they involve various important classes such as LCUs.

§.§.§ Construction of a refined effective Hamiltonian

Let us compose a refined effective Hamiltonian, which accurately reproduces the dynamics of $e^{-i\ms{H}_\eff^{4l_\imax}t}$.
To come to the point, such a Hamiltonian equipped with suitability for efficient block-encoding is given by

\begin{equation}\label{Eq6B:def_effective_H_pbc}
\ms{H}_{\eff,\mr{pbc}}^{4l_\imax} = \ms{H}_\eff^{4l_\imax} + \sum_{(l,m) \in \partial F^{4l_\imax}} (\ket{l}\bra{l \oplus m} \otimes H_{-m} + \text{h.c.}),
\end{equation}

with $\partial F^{4 l_\imax} = \{ (l,m) \, | \, 4l_\imax-m_\imax+1 \leq l \leq 4l_\imax, \, 4l_\imax - l +1 \leq m \leq m_\imax \}$.
The integer $l \oplus m \in D^{4l_\imax}$ is defined modulo $8l_\imax$.
The additional terms in $\ms{H}_{\eff,\mr{pbc}}^{4l_\imax}$ connects the boundary of $D^{4l_\imax}$ so that the hopping terms, $\ket{l}\bra{l+m}$, induced by $H_{-m}$ can be translation-invariant under the shift of $\ket{l}$.
The subscript “pbc" comes from the periodic boundary conditions (PBC) for the hopping terms.

In our algorithm, we employ the time evolution $e^{-i \ms{H}_{\eff,\mr{pbc}}^{4l_\imax} t}$ instead of the original one $e^{-i \ms{H}_{\eff}^{4l_\imax} t}$, which is derived by Floquet theory.
In other words, we organize the amplification protocols with the refined effective Hamiltonian $\ms{H}_{\eff,\mr{pbc}}^{4l_\imax}$ as

\begin{eqnarray}
\ms{U}_\mr{amp1,pbc}^{l_\imax} (t) &=& (\ms{U}_\mr{ini}^{4l_\imax})^\dagger e^{-i \ms{H}_\mr{LP}^{4l_\imax}t} e^{-i \ms{H}_\mr{eff,pbc}^{4l_\imax}t} \ms{U}_\mr{ini}^{l_\imax}, \nonumber \\
&& \label{Eq6B:def_amp1_pbc} \\
\ms{U}_\mr{amp2,pbc}^{l_\imax} (t) &=& - \ms{U}_\mr{amp1,pbc}^{l_\imax}\ms{R} (\ms{U}_\mr{amp1,pbc}^{l_\imax})^\dagger \ms{R} \ms{U}_\mr{amp1,pbc}^{l_\imax}, \nonumber \\
&& \label{Eq6B:def_amp2_pbc}
\end{eqnarray}

and apply them to the initial state $\ket{0}\ket{\psi(0)}$.
The refined effective Hamiltonian is valid in a sense that it can also provide the exact time-evolved state as

\begin{equation}\label{Eq6B:approximate_amp2_pbc_result}
\ms{U}_\mr{amp2,pbc}^{l_\imax} (t) \ket{0}\ket{\psi(0)} = \ket{0}\ket{\psi(t)} + \order{\varepsilon}.
\end{equation}

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

(Refined effective Hamiltonian)

We choose $l_\imax \in \order{\gamma t + \log (1/\varepsilon)/\log \log (1/\varepsilon)}$ by Eq. (<ref>).
Then, the refined effective Hamiltonian $\ms{H}_{\eff,\mr{pbc}}^{4l_\imax}$ well reproduces the dynamics by the original one $\ms{H}_\eff^{4l_\imax}$ as

\begin{eqnarray}
\norm{\bra{0} \ms{U}_\mr{amp1,pbc}^{l_\imax} (t) \ket{0} \ket{\psi(0)} - \frac12 \ket{\psi(t)}} &\leq& \frac{\varepsilon}{3}, \label{Eq6B:thm_amp1_pbc}\\
\norm{ \ms{U}_\mr{amp2,pbc}^{l_\imax} (t) \ket{0} \ket{\psi(0)} - \ket{0}\ket{\psi(t)}} &\leq& \varepsilon, \label{Eq6B:thm_amp2_pbc}
\end{eqnarray}

for an arbitrary initial state $\ket{\psi(0)} \in \mcl{H}$.

We deliver its detailed proof in Appendix <ref>, and instead we briefly explain why the refined effective Hamiltonian is valid.
The proof relies mainly on the Lieb-Robinson bound, stated by Theorem <ref>.
Focusing on the amplification protocols $\ms{U}_\mr{amp1}^{l_\imax} (t)$ and $\ms{U}_\mr{amp1,pbc}^{l_\imax} (t)$, the resulting difference caused by them originates from their actions on the uniform initial state $\ket{a^{l_\imax}}\ket{\psi(0)}$, as Eq. (<ref>).
As we can see from Fig. <ref> (b), the Lieb-Robinson bound dictates that the dynamics during $\ms{U}_\mr{amp1}^{l_\imax} (t)$ is almost closed within $- 2l_\imax \lesssim l \lesssim 2l_\imax$ in the Fourier indices.
In contrast, the support of the additional terms of the refined effective Hamiltonian in Eq. (<ref>) is located at $\{ \ket{l} \}_l$ for $l \simeq \pm 4l_\imax$, which is sufficiently far from the Fourier indices relevant for the dynamics.
As a result, the refined protocol $\ms{U}_\mr{amp1,pbc}^{l_\imax} (t)$ transforms the state $\ket{a^{l_\imax}}\ket{\psi(0)}$ in almost the same way as the original one $\ms{U}_\mr{amp1}^{l_\imax} (t)$, and hence it also outputs the target state $\ket{\psi(t)}$ as Eq. (<ref>).
Discussion similar to Section <ref> ensures that the oblivious amplitude amplification under the refined effective Hamiltonian, represented by $\ms{U}_\mr{amp2,pbc}^{l_\imax} (t)$, provides $\ket{\psi(t)}$ with a sufficiently small error $\order{\varepsilon}$ as Eq. (<ref>).

§.§.§ Block-encoding of a refined effective Hamiltonian

The benefit of employing the refined effective Hamiltonian $\ms{H}_{\eff,\mr{pbc}}^{4l_\imax}$ instead of the original one $\ms{H}_{\eff}^{4l_\imax}$ is reduction of resources for implementing block-encoding.
We construct an oracle unitary gate $\ms{O}_\eff^{4l_\imax}$ and an oracle state $\ket{G_\eff^{4l_\imax}}$ for it.
With the help of the additional terms in Eq. (<ref>), the Hamiltonian is written in the following form;

\begin{eqnarray}
\ms{H}_{\eff,\mr{pbc}}^{4l_\imax} &=& \sum_{|m| \leq m_\imax} \mr{Add}_m^{4l_\imax} \otimes H_m - \ms{H}_\mr{LP}^{4l_\imax}, \\
\mr{Add}_m^{4l_\imax} &=& \sum_{l \in D^{4l_\imax}} \ket{l \oplus m}\bra{l}. \label{Eq6B:def_adder}
\end{eqnarray}

The unitary gate $\mr{Add}_m^{4l_\imax}$ is a full quantum adder that translates an index $l$ by $m$ modulo $8l_\imax$.
It can be implemented by $\order{\log l_\imax}$ 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 $n_a$-qubit system prepared for the block-encoding of $H_m$ [See Eq. (<ref>)], and the system $b$ is the one for $\ms{H}_\mr{LP}^{4l_\imax}$ having $8l_\imax$ dimension [See Section <ref>].
The other systems $c$ and $d$ have the dimension $2m_\imax +1$ and $2$ respectively, yielding the number of qubits $n_{cd} \in \order{\log m_\imax}$ (it is a constant here).
We define the oracle unitary gate $\ms{O}_\eff^{4l_\imax}$ by

\begin{eqnarray}
&& \ms{O}_\eff^{4l_\imax} = \nonumber \\
&\quad& \quad \ket{0}\bra{0}_d \otimes \sum_{|m| \leq m_\imax} \ket{m}\bra{m}_c \otimes I_b \otimes \mr{Add}_m^{4l_\imax}\otimes O_m \nonumber \\
&\quad& \qquad  - \ket{1}\bra{1}_d \otimes I_c \otimes I_a \otimes \ms{O}_\mr{LP}^{4l_\imax}.  \label{Eq6B:Oracle_gate_effective_H}
\end{eqnarray}

We also provide the oracle state $\ket{G_\eff^{4l_\imax}}$ by

\begin{eqnarray}
\ket{G_\eff^{4l_\imax}} &=& \sum_{m} \sqrt{\frac{\alpha_m}{\alpha+8l_\imax \omega}} \ket{0}_d \ket{m}_c \ket{a^{4l_\imax}}_b \ket{G_m}_a \nonumber  \\
&\qquad& \quad + \sqrt{\frac{8l_\imax \omega}{\alpha+8l_\imax \omega}} \ket{1}_d \ket{0}_c \ket{a^{4l_\imax}}_b \ket{0}_a^{\otimes n_a}. \nonumber \\
&& \label{Eq6B:Oracle_state_effective_H}
\end{eqnarray}

In the above formula, we implicitly include ancillary qubits to efficiently implement the oracles [e.g. the comparator for $\ms{O}_\mr{LP}^{4l_\imax}$ in Eq. (<ref>)].
The number of such ancillary qubits is smaller than $\order{\log l_\imax}$, and hence we neglect it below.
We can confirm that they provide the block-encoding of $\ms{H}_\eff^{4l_\imax}$ as follows;

\begin{equation}\label{Eq6B:block_encode_H_eff_pbc}
\braket{G_\eff^{4l_\imax}|\ms{O}_\eff^{4l_\imax}|G_\eff^{4l_\imax}} = \frac{\ms{H}_{\eff,\mr{pbc}}^{4l_\imax}}{\alpha + 8 l_\imax \omega},
\end{equation}

where we use the block-encoding for $H_m$ and $\ms{H}_\mr{LP}^{4l_\imax}$ represented by Eqs. (<ref>) and (<ref>).

We evaluate the resources for implementing the oracles.
The oracle unitary gate $\ms{O}_\eff^{4l_\imax}$ given by Eq. (<ref>) requires elementary gates at-most $\order{m_\imax (\log l_\imax + C) + \log l_\imax}$, since the unitary operators $\mr{Add}_m^{4l_\imax}$, $O_m$, and $\ms{O}_\mr{LP}^{4l_\imax}$ respectively yield $\order{\log l_\imax}$, $C$, $\order{\log l_\imax}$ gates.
On the other hand, the oracle state $\ket{G_\eff^{4l_\imax}}$ can be prepared by

\begin{equation}
\ket{G_\eff^{4l_\imax}} = G_\eff^{4l_\imax} \ket{w}_d (G_\mr{coef}\ket{0}_c) \ket{a^{4l_\imax}}_b \ket{0}_a^{\otimes n_a}, \\
\end{equation}

with using the quantum circuit $G_\eff^{4l_\imax}$ and the state $\ket{w}_d$ defined by

\begin{eqnarray}
G_\eff^{4l_\imax}  &=& I_d \otimes \sum_{|m| \leq m_\imax} \ket{m}\bra{m}_c \otimes I_b \otimes (G_m)_a, \\
\ket{w}_d &=& \sqrt{\frac{\alpha}{\alpha + 8l_\imax \omega }} \ket{0}_d + \sqrt{\frac{8l_\imax \omega}{\alpha+8l_\imax \omega}} \ket{1}_d. \nonumber \\
\end{eqnarray}

The quantum circuit $G_\eff^{4l_\imax}$ can be composed of $\order{m_\imax C}$ elementary gates.
The cost of the state preparation $\ket{w}_d$ is negligible compared to others since it requires only a single-qubit rotation.
The cost for preparing $\ket{a^{4l_\imax}}_b$ is at-most $\order{\log l_\imax}$ complexity and elementary gates.
To summarize, the cost for the oracles $\ms{O}_\eff^{4l_\imax}$ and $\ket{G_\eff^{4l_\imax}}$, which embody the refined effective Hamiltonian $\ms{H}_{\eff,\mr{pbc}}^{4l_\imax}$ via Eq. (<ref>), is at-most $\order{m_\imax (C+\log l_\imax)}$ elementary gates.
The additional gates per query, other than the block-encoding, amounts to $\order{n_a + \log (l_\imax)}$.

As a result, preparing a unitary circuit corresponding to $e^{-i \ms{H}_\eff^{4l_\imax} t}$ with an allowable error up to $\order{\varepsilon}$ is summarized as follows;

* Number of ancillary qubits;

\begin{equation}\label{Eq6B:number_ancillary}
n_a (l_\imax) = n_a + \order{\log l_\imax + \log m_\imax}.
\end{equation}

* Scaling of query complexity $q(l_\imax)$;

\begin{equation}
(\alpha + l_\imax \omega) t + \frac{\log (1/\varepsilon)}{\log (e + \{ (\alpha + l_\imax \omega) t \}^{-1} \log (1/\varepsilon))}.
\end{equation}

* Number of overall gates;

\begin{equation}
\order{q(l_\imax) m_\imax \{ C + n_a + \log l_\imax \} }
\label{Eq6B:gate_count_effective_H}
\end{equation}

Since $m_\imax$ is supposed to be a $\order{1}$ constant, we will omit it from the cost in the rest of the paper.

§ ALGORITHM AND ITS COMPUTATIONAL COST

This section provides the main result of this paper; we compose the efficient quantum algorithm for simulating the time-evolved state $\ket{\psi(t)}=U(t)\ket{\psi(0)}$ under the Hamiltonian $H(t+T)=H(t)$.
We note that we take two different approaches depending on the time scale of interest.
The first case is the adiabatic-like case, in which we are interested in $\order{1}$-period dynamics with $t/T \in \order{1}$.
We call it “adiabatic-like" since long-time dynamics during $0 \leq t \leq 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 long-time case, where we are interested in multiple-period dynamics at $t/T \in \Omega (1)$.
In that case, the the period $T$ is not so large and we often consider long-time dynamics at $t \gg T$, exemplified by laser-irradiated materials.

Ancillary qubits     Query complexity
(for $H$ and $H(t)$)      0     $\alpha t \cdot  (\alpha t/\varepsilon)^{1/p}$,   ($p$-th order)
Time-independent $H$
(Qubitization <cit.>)      $n_a+\order{1}$ (Independent of $t$, $\varepsilon$)    $ \alpha t + \frac{\log (1/\varepsilon)}{\log (e + (\alpha t)^{-1} \log (1/\varepsilon) )}  $
Time-periodic $H(t)$
(Adiabatic, Theorem <ref>)      $n_a + \order{\log (\gamma t) + \log \log (1/\varepsilon)}$     $\alpha t + \frac{\log (1/\varepsilon)}{\log (e + \{ \alpha t + o(\log (1/\varepsilon))\}^{-1} \log (1/\varepsilon))}$
Time-periodic $H(t)$
(Generic, Theorem <ref>)      $\quad n_a + \order{\log (\gamma/\omega) + \log \log (\omega t/\varepsilon)} \quad $     $\quad \alpha t + \frac{ \omega t \log (\omega t/\varepsilon)}{\log ( e + \{ \alpha/\omega + o(\log (\omega t/\varepsilon)) \}^{-1} \log (\omega t/\varepsilon))} \quad $

Time-dependent $H(t)$
(Dyson series <cit.>)      $n_a + \order{\log \{ (\gamma \omega t / \alpha + \alpha t) / \varepsilon \} }$     $\alpha t \frac{\log(\alpha t/\varepsilon)}{\log \log (\alpha t /\varepsilon)}$

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 them by those which have the similar scales (e.g. $\norm{H(t)} \leftrightarrow \alpha$, $\norm{\dv{t}H(t)} \leftrightarrow \omega \gamma$). We also note that the query complexity is measured by the oracles, but their definitions depends on the algorithms. See Eq. (<ref>) [for Trotterization], Eq. (<ref>) [for qubitization], Eq. (<ref>) and Eq. (<ref>) [for our algorithm], and Eq. (<ref>) [for the truncated Dyson-series algorithm].

§.§ Adiabatic-like cases

We consider the adiabatic-like cases, where long-time dynamics over $\order{1}$-periods is of interest.
For simplicity, we first consider the dynamics within one period at $t \in [0,T]$.
The algorithm in this case is composed of the following steps;

* Determine the truncation order of the Fourier index $l_\imax$ by Eq. (<ref>), for the given time $t$.

* Compose two unitary gates $\ms{U}_\mr{LP}^{4l_\imax}(t)$ and $\ms{U}_\eff^{4l_\imax}(t)$, corresponding to the time evolution operators $e^{-i \ms{H}_\mr{LP}t}$ and $e^{-i \ms{H}_{\eff,\mr{pbc}}^{4l_\imax} t}$ by the qubitization technique; We can find such operators satisfying

\begin{eqnarray}
&& \ms{U}_\mr{LP}^{4l_\imax} (t) \ket{0}^{\otimes n_a (4l_\imax)} \ket{l} \ket{\psi} \nonumber \\
&& \quad = \ket{0}^{\otimes n_a (4l_\imax)} e^{-il\omega t}\ket{l} \ket{\psi} + \order{\varepsilon}, \\
&& \ms{U}_\mr{eff}^{4l_\imax} (t) \ket{0}^{ n_a (4l_\imax)} \ket{l} \ket{\psi}
\nonumber \\
&&  \quad = \ket{0}^{\otimes  n_a (4l_\imax) }e^{-i \ms{H}_{\eff,\mr{pbc}}^{4l_\imax}t} \ket{l}\ket{\psi} + \order{\varepsilon},
\end{eqnarray}

for arbitrary inputs $l \in D^{4l_\imax}$ and $\ket{\psi} \in \mcl{H}$.

* Construct the quantum circuit for the amplification protocols with embedding $\ms{U}_\mr{LP}^{4l_\imax}(t)$ and $\ms{U}_\eff^{4l_\imax}(t)$ as $e^{-i \ms{H}_\mr{LP}t}$ and $e^{-i \ms{H}_{\eff,\mr{pbc}}^{4l_\imax} t}$ [See Section <ref>];

\begin{eqnarray}
\bar{\ms{U}}_\mr{amp1}^{l_\imax} &=& (\ms{U}_\mr{ini}^{4l_\imax})^\dagger \ms{U}_\mr{LP}^{4l_\imax} \ms{U}_\eff^{4l_\imax} \ms{U}_\mr{ini}^{l_\imax}, \label{Eq7A:amp1_adiabatic}\\
\bar{\ms{U}}_\mr{amp2}^{l_\imax} &=& - \bar{\ms{U}}_\mr{amp1}^{l_\imax} \ms{R} (\bar{\ms{U}}_\mr{amp1}^{l_\imax})^\dagger \ms{R} \ms{U}_\mr{amp1}^{l_\imax}. \label{Eq7A:amp2_adiabatic}
\end{eqnarray}

* Prepare the initial state $\ket{0}^{\otimes n_a(4l_\imax)} \ket{0} \ket{\psi(0)}$ and apply the unitary operation $\bar{\ms{U}}_\mr{amp2}^{l_\imax}(t)$ to this state, which results in

\begin{eqnarray}
&& \bar{\ms{U}}_\mr{amp2}^{l_\imax} (t) \ket{0}^{\otimes n_a (4l_\imax)} \ket{0}\ket{\psi(0)} \nonumber \\
&& \qquad  = \ket{0}^{\otimes n_a (4l_\imax)} \ket{0}\ket{\psi(t)} + \order{\varepsilon}.
\end{eqnarray}

Step 1. determines the dimension of the truncated Floquet-Hilbert space as $\bbC^{8l_\imax} \otimes \mcl{H}$ so that it can accurately reproduce $\ket{\psi(t)}$ with the precision $1-\order{\varepsilon}$, as discussed in Section <ref>.
In Step 2., we employ the qubitization technique to realize the time evolution in the truncated Floquet-Hilbert space, as discussed in Section <ref>.
We note that $\ms{U}_\mr{LP}^{4l_\imax} (t)$ and $\ms{U}_\mr{eff}^{4l_\imax} (t)$ can exploit a common $\order{\log (l_\imax)}$-qubit ancillary system.
It is sufficient to prepare an $\{ n_a(4l_\imax) \}$-qubit ancillary state $\ket{0}^{\otimes n_a(4l_\imax)}$ [See Eq. (<ref>)].
Steps 3. and 4. execute the amplification of the time-evolved state $\ket{\psi(t)}$ as discussed in Section <ref>.
By projecting the ancillary state to $\ket{0}^{\otimes n_a (4l_\imax)}\ket{0}$, we succeed in preparing the target state $\ket{\psi(t)}$ with precision $1-\order{\varepsilon}$, where the success probability is $1-\order{\varepsilon}$.
These steps do not change the scaling of the required resource from Eqs. (<ref>) and (<ref>), reflecting that the unitary operations $\ms{U}_\mr{ini}$ and $\ms{R}$ yield at-most $\order{\log l_\imax}$ gates.

We finally determine the resource for simulating time-periodic Hamiltonians in adiabatic-like regimes.
The ancillary system should involve the degree of freedom for Fourier indices $\{ \ket{l} \}_{l \in D^{4l_\imax}}$ and $n_a(4l_\imax)$ qubits for qubitization.
The total number of ancillary qubits is

\begin{equation}
n_a(4l_\imax) + \lceil \log_2 (8l_\imax) \rceil \in n_a + \order{\log (\gamma t) + \log \log (1/\varepsilon)}.
\end{equation}

The numbers of the oracles $O_m$ and $G_\mr{coef}$ are dominated by that for the qubitization, designated by Eq. (<ref>).
Considering the assumption $\omega t \in \order{1}$, the query complexity scales as

\begin{eqnarray}
&& (\alpha + \gamma) t + \frac{\log (1/\varepsilon)}{\log (e + (\gamma t)^{-1} \log (1/\varepsilon))} \nonumber \\
&& \quad + \frac{\log (1/\varepsilon)}{\log (e + \{(\alpha + \gamma) t + o(\log (1/\varepsilon))\}^{-1} \log (1/\varepsilon))}, \nonumber \\
&& \in \order{\alpha t + \frac{\log (1/\varepsilon)}{\log (e + \{\alpha t + o(\log (1/\varepsilon)) \}^{-1} \log (1/\varepsilon))} }, \nonumber \\
\end{eqnarray}

where the $o(\log (1/\varepsilon))$ term scales as Eq. (<ref>) originating from the scaling of $l_\imax$.
The scaling dependent on $\gamma$ can be absorbed into that of $\alpha$, due to the relation $\gamma \leq \alpha$.
We summarize the results in Theorem <ref> and Table <ref>.
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 $\order{\log l_\imax}$ gates.
By the $\order{1}$-times repetition of this procedure for $t \in [0,T]$, the same result for $\order{1}$-period dynamics is obtained.

§.§ Generic long-time cases

We next consider the other generic cases where we are interested in long-time dynamics over multiple periods as $\omega t \in \Omega (1)$.
In this regime, we take a different strategy from that for the adiabatic cases; We split the time $t$ by $t = (n+\delta) T$  with $n \in \bbN$ and $\delta \in [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(\delta T)$.
The algorithm is composed of the following steps;

* Split the time by $t=(n+\delta) T$ with $n \in \bbN$ and $\delta \in [0,1)$. Determine the truncation order $l_\imax^T$ by substituting $T$ and $\varepsilon/n$ into $t$ and $\varepsilon$ of Eq. (<ref>);

\begin{eqnarray}
l_\imax^T &\in& \Theta \left( \gamma T + \frac{\log (n/\varepsilon)}{\log (e + (\gamma T)^{-1} \log (n/\varepsilon))}\right) \nonumber \\
&=& \Theta \left( \gamma / \omega  + \frac{\log (\omega t/\varepsilon)}{\log (e + (\gamma/\omega)^{-1} \log (\omega t/\varepsilon) )}\right). \nonumber \\
\end{eqnarray}

* Compose a unitary gate $\bar{\ms{U}}^{l_\imax^T}_\mr{amp2}(T)$, which satisfies

\begin{eqnarray}
&& \bar{\ms{U}}_\mr{amp2}^{l_\imax^T} (T) \ket{0}^{\otimes n_a (4l_\imax^T)}\ket{0}\ket{\psi(0)} \nonumber \\
&& \qquad  = \ket{0}^{\otimes n_a (4l_\imax^T)} \ket{0}\ket{\psi(T)} + \order{\varepsilon/n},
\end{eqnarray}

by the qubitization and the amplitude amplification protocols (Follow Steps 2.-4. of Section <ref> wih substituting $T$, $\varepsilon/n$, and $l_\imax^T$ into $t$, $\varepsilon$, and $l_\imax$ respectively).

* Apply the unitary gate  $\bar{\ms{U}}^{l_\imax^T}_\mr{amp2}(T)$ to the initial state $n$ times, which results in

\begin{eqnarray}
&& [\bar{\ms{U}}_\mr{amp2}^{l_\imax^T} (T)]^n \ket{0}^{\otimes n_a (4l_\imax^T)} \ket{0}\ket{\psi(0)} \nonumber \\
&& \qquad  = \ket{0}^{\otimes n_a (4l_\imax^T)}\ket{0}\ket{\psi(nT)} + \order{\varepsilon}.
\end{eqnarray}

* Prepare the unitary gate $\bar{\ms{U}}^{l_\imax^T}_\mr{amp2}(\delta T)$ by substituting $\delta T$, $\varepsilon$, and $l_\imax^T$ into $t$, $\varepsilon$, and $l_\imax$ respectively in Steps 2.-4. of Section <ref>. Applying it once to the above state results in

\begin{eqnarray}
&& \bar{\ms{U}}^{l_\imax^T}_\mr{amp2}(\delta T) [\bar{\ms{U}}_\mr{amp2}^{l_\imax^T} (T)]^n \ket{0}^{\otimes n_a (4l_\imax^T)} \ket{0}\ket{\psi(0)} \nonumber \\
&& \qquad = \ket{0}^{\otimes n_a (4l_\imax^T)} \ket{0}\ket{\psi(t)} + \order{\varepsilon}.
\end{eqnarray}

for arbitrary initial states $\ket{\psi(0)}$.

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 $\order{\varepsilon /n}$ so that we can obtain the time-evolved state $\ket{\psi(nT)}$ with an error up to $\order{\varepsilon}$ after the $n$-times repetition in Step 3.
The cost for the unitary gate $\bar{\ms{U}}^{l_\imax^T}_\mr{amp2}(T)$ is dominated only by that of $\ms{U}_\eff^{4l_\imax^T}$, which reproduces the time evolution $\exp ( -i \ms{H}_{\eff,\mr{pbc}}^{4l_\imax^T} T)$ by the qubitization.
This comes from the fact that the time evolution $e^{-i \ms{H}_\mr{LP}^{4l_\imax^T} T}$ is trivial due to

\begin{equation}
e^{-i \ms{H}_\mr{LP}^{4l_\imax^T} T} = \sum_{l \in D^{4l_\imax^T}} e^{i l \omega T} \ket{l}\bra{l} \otimes I = 1,
\end{equation}

and we do not need $\ms{U}_\mr{LP}^{4l_\imax^T}$ in Eq. (<ref>).
As a result, the query complexity for $n$-times implementation of $\bar{\ms{U}}^{l_\imax^T}_\mr{amp2}(T)$ during Step 3. is  proportional to

\begin{eqnarray}
&& n \left\{ (\alpha + l_\imax^T \omega) T + \frac{\log (n/\varepsilon)}{\log ( e + \{ (\alpha + l_\imax^T \omega)T \}^{-1} \log (n/\varepsilon) )} \right\} \nonumber \\
&& \in \order{\alpha t + \frac{\omega t \log (\omega t /\varepsilon)}{\log (e+ \{ \alpha/\omega + o(\log (\omega t/\varepsilon))\}^{-1} \log (\omega t/\varepsilon))}}. \nonumber \\
\label{Eq7B:cost_amp2_long_time}
\end{eqnarray}

We use the fact $n \in \order{\omega t}$ in the above relation.
Step 4. realizes the remaining micromotion $U(\delta T)$ for the duration $\delta T$ with an error up to $\order{\varepsilon}$.
We remark that the choice of the truncation order $l_\imax^T$ so far is sufficient to achieve the precision $1-\order{\varepsilon}$ for $U(\delta T)$ due to $l_\imax^T > l_\imax^{\delta T}$.
This implies that we can reuse the ancillary state $\ket{0}^{\otimes n_a (4l_\imax^T)}$ for the qubitization approach to the time-evolution $U(\delta T)$ in Step 4.
The cost for implementing $\bar{\ms{U}}^{l_\imax^T}_\mr{amp2}(\delta T)$ once is obtained by setting $t=\delta T$ in adiabatic-like cases.
It is smaller than the cost for implementing the time-evolution $U(T)$, and does not affect the scaling of the computational 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 $n_a(4l_\imax)+\lceil \log_2 (8l_\imax^T) \rceil$, which is bounded by

\begin{equation}\label{Eq7B:ancillary_long_time}
n_a + \order{\log (\gamma /\omega) + \log \log (\omega t /\varepsilon)}.
\end{equation}

The query complexity is dominated by Eq. (<ref>).
We summarize them in Theorem <ref> and Table <ref>.

§.§ Comparison with other algorithms

Let us compare our algorithm on time-periodic Hamiltonians with other quantum algorithms for Hamiltonian simulation, based on Table <ref>.
We pick up the qubitization technique <cit.> for time-independent Hamiltonian $H$ and the truncated Dyson-series algorithm <cit.> for generic time-dependent Hamiltonian $H(t)$, whose resources have the best scaling in $t$ and $1/\varepsilon$ 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 be divided into $H(t)=\sum_{r=1}^\Gamma H_r(t)$, where every term in $H_r(t)$ commutes with one another at every time, the first-order Trotterization roughly approximates the time-evolution $U(t)$ by

\begin{equation}
U(t) = \prod_{l=0}^{M-1} \left[ \prod_{r=1}^\Gamma \left( \mcl{T} e^{-i \int_{t_l}^{t_{l+1}} d t^\prime H_r(t^\prime)} \right) \right] + \order{\frac{(\alpha t)^2}{M}},
\end{equation}

with $M$ partitions of the time $t$ as $t_l=l t/M$ <cit.>.
Reflecting that the Trotterization error polynomially decreases in the partition number $M$, we need queries to a layered quantum circuit

\begin{equation}\label{Eq6C:oracle_trotter}
\prod_{r=1}^\Gamma \left( \mcl{T} e^{-i \int_{t_l}^{t_{l+1}} d t^\prime H_r(t^\prime)} \right)
\end{equation}

for $\order{(\alpha t)^2/\varepsilon}$ times to achieve the allowable error $\varepsilon$.
While the coefficient $\alpha$ and the powers in $t$ and $1/\varepsilon$ can be improved by higher-order Trotterization <cit.>, its resource increases polynomially in $1/\varepsilon$.
Our algorithm has better scaling of elementary gates than Trotterization in that the resource increases logarithmically in $1/\varepsilon$.

The comparison with the qubitization and the truncated Dyson-series algorithm is instructive for evaluating the efficiency of our algorithm due to the inclusion relation,

\begin{eqnarray}
&& [\text{Set of time-independent $H$}] \nonumber \\
&& \quad \subset [\text{Set of time-periodic $H(t)=H(t+T)$}] \nonumber \\
&& \quad \subset [\text{Set of generic time-dependent $H(t)$}].
\end{eqnarray}

The qubitization technique achieves the least number of ancillary qubits in a sense that it is independent of $t$ and $1/\varepsilon$.
It also has the best query complexity, in which the optimal scaling both in $t$ and $1/\varepsilon$ appears in an additive way as Eq. (<ref>).
The number of ancillary qubits and the query complexity for the qubitization provide the best bound for time-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 with the time-dependent Schrödinger equation.
The cost reduction compared to them can be interpreted as improvement of efficiency due to the time-periodicity.

Before going to the comparison respectively for the adiabatic-like regime and the generic long-time regime, we remark several points in common.
First, we replace some parameters by those which have similar scales to simply compare these algorithms.
For instance, the parameters $\alpha$ and $\gamma$ respectively give the energy scales of the whole Hamiltonian $H(t)$ and the time-dependent terms $H(t)-H_0$ according to Eqs. (<ref>) and (<ref>).
We substitute $\alpha$ and $\omega \gamma$ respectively for $\sup_t (\norm{H(t)})$ and $\sup_t (\norm{\dv{t} 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. (<ref>).
In contrast, the truncated Dyson-series algorithm for generic cases employs the oracle

\begin{equation}\label{Eq6C:oracle_dyson}
\mr{Ham-}O = \sum_{l=1}^M \ket{l} \bra{l} \otimes O(t_l),
\end{equation}

where $l$ labels the discretized time $t_l$ with the partition number $M \in \order{1/\varepsilon}$ <cit.>.
The oracle $\mr{Ham-}O$ includes multiple implementation of $O(t_l)$, which is an oracle for a static instantaneous Hamiltonian $H(t_l)$ (some specific cases such as LCUs and sparse-access matrices with integrable time-dependency can be simplified <cit.>).
Our algorithm for time-periodic Hamiltonians uses the oracles $\{ O_m,\ket{G_m}\}$, which gives block-encoding of each Fourier component $H_m$ as Eq. (<ref>), and the oracle for the coefficients, $G_\mr{coeff}$ given by Eq. (<ref>).
The query complexity in our algorithm is roughly measured by the oracle for a static operator $H_m$ since the latter one $G_\mr{coeff}$, a quantum gate on at-most $\order{1}$ qubits, is usually negligible.
Therefore, note that our algorithm adopts essentially the same measure for the query complexity as the query complexity, while the truncated Dyson-series algorithm counts it by rather complicated oracles involving discretized time.

§.§.§ Adiabatic-like cases

We assess the cost for simulating time-periodic Hamiltonians in the adiabatic-like cases based on Table <ref>.
First, we compare the number of ancillary qubits with those for other algorithms.
What should be noted in our algorithm is its scaling in the inverse error $1/\varepsilon$.
Our algorithm requires $\order{\log \log (1/\varepsilon)}$ additional qubits, whose number lies just in the middle of that for the qubitization ($\varepsilon$-independent) and that for the truncated Dyson-series algorithm [$\order{\log (1/\varepsilon)}$].
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 $\order{\log (1/\varepsilon)}$ ancillary qubits via the oracle Eq. (<ref>).
In contrast, we introduce the truncation order $l_\imax$ for the Floquet-Hilbert space.
The error by this cutoff scales as $\order{(\gamma t /l_\imax)^{l_\imax}}$, whose decay is faster than the exponential function $e^{-\order{l_\imax}}$.
This leads to the reduction of ancillary qubits from $\order{\log (1/\varepsilon)}$ to $\order{\log \log (1/\varepsilon)}$.

Next, we compare the query complexity.
When we fix the inverse error $1/\varepsilon$ and increase time $t$, it scales as $\order{\alpha t}$.
The linear increase in time $t$ implies the optimality of our algorithm in time.
When we consider $1/\varepsilon \to \infty$ for given time, the query complexity scales as

\begin{equation}
\frac{\log (1/\varepsilon)}{\log \log \log (1/\varepsilon)},
\end{equation}

where we substitute the form of the $o(\log (1/\varepsilon))$ term, Eq. (<ref>), into Eq. (<ref>).
This scaling is nearly-optimal in $1/\varepsilon$.
Importantly, their contributions affect the query complexity in an additive way as $\alpha t + o (\log (1/\varepsilon))$.
This means that our algorithm deals with time-periodic systems with the cost sufficiently close to qubitization [See Eq. (<ref>)], and saves resource compared to the truncated Dyson-series algorithm requiring the productive query complexity $\alpha t \times o(\log (\alpha t /\varepsilon))$.
As long as we suppose polynomial accuracy $\varepsilon \in 1/\poly{N}$, implying $\log (1/\varepsilon) \lesssim \alpha t$, the query complexity can achieve the best scaling for time-independent systems.

§.§.§ Generic long-time cases

Here, we evaluate the computational resources for generic long-time regimes.
We suppose that the frequency $\omega$ is constant and hence small compared to the whole-system energy scales $\alpha$ and $\gamma$, which usually increase polynomially in the system size $N$.
This assumption comes from the fact that $\omega$ typically represents the frequency of external drives (e.g. light), which is size-independent.
More importantly, when the frequency $\omega$ is $\poly{N}$ so that it becomes comparable to $\alpha$ and $\gamma$, the dynamics can be efficiently simulated by the methods for time-independent Hamiltonians with the help of high-frequency expansions <cit.>, as we prove in Appendix <ref>.
Therefore, it is sufficient to consider the computational resources under $\omega \in \order{N^0}$ and $\alpha,\gamma \in \poly{N}$.

The number of ancillary qubits is similar to that for the adiabatic cases.
We note that its scaling in the time [$\order{\log \log t}$] overwhelms that of the truncated Dyson-series algorithm [$\order{\log t}$] in addition to the inverse error $1/\varepsilon$.
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/\varepsilon$.

We discuss the query complexity given by Theorem <ref>.
Under the fixed allowable error $\varepsilon$, the scaling of the query complexity in sufficiently large time $t$ is given by

\begin{equation}
\alpha t + \omega t \frac{\log (\omega t /\varepsilon)}{\log \log \log (\omega t /\varepsilon)} \leq \{\alpha + \log (1/\varepsilon) \} t + \omega t \log (\omega t).
\end{equation}

Due to the second term in the right hand side, the query complexity has at least nearly-optimal dependence in $t$, which involves logarithmic correction $\log (\omega t)$.
However, reflecting the assumptions $\alpha \in \poly{N}$ and $\omega \in \order{N^0}$, the $\order{\omega t \log (\omega t)}$ term can be non-negligible compared to $\order{\alpha t}$ only when the time reaches $t \sim e^{\alpha/\omega}/\omega \in \order{e^{\poly{N}} T}$.
Therefore, unless we consider unpractical time scales $t \in e^{\Theta (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/\varepsilon$, the query complexity scales as

\begin{equation}
\omega t \frac{\log (1/\varepsilon)}{\log \log \log (1/\varepsilon)}.
\end{equation}

Thus, our algorithm achieves nearly-optimal scaling in $1/\varepsilon$ with slight deviation from the optimal dependence by $\log \log (1/\varepsilon) / \log \log \log (1/\varepsilon)$.

When it comes to the combined scaling of the query complexity in $t$ and $1/\varepsilon$, we note that it has an additive form as

\begin{equation}
\alpha t + \omega t \times o(\log (1/\varepsilon)),
\end{equation}

where we neglect $\order{\omega t \log (\omega t)}$.
Although linear increase in $t$ couples with logarithmic increase in $1/\varepsilon$, the coefficient $\omega t \in \order{N^0} t$ is much smaller than the whole energy-time scale $\alpha t \in \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. (<ref>)] and saves much cost compared to the truncated Dyson-series algorithm, where the query complexity scales as $\poly{N} t \times o(\log (\alpha t/\varepsilon))$.

§ ILLUSTRATIVE EXAMPLES

In this section, we briefly discuss some potential applications of the algorithm.
We expect that 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 can be applied to quantum phase estimation for instance.
We suggest the simplest examples for both applications below.

§.§ 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 matter physics and quantum chemistry.
We pick up an $N$-site Fermi-Hubbard model under laser light as the simplest case;

\begin{eqnarray}
H(t) &=& H_\text{Hub} + H_\mr{ext}(t), \\
H_\text{Hub} &=& \sum_{k } \sum_{\sigma= \uparrow, \downarrow} \epsilon_{k} \hat{n}_{k\sigma} + U \sum_{x} \hat{n}_{x\uparrow} \hat{n}_{x\downarrow}. \label{Eq8:Fermi_Hubbard}
\end{eqnarray}

Here, $\hat{n}_{k\sigma}=\hat{c}_{k\sigma}^\dagger \hat{c}_{k\sigma}$ and $\hat{n}_{x\sigma}=\hat{c}_{x\sigma}^\dagger \hat{c}_{x\sigma}$ are number operators of electrons in the momentum and real spaces respectively, generated by fermionic annihilation operators $\hat{c}_{k\sigma}$ and $\hat{c}_{x\sigma}$.

The time-periodic term $H_\mr{ext}(t)$ represents the coupling with light.
When we shine linearly-polarized light with the frequency $\Omega$, it is given by

\begin{equation}
H_\mr{ext}(t) = \sin (\Omega t) \sum_{x, \sigma} V_x \hat{n}_{x\sigma},
\end{equation}

which results in $T=2\pi/\Omega$ and $m_\imax=1$.
The Fourier components $H_m$ are given by $H_0=H_\mr{Hub}$ and $H_{\pm 1} = (\pm i) \sum_{x,\sigma} (V_x/2) \hat{n}_{x\sigma}$.
To evaluate the cost of simulating $\ket{\psi(t)}$, we compose the oracles for them.
We employ a unitary operation, called fermionic fast Fourier transform (FFFT) <cit.>, which transforms the basis in the momentum space $k_i$ to that in the real space $x_i$ as $\hat{n}_{k_i \sigma} = \mr{FFFT}^\dagger \hat{n}_{x_i \sigma} \mr{FFFT}$.
We map the fermionic system to a spin system by Jordan-Wigner transformation as $\hat{n}_{x\sigma}=(1+Z_{x\sigma})/2$ (“$x\sigma$" denotes an index for qubits).
By neglecting constant terms and a conserved particle number $\sum_{x\sigma} \hat{n}_{x\sigma}$, the block-encoding for $\{H_m\}$ can be constructed by the technique for LCUs, as Eq. (<ref>).
Assuming $\epsilon_{k}, U, V_x \geq 0$ without loss of generality, this leads to

\begin{eqnarray}
O_0 &=& \ket{0}\bra{0} \otimes \sum_{x, \sigma} \ket{x,\sigma}\bra{x,\sigma} \otimes \mr{FFFT}^\dagger Z_{x\sigma} \mr{FFFT} \nonumber \\
&& \quad + \ket{1}\bra{1} \otimes \sum_{x,\sigma} \ket{x,\sigma}\bra{x,\sigma} \otimes Z_{x\uparrow} Z_{x\downarrow}, \label{Eq8:Oracle_gate_Fermi_Hubbard} \\
\ket{G_0} &=& \sum_{x,\sigma} \left( \sqrt{\frac{\varepsilon_{x}}{\alpha_0}}\ket{0} +  \sqrt{\frac{U}{2\alpha_0}}\ket{1} \right) \ket{x,\sigma}, \label{Eq8:Oracle_state_Fermi_Hubbard} \\
O_{\pm 1} &=& I \otimes \sum_{x,\sigma} \ket{x,\sigma}\bra{x,\sigma} \otimes (\pm i Z_{x\sigma}), \label{Eq8:Oracle_gate_external}  \\
\ket{G_{\pm 1}} &=& \ket{0} \sum_{x,\sigma} \sqrt{\frac{V_x}{ 2 \alpha_{\pm 1}}} \ket{x,\sigma}, \label{Eq8:Oracle_state_external}
\end{eqnarray}

with $\alpha_0 = 2 \sum_{x} \epsilon_{x} + UN$ and $\alpha_{\pm 1}= \sum_{x} V_x$.
In the above oracles, the unitary gate $O_0$ requires much cost due to $\order{1}$-times usage of FFFT, which can be implemented with at-most $\order{N \log N}$-depth quantum circuits composed of adjacent two-qubit gates <cit.>.
The depth needed for each oracle is at-most $\order{N \log N}$.
The number of ancillary qubits for them to express $\{ \ket{0}\ket{x,\sigma}, \ket{1}\ket{x,\sigma} \}$ amounts to $n_a \in \order{\log N}$.
The energy scales of the whole Hamiltonian, $\alpha$, and that of time-dependent terms, $\gamma$, are respectively described as follows [See Eqs. (<ref>) and (<ref>)];

\begin{equation}
\alpha =  2 \sum_{x} \epsilon_{x} + UN + 2\sum_{x} V_x, \quad \gamma \leq \sum_{x} V_x.
\end{equation}

When we define the characteristic local energy scale by $\alpha_\mr{loc}=\max (\epsilon_k,U,V_x)$, they are approximately described by $\alpha, \gamma \in \order{\alpha_\mr{loc} N}$.

Finally, if we are interested in the time-evolved state $\ket{\psi(t)}$ over multiple periods, the following resources are required to achieve the precision $1-\order{\varepsilon}$;

* Number of ancillary qubits;

\begin{equation}\label{Eq8A:ancillary_fermi_hubbard}
\order{\log ( \alpha_\mr{loc} N/\Omega ) + \log \log (\Omega t /\varepsilon)}.
\end{equation}

* Overall gate complexity;

\begin{equation}\label{Eq8A:complexity_fermi_hubbard}
\order{\left\{ \alpha_\mr{loc} Nt + \Omega t \log (\Omega t/\varepsilon) \right\} N \log N}.
\end{equation}

The above results are based on Eqs. (<ref>) and (<ref>), where we neglect $\order{\log (\gamma /\omega)+\log \log (\omega t /\varepsilon)}$ additional gates for each query to the oracles.
According to classical numerical calculations <cit.>, the above model is expected to host high-harmonic generation, where intense oscillation with the frequency $n\Omega$ ($n=2,3,\hdots$) arises in response to laser light with the frequency $\Omega$.
Our algorithm allows to efficiently identify such a nontrivial nonequilibrium phenomenon with guaranteed accuracy by simulating the Fourier spectrum of some observables $\braket{\psi(t)|O|\psi(t)}$ (e.g. electric current) \cite{hhg}.

Another interesting example is a discrete time crystal as a phase of matter inherent in nonequilibrium, where time-translation symmetry is spontaneously broken \cite{Khemani2016-dq,Else2016-mg,Khemani2019-pf}.
When we choose $H_ext(t)$ by uniform circularly-polarized ac field represented by
\begin{equation}
H_\mr{ext}(t) = V \sum_{x} \left( e^{i\Omega t} \hat{c}_{x\uparrow}^\dagger \hat{c}_{x\downarrow} + e^{-i \Omega t} \hat{c}_{x\downarrow}^\dagger \hat{c}_{x\uparrow} \right),
\end{equation}
the Fermi-Hubbard model becomes a potential platform for a time-crystalline phenomenon protected by Floquet dynamical symmetry \cite{Chinzei2020dtc}.
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. (\ref{Eq8A:ancillary_fermi_hubbard}) and (\ref{Eq8A:complexity_fermi_hubbard}).

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 exponentially-decaying Fourier components $H_m ≲e^-|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, other than the above examples (See also Appendix \ref{Asubsec3:example_expo}).

\subsection{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 \cite{Albash2018RevModPhys}, 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 \cite{Aspuru-Guzik2005-nr,Du2010-hl,Rattew2022-uh}.
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$, denoted by $|ψ̅_0⟩$, is supposed to be easily prepared on quantum circuits, and we assume that the eigenstate of $H̅_1$, which is continuously connected to $|ψ̅_0⟩$, corresponds to the target state $|ψ̅_1⟩$.
As the simplest case, we organize such a time-periodic Hamiltonian by
\begin{equation}
H(t) = \bar{H}_0 (1-\sin \omega t) + \bar{H}_1 \sin \omega t ,
\end{equation}
with satisfies $H(0)=H̅_0$ and $H(T/4)= H̅_1$.
Our algorithm for the adiabatic-like cases in Section \ref{Subsec:algorithm_adiabatic} enables us to efficiently execute the adiabatic state preparation.
Since the Fourier components of $H(t)$ is simply given by $H_0 = H̅_0$ and $H_±1= ±(H̅_0 - H̅_1)/2i$ with $m_=1$, the following oracles are necessary;
\begin{eqnarray}
O_0 &=& I \otimes \bar{O}_0, \quad \ket{G_0} = \ket{0} \ket{\bar{G}_0}, \\
O_{\pm 1} &=& \ket{0}\bra{0} \otimes ( \mp i \bar{O}_0 ) + \ket{1}\bra{1} \otimes (\pm i \bar{O}_1), \\
\ket{G_{\pm 1}} &=& \frac{\sqrt{\bar{\alpha}_0}\ket{0} \ket{\bar{G}_0}+ \sqrt{\bar{\alpha}_1}\ket{1} \ket{\bar{G}_1}}{\sqrt{\bar{\alpha}_0+\bar{\alpha}_1}}.
\end{eqnarray}
Here, $(O̅_0,|G̅_0⟩,α̅_0)$ and  $(O̅_1,|G̅_1⟩,α̅_1)$ respectively provides the block-encoding of the time-independent Hamiltonians $H̅_0$ and $H̅_1$ as Eq. (\ref{Eq2B:block_encode}).
The oracles $O_±1$ and $|G_±1⟩$ gives the block-encoding of $H_±1$ as $⟨G_±1 |O_±1|G_±1|=⟩ 2 H_±1/(α̅_0 + α̅_1)$.
According to the standard adiabatic theorem \cite{Albash2018RevModPhys}, the duration required for approximating the target state $|ψ̅_1⟩$ by the adiabatically-evolved state $U(T/4) |ψ̅_0⟩$ with accuracy $1-ε$ roughly amounts to $T ≳H̅_0 - H̅_1 / εΔ^2$, where $Δ$ denotes the minimal gap upon the instantaneous eigenstate of $H(t)$ continuously connecting $|ψ̅_0⟩$ and $|ψ̅_1⟩$.
Therefore, in the adiabatic state preparation based on our algorithm, the query complexity counted by the oracles of the two static Hamiltonians $H̅_0$ and $H̅_1$ amounts to
\begin{equation}
\frac{\bar{\alpha}_0 + \bar{\alpha}_1}{\varepsilon\Delta^2} \norm{\bar{H}_0 - \bar{H}_1} + o(\log (1/\varepsilon)).
\end{equation}
Although the second term is buried by the first one polynomially increasing in $1/ε$, the cost has better scaling in $1/ε$ compared to the cases where similar schedules for $H(t)$ are tackled 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 more sophisticated scheduling, such as local adiabatic interpolation \cite{Roland2002adiabatic}, boundary cancellation \cite{Lidar2009-bv,Wiebe2012-xx} and quasi-adiabatic processes \cite{Wan2020-at}.
In particular, the last one has achieved the query 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.

% Section: Discussion and Conclusion
\section{Discussion and Conclusion}\label{Sec:Discussion}

We conclude our paper with summarizing the results.
In this paper, we focus on time-dependent systems with time-periodicity, and organize an efficient implementation of their time-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 truncated Floquet-Hilbert space extract the time-evolved state $|ψ(t)⟩$ 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 best algorithm for time-independent systems \cite{Low2019-qubitization}, 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 than the time-evolution operator $U(t)= Texp(-i ∫_0^t H(t^') t^
In our algorithm, we exploit the qubitization technique to implement the exponential function, $e^-i H_^l_t$, in the effective Hamiltonian $H_^l_$.
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 \ref{Sec:truncated_Floquet_Hilbert} is valid also when a time-periodic Hamiltonian $H(t)$ is non-hermitian.
This suggests that our algorithm may be available also for solving time-periodic linear differential equations \cite{Fang2022-wq}, exemplified by dissipative quantum may-body systems \cite{Breuer2002-xj,Rivas2012-jb}.
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.

Another significant direction is to seek for useful tasks that can be efficiently tackled with time-periodic Hamiltonians.
As we have shown throughout the paper, time-periodic Hamiltonians can be simulated more efficiently generic time-dependent Hamiltonians with computational resources close to those for time-independent Hamiltonians.
This means that 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.

\section*{Acknowledgment}

We thank T. Kuwahara for fruitful discussion, especially on the Lieb-Robinson bound presented in Section \ref{Subsec:Lieb_Robinson_Floquet_Hilbert}.
This work is supported by MEXT Quantum Leap Flagship Program (MEXTQLEAP) Grant No. JPMXS0118067394, JPMXS0120319794, and JST COI-NEXT program Grant No. JPMJPF2014.

\bibliography{bibliography.bib}

% Appendix
\appendix
\begin{center}
\bf{\large Appendix}
\end{center}

\section{Short explanation of Eq. (\ref{Eq2A:psi_t_infinite})}
Our algorithm is mainly based on the fact that the dynamics in the Floquet-Hilbert space,
\begin{equation}
\ket{\psi^\infty(t)} = \lim_{l_\imax \to \infty} \sum_{l \in D^{l_\imax}} e^{-il\omega t}\bra{l} e^{-i \ms{H}_\eff^{l_\imax} t} \ket{0}\ket{\psi(0)},
\end{equation}
exactly corresponds to the solution $|ψ(t)⟩$ \cite{Levante1995-tx}.
Here, we roughly explain this relation for this paper to be self-contained.
We differentiate the above formula in $t$;
\begin{eqnarray}
&& i \dv{t} \ket{\psi^\infty(t)} \nonumber \\
&=& \sum_{l \in \bbZ} e^{-il\omega t} \bra{l} (\ms{H}_\eff+l\omega) e^{-i \ms{H}_\eff t} \ket{0}\ket{\psi(0)} \nonumber \\
&=& \sum_{l \in \bbZ} \sum_{m \in \bbZ} e^{-il\omega t} H_{-m} \bra{l+m} e^{-i \ms{H}_\eff t} \ket{0}\ket{\psi(0)} \nonumber \\
&=& \left( \sum_{m \in \bbZ} H_{-m} e^{im\omega t} \right)
\left( \sum_{l \in \bbZ}  e^{-il\omega t} \bra{l} e^{-i \ms{H}_\eff t} \ket{0}\ket{\psi(0)}\right) \nonumber \\
&=& H(t) \ket{\psi^\infty(t)}.
\end{eqnarray}
Noting that $|ψ^∞(0)⟩=|ψ(0)⟩$, the uniqueness of the solution of the Schr\"{o}dinger equation suggests the relation $|ψ^∞(t)⟩=|ψ(t)⟩$.

% Section
\section{Proof of the theorems for formulation}\label{Asec:Proof_theorems}

% Subsection
\subsection{Amplitude amplification by symmetry}\label{Asubsec2:amplification_symmetry}
In Section \ref{Subsec:amplification_symmetry}, we discuss the amplification of the time-evolved state $|ψ(t)⟩$ with exploiting the translation symmetry of generic time-periodic Hamiltonians.
With the usage of the protocol,
\begin{equation}
\ms{U}_\mr{amp1}^{l_\imax} (t) = (\ms{U}_\mr{ini}^{4l_\imax})^\dagger e^{-i t \sum_l l \omega \ket{l}\bra{l}} e^{-i \ms{H}_\eff^{4l_\imax} t} \ms{U}_\mr{ini}^{l_\imax},
\end{equation}
we see that it generates the time-evolved state $|ψ(t)⟩$ with $1$ amplitude as
\begin{equation}
\bra{0} \ms{U}_\mr{amp1}^{l_\imax} (t) \ket{0} \ket{\psi(0)} \simeq \frac{1}{2} (\ket{\psi(t)}+\order{\varepsilon}).
\end{equation}
We rigorously prove this statement, which is summarized as Theorem \ref{Thm5:amplification_symmetry}.

In order to show Theorem \ref{Thm5:amplification_symmetry}, we begin with discussing the approximate translation symmetry.
To be precise, we evaluate how much error appears in the approximation,
\begin{equation}
\bra{l} e^{-i\ms{H}_\eff^{4l_\imax}t} \ket{l^\prime} \simeq e^{il^\prime \omega t}\bra{l-l^\prime} e^{-i\ms{H}_\eff^{4l_\imax}t} \ket{0},
\end{equation}
which we refer to as Eq. (\ref{Eq5A:approximate_symmetry}) in the main text.
The exact upper bound on this error is given by the following theorem.

\begin{lemma}\label{LemmaA2:translation_symmetry}
\textbf{(Approximate translation symmetry)}

We choose the truncation order $l_\imax \in \Theta (\gamma t + \log (1/\varepsilon)/\log\log (1/\varepsilon))$ by Eq. (\ref{Eq4B:l_max_choice}).
For indices $l^\prime \in D^{l_\imax}$ and $l \in D^{4l_\imax}$, the inequality
\begin{eqnarray}
&& \norm{\braket{l |e^{-i \ms{H}_\eff^{4l_\imax} t}| l^\prime} - e^{i l^\prime \omega t} \braket{l \ominus l^\prime |e^{-i \ms{H}_\eff^{4l_\imax} t}| 0}} \nonumber \\
&& \qquad \leq 8 \frac{(\gamma t)^{n(l^\prime,l)}}{n(l^\prime,l)!} \label{EqA2:lemma_translation_symmetry}
\end{eqnarray}
is satisfied with
\begin{equation}\label{EqA2:Dyson_order_minimum}
n(l^\prime,l) = \left\lceil \frac{(8l_\imax - 2 m_\imax - |l| - |l^\prime|)}{m_\imax}\right\rceil.
\end{equation}
Here, $l \ominus l^\prime \in D^{4l\imax}$ denotes the difference $l-l^\prime$ defined modulo $8l_\imax$.
\end{lemma}

\textbf{Proof.---} We consider dynamics under a perturbed Hamiltonian $H_pert(t)$, defined by
\begin{eqnarray}
\ms{H}_\mr{pert}(t) &=& \ms{H}_\eff^{4l_\imax} + \ms{H}_b (t), \\
\ms{H}_b(t) &=& \sum_{(l,m) \in \partial F^{4l_\imax}} (\ket{l}\bra{l \oplus m} \otimes e^{i8l_\imax \omega t}H_{-m} + \text{h.c.}), \nonumber \\
\end{eqnarray}
with $∂F^4l_ = { (l,m)   |   4l_- m_+1 ≤l ≤4l_,   4l_- l  +1 ≤m ≤m_}$.
In the interaction picture based on $H_0 = ∑_l ∈D^4l_ |l⟩⟨0| ⊗(H_0 -lω)$, the perturbed Hamiltonian $H_pert,I(t)$ is exactly translation symmetric under arbitrary shift $|l⟩ →|l ⊕m⟩$ as
\begin{equation}\label{EqA3:H_eff_pert_int}
\ms{H}_{\mr{pert},I}(t) = \sum_{m \neq 0; |m| \leq m_\imax} \mr{Add}_m^{4l_\imax} \otimes e^{- im\omega t} H_m^I (t).
\end{equation}
In the original frame, this exact symmetry implies the satisfaction of
\begin{equation}
\braket{l|\ms{U}_\mr{pert}(t)|l^\prime} = e^{i l^\prime \omega t} \braket{l \ominus l^\prime |\ms{U}_\mr{pert}(t)|0},
\end{equation}
where $U_pert(t)$ indicates the time evolution operator under the perturbed Hamiltonian $H_pert(t)$.
Upon this relation, we evaluate the upper bound on the error by a triangle inequality,
\begin{eqnarray}
\text{[l.h.s. of Eq. (\ref{EqA2:lemma_translation_symmetry})]} &\leq& \norm{\braket{l |e^{-i \ms{H}_\eff^{4l_\imax} t} - \ms{U}_\mr{pert}(t)|l^\prime}} \nonumber \\
&\quad& + \norm{\braket{l \ominus l^\prime |e^{-i \ms{H}_\eff^{4l_\imax} t} - \ms{U}_\mr{pert}(t)|0}}. \nonumber \\
&& \label{EqA2:Lemma_symmetry_proof_1}
\end{eqnarray}

We begin with computing the bound of the first term.
Using the Dyson series expansion in the interaction picture, similarly to Eqs. (\ref{Eq4A:def_reference_hamiltonian})-(\ref{Eq4A:Dyson_transition_path}), it can be bounded by
\begin{eqnarray}
&& \norm{\braket{l |e^{-i \ms{H}_\eff^{4l_\imax} t} - \ms{U}_\mr{pert}(t)|l^\prime}}  \nonumber \\
&& \qquad \leq  \sum_{n=0}^\infty \int_0^t \dd t_n \hdots \int_0^{t_2} \dd t_1 \nonumber \\
&& \qquad \qquad \norm{\bra{l} \prod_{i=1}^n \ms{H}_{\mr{pert},I}(t_i) - \prod_{i=1}^n \ms{H}_I^{4l_\imax}(t_i)\ket{l^\prime}}. \nonumber \\
&& \label{EqA2:Lemma_symmetry_proof_2}
\end{eqnarray}
As discussed in the proof of Theorem \ref{Thm4:Bound_transition_rate}, each integrand is decomposed into contributions from the transition amplitudes via paths $|l^'⟩ →|l_1⟩ →→|l⟩$.
The difference between the Hamiltonians $H_pert(t)$ and $H_^4l_$ arises only when the path goes across $D^4l_ \D^4l_- m_$, which is the support of $H_b(t)$.
For low order terms with $n < n (l^', l)$, defined by Eq. (\ref{EqA2:Dyson_order_minimum}), such a nontrivial path is absent, and we obtain
\begin{eqnarray}
&& \norm{\braket{l |e^{-i \ms{H}_\eff^{4l_\imax} t} - \ms{U}_\mr{pert}(t)|l^\prime}}  \nonumber \\
&& \qquad \leq  \sum_{n=n(l^\prime,l)}^\infty \frac{t^n}{n!} \left\{ \left(\sup_t (\norm{\ms{H}_{\mr{pert},I}(t)}) \right)^n + \gamma^n \right\} \nonumber \\
&& \label{EqA2:Lemma_symmetry_proof_3}
\end{eqnarray}
Since the Hamiltonian $H_pert,I(t)$ given by Eq. (\ref{EqA3:H_eff_pert_int}) is a circulant matrix (i.e. invariant under the shift of $|l⟩$) \cite{Gray2006-da}, its operator norm can be bounded by
\begin{eqnarray}
\norm{\ms{H}_{\mr{pert},I}(t)} &\leq& \sup_{t^\prime} \left( \norm{\sum_{m \neq 0} e^{-i m \omega t} H_m^I(t) e^{im \omega t^\prime} }\right) \nonumber \\
&=& \sup_{t^\prime} \left( \norm{H(t-t^\prime) -H_0} \right) \nonumber \\
&=& \gamma.
\end{eqnarray}
As a result, Eq. (\ref{EqA2:Lemma_symmetry_proof_3}) is further bounded by
\begin{equation}
\norm{\braket{l |e^{-i \ms{H}_\eff^{4l_\imax} t} - \ms{U}_\mr{pert}(t)|l^\prime}} \leq 2 \sum_{n=n(l^\prime, l)}^\infty \frac{(\gamma t)^n}{n!}.
\end{equation}