subsecref subsecname = thmref thmname = theorem lemref lemname = lemma fig
refcmd=Figure LABEL:#1 tab refcmd=Table LABEL:#1 eq refcmd=(LABEL:#1) subsec
refcmd=Subsection LABEL:#1 sec refcmd=Section LABEL:#1 app refcmd=Appendix
LABEL:#1 part refcmd=Part LABEL:#1 secapp refcmd=Appendix LABEL:#1 subsecapp
refcmd=Appendix Subsection LABEL:#1

# Influence of Trotterization error on single-particle tunneling

Anton V. Khvalyuk<EMAIL_ADDRESS>LPMMC, University Grenoble-
Alpes, France  Kostyantyn Kechedzhi  Google Research, Mountain View, CA, USA
Vadim S. Smelyansky  Google Research, Mountain View, CA, USA  Lev. B Ioffe
Google Research, Mountain View, CA, USA

(August 28, 2024)

###### Abstract

Simulation of the single-particle tunneling problem by means of the Suzuki-
Trotter approximation (STA) is analyzed. The target system describes a
particle hopping across a chain of sites with position-dependent potential
profile. The latter is assumed to be smooth and posses several local minima
separated by a potential barrier, arranging a tunneling problem between the
localized states in different minima. The STA error is found to manifest
itself in three ways: _i)_ perturbative energy shifts, _ii)_ nonperturbartive
renormalization of the tunneling rates, and _iii)_ perturbative leakage of
total probability to other states. In the general case, the first type of
error is the most essential, as the emerging detuning of the tunneling
resonance has to be compared with exponentially small tunneling rates. In
absence of detuning (e.g. if the resonance is protected by symmetry), STA is
found to cause exponential enhancement of the tunneling rates. The last type
of error classifies the overall defect in the wave function and thus
delineates the region of sufficiently weak distortion of the wave function due
to STA. The conducted analysis confirms the naive criteria of applicability
$\max\\{T,P\\}\ll\delta t^{-1}$ (with $T,P$ being the typical scales of
kinetic and potential terms, respectively), while also revealing the structure
of error and its actual behavior with system parameters. Analysis of the case
of large Trotter step is also performed, with the main result being the
reconstruction of the low-energy spectrum due to coupling between states with
energy difference close to $2\pi/\delta t$. The connection of the obtained
results with the rigorous upper error bounds on the STA error is discussed,
with particular emphasis on the reasons for the fact that these rigorous
bounds are not always saturated. We also point out that the proposed problem
can be directly implemented on existing quantum devices [1]. In particular, we
give a detailed description of an experimental design that demonstrates the
described physics.

##  I Introduction

Simulating many-body quantum systems is a notoriously difficult problem that
is often beyond the capabilities of even the most powerful classical
supercomputers. In addressing this challenge, developing well-controlled
quantum system that can further be used to simulate a wide range of other
quantum systems has been a promising approach [2]. One example of such devices
is given by superconducting systems [3, 4] that realize a 2D grid of two-level
systems (qubits) with the possibility to apply any unitary operation involving
one or two spatially adjacent qubits. These unitary operations (so-called
gates) can be applied simultaneously for nonintersecting sets of qubits, which
in the main idea behind all applications. In particular, such devices are a
natural platform for simulating 1D and 2D spin-1/2 systems, as several
prominent examples demonstrate [1, 5], but other computationally expensive
tasks can also executed [3]. However, several main obstacles stand in the way
of practical applications of these devices [1]: _i)_ dephasing and decoherence
of many-body state caused by noise from coupling to external devices, _ii)_
errors in unitary gates and in readout of the final state, and _iii)_
representation of Hamiltonian evolution in terms of local unitary gates, with
the latter being the topic of the present paper.

The primary challenge in simulating the evolution of a system with Hamiltonian
$H$ is that the evolution operator $\hat{U}=e^{-itH}$ is generally truly
nonlocal unitary operator, whereas any unitary gate available on existing
platforms is local, so the exact evolution can only be restored with infinite
number of gates. However, numerous algorithms [6] have been proposed to tackle
this problem. All such algorithms essentially represent a cleverly designed
sequence of local one- and two-qubit gates whose combined action
$U_{\text{appr}}$ approaches that of the target evolution operator $\hat{U}$
as the number of elementary steps in the sequence $M$ is increased. The
convergence towards $\hat{U}$ with the increase $M$ is established by deriving
a rigorous upper bound on the error measured, e.g., as the spectral norm
distance $\left\|U-U_{\text{appr}}\right\|$. These schemes can also be viewed
as quantum-mechanical analogs of Hamiltonian-blind symplectic integrators of
Hamiltonian equations of motion that guarantees an upper bound on the error of
the result. Ref [6] provides a concise review of existing approaches and the
associated expressions for the upper bounds. Although it can be proved that
any single algorithm cannot perform accurate simulation of _all_ Hamiltonians
with number of gates that growth slower than linearly with the required
simulation time and\or system size [7, 8], the available algorithms are
capable of coming arbitrarily close to the optimal linear scaling [6].

It is often the case, however, that the analytical upper bounds are not
saturated by the actual evolution [6, 9, 10], so the natural question of the
true value of the error emerges. Moreover, the total error in the wave
function (measured e.g. as the square norm distance) does not necessarily
translate to a comparable error in physical observables, as those might be
largely insensitive to the sectors of the Hilbert space that host the largest
part of the error. On the other hand, the idea of quantum simulations is
beneficial for analysis of behavior of the systems that are not amenable for
classical simulation, in which case it is also crucial to gain qualitative
understanding of the errors brought in by the simulation routine. Detailed
error analysis thus also serves the purpose of reducing of the required number
of operations by accurately reproducing only the evolution of the part of the
Hilbert space that plays the most important role in the target physical
observables.

The present paper is focused on the detailed analysis of the error induced by
the Suzuki-Trotter algorithm [11], also known as the product formula approach
[6] or Trotterization. This algorithm has two practical advantages: _i)_ it
requires no qubits in excess of those directly storing the configuration of
the target spin system, and _ii)_ it approximately conserves the locality of
the underlying Hamiltonian, which is also to the benefit of its theoretical
analysis. The particular emphasis is made on the distortion of the tunneling
processes, as those appear to play a role in quantum algorithms for
combinatorial problems [12, 13, 14], although here we only focus on a single-
particle model with sufficiently simple dynamics.

One particularly important question about the Suzuki-Trotter algorithm is the
interpretation of the distortion brought in by this algorithm. The latter
reproduces the evolution by performing a sequence of discrete steps, the thus
introducing a time scale of the order of the Trotter step duration $\delta t$.
As a result, it is tempting to regard the distortion of the evolution as
perturbation at high frequencies ~$2\pi n/\delta t,\,\,\,n\in\mathbb{N}$, so
one expects heating, assisted tunneling, and other irradiation-induced effects
to occur. On the other hand, for a time-independent target Hamiltonian, the
Suzuki-Trotter evolution is exactly periodic with period $\delta t$, so the
Floquet theorem actually guarantees the existence the effective time-
independent Hamiltonian $H_{\text{eff}}$ that exactly corresponds to evolution
during a single Trotter step. $H_{\text{eff}}$ is thus conserved during the
Suzuki-Trotter evolution, which at first glance is in contradiction with
heating. For a general setting, this issue has already been addressed in the
literature [15], but in this work we intentionally avoid the problem
altogether by restricting ourselves to a simple model that does not feature
many-body physics.

The paper is organized as follows: the following II describes the model
Hamiltonian and the associated Suzuki-Trotter algorithm. It then formulates
more precisely two flavors of the tunneling problem to be discussed and
provides a concise review in the semiclassical picture of the system’s
dynamics. III then exposes the detailed description of the error introduced by
the simulation algorithm. It describes the three most important effects:
resonance detuning due to the energy shift (III.4), enhancement of the
tunneling rate (Subsections III.5 and III.6), and probability defect due to
leak from the target initial state (III.7). It also contains description of
certain effects of sufficiently large Trotter steps (III.8). Finally, IV
summarizes the key qualitative findings of the work (IV.1), conducts an
instructive comparison with the rigorous upper bounds known from the
literature (IV.2), provides a conceptual description of the experimental setup
that demonstrates the aforementioned effects (IV.3), and outlines certain
directions of further research (IV.4). The paper is supplemented by three
Appendices describing the technical details of the calculations and the design
of the experiment.

##  II Solvable model of Trotterized evolution

###  II.1 Hamiltonian

We consider the Hamiltonian of inhomogeneous chain of spins $1/2$ with $XY$
nearest-neighbor interaction:

$H=\sum_{i=1}^{L}h_{i}S_{i}^{z}-\sum_{i=1}^{L-1}J\left(S_{i}^{x}S_{i+1}^{x}+S_{i}^{y}S_{i+1}^{y}\right),$
(1)

where $h_{i}$ is position-dependent, $J>0$ is a constant. The spatial profile
of $h_{i}$ is assumed to be smooth, as explained below in more detail.

The Hamiltonian (1) conserves the $z$-projection of the total spin
$S_{\text{total}}^{z}=\sum_{i}\sigma_{i}^{z}/2$, which allows one to map the
problem to that of non-interacting fermions by means of the Jordan-Wigner
transformation [16]:

$S_{i}^{z}=a_{i}^{\dagger}a_{i}-\frac{1}{2},\,\,\,\,S_{i}^{x}S_{i+1}^{x}+S_{i}^{y}S_{i+1}^{y}=\frac{a_{i}^{\dagger}a_{i+1}+a_{i+1}^{\dagger}a_{i}}{2},$
(2)

where $a_{i}^{\dagger},a_{i}$ are fermionic creation and annihilation
operators for site $i$. The vacuum state with no fermions corresponds to the
state with all spins down
$\left|\Omega\right\rangle=\otimes_{i}\left|\downarrow_{i}\right\rangle$, and,
more generally, the subspace with $N$ fermions reproduces the subspace with
total $z$-projection $S_{\text{total}}^{z}=\sum_{i}S_{i}^{z}=N-L/2$.

###  II.2 Approximate evolution by Suzuki-Trotter algorithm

Figure 1: Circuits of the three possible orderings of non-commuting operators
in Suzuki-Trotter formula. Horizontal black lines correspond to qubits, and
rectangular blocks represent gates executed on the corresponding qubits at a
given moment of time (horizontal axis). Each circuit represents one step of
Suzuki-Trotter algorithm according to Eq. (7): _left_ :
$\left\\{P,K_{\text{odd}},K_{\text{even}}\right\\}$, _center_ :
$\left\\{K_{\text{odd}},P,K_{\text{even}}\right\\}$, _right_ :
$\left\\{K_{\text{odd}}+\left(1-\alpha\right)P,K_{\text{even}}+\alpha
P\right\\}$ for $\alpha=1/2$. Written inside or near a blocks are the
operators whose exponent is being applied within this block. Below each
vertical cross-section of the circuits written are the parts of the
Hamiltonian applied by the corresponding layer of the circuit. For all three
circuits, the central parts of the Suzuki-Trotter formula (7) have been
coalesced as $e^{-i\delta tA_{n}/2}e^{-i\delta tA_{n}/2}=e^{-i\delta tA_{n}}$.
In all three circuits, obvious optimizations are possible for multiple Trotter
steps (such as, e.g. merging the two adjacent single-qubits gates into a
single one, and so on). For the right circuit, the scheme inside the block
indicates that to reproduce the target two-qubit operator a sequence of
several one- and two-qubit gates is applied on the hardware. Note also the
$\delta tP/2$ gates (green) applied at both end qubits on the left circuit – a
consequence of open boundary conditions.

The problem of simulating the evolution under the Hamiltonian $H$ during time
$t$ amounts to approximating the unitary $U=\exp\left\\{-itH\right\\}$ with
sufficient accuracy by a product of local one- and two-qubit gates available
on the quantum hardware. This a nontrivial task due to the noncommutativity of
various local terms in the Hamiltonian (1). Indeed, Eq. (1) can be rewritten
as

$H=K_{\text{even}}+K_{\text{odd}}+P,$ (3)
$K_{\text{even}(\text{odd})}=\sum_{k:\text{even}(\text{odd})}H_{k,k+1},\,\,\,P=\sum_{k=1}^{L}H_{k},$
(4)
$H_{k,k+1}=-J\left(S_{k}^{x}S_{k+1}^{x}+S_{k}^{y}S_{k+1}^{y}\right),\,\,\,H_{k}=h_{k}S_{k}^{z}.$
(5)

The operators within each term in Eq. (4) commute with each other, but the
three terms $K_{\text{even}},K_{\text{odd}},P$ are mutually non-commuting. As
a result, on quantum hardware one can apply only the evolution operators
corresponding to $aK_{\text{even}}$, $bK_{\text{odd}}$ and $cP$ with some real
coefficients $a,b,c$. The operators of the form $aK_{\text{even}}+bP$ and
$cK_{\text{odd}}+dP$ with real $a,b,c,d$ can be further represented as a
collection of unitaries acting on nonintersecting pairs of qubits. Although at
the hardware level each such unitary is typically implemented as a combination
of one- and two-qubit gates, such grouping allows one to remove part of the
error from noncommutativity of $K_{\text{even}(\text{odd})}$ and $P$. On the
other hand, such grouping might potentially increase the share number of gates
depending on the compilation algorithm, so we will consider both cases: with
and without such grouping.

The second order Suzuki-Trotter formula [11, 7] provides the approximation for
the evolution operator $U=e^{-itH}$ for $H=A_{1}+..+A_{n}$, where each of
$A_{k}$ can be exponentiated exactly. The whole evolution is partitioned into
identical steps of length $\delta t$ (assuming that $t/\delta t$ is an
integer):

$U_{\text{appr}}\left(t\right)=\left\\{U_{\text{appr}}\left(\delta
t\right)\right\\}^{t/\delta t},$ (6)

with each step performed according to the approximate evolution operator

$U_{\text{appr}}\left(\delta t\right)=e^{-iA_{n}\delta t/2}...e^{-iA_{1}\delta
t/2}\,\times\,e^{-iA_{1}\delta t/2}...e^{-iA_{n}\delta t/2},$ (7)

where $A_{k},k=1,..,n$ correspond to the set of applicable operators in any
fixed order, and the second term in Eq. (7) features product of operators in
reverse direction relative to that of the first term. In the limit $\delta
t\rightarrow 0$, the approximate evolution operator
$U_{\text{appr}}\left(t\right)$ approaches the target one $e^{-itH}$, while
finite $\delta t$ causes a distortion of the evolution, which can be
characterized by the effective Hamiltonian:

$U_{\text{appr}}\left(t\right)=e^{-itH_{\text{eff}}},\,\,\,H_{\text{eff}}=-\frac{1}{i\delta
t}\text{Ln}\,U_{\text{appr}}\left(\delta t\right),$ (8)

so the actual evolution after $t/\delta t$ steps is described by
$U_{\text{appr}}\left(t\right)=\exp\left\\{-iH_{\text{eff}}t\right\\}$. The
result of Ref. [11] is that

$H_{\text{eff}}-H=O\left(\delta t^{2}\right).$ (9)

The Suzuki-Trotter algorithm can be viewed as the quantum-mechanical
counterpart of the symplectic Runge Kutta method [17] for classical
Hamiltonian equations of motion, thus also sharing certain features, such as
naturally occurring energy conservation (i.e., symplecticity), locality and
existing of the effective semiclassical Hamiltonian.

The appropriate values of $\delta t$ are determined from the naive comparison
of the characteristic energy scales:

$\left|J\delta t\right|,\left|P\delta t\right|\ll 1,$ (10)

where $P=\max_{i}h_{i}-\min_{i}h_{i}$ is the magnitude of the potential term,
and $J$ is the order of magnitude of the kinetic term. However, if the problem
contains other relevant energy scales, one should also require these scales to
be small compared to $\delta t^{-1}$ , which might lead to a much stronger
requirement on the step size $\delta t$, as it is shown below.

Eq. (7) admits different splittings of the total Hamiltonian (3) into allowed
operators, with three conceptually distinct options:
$\left\\{K_{\text{even}},K_{\text{odd}},P\right\\}$,
$\left\\{K_{\text{even}},P,K_{\text{odd}}\right\\}$, and
$\left\\{K_{\text{even}}+\alpha
P,K_{\text{odd}}+\left(1-\alpha\right)P\right\\}$ for some $\alpha\in[0,1]$.
All other groupings differ from these three either by redefining the
“beginning” of each Trotter step to be at a different gate, or by
interchanging $K_{\text{even}}$ and $K_{\text{odd}}$. Moreover, the structure
of and the effect from the resulting error term is similar for both cases, as
it will be shown below, so the choice of the particular design boils down to
hardware considerations. The visualization of the algorithm is shown on 1.

Figure 2: Sketch of the potential profile (filled blue curve), the resulting
eigenlevels (horizontal dashed lines) and wave functions (filled curves) for
the two tunneling problems in question. For a given energy, the potential
defines a set of classically allowed and forbidden regions. If the allowed
region is bounded, it hosts a series of discrete levels with energy-dependent
level spacing $\Delta$ that is determined by the particular geometry of the
allowed region, typically being inversely proportional to the size of the
allowed region at a given energy. _Left_ : double-well potential with two sets
of discrete energy levels. Due to finite height of the potential barrier, the
resonant states in the two wells are hybridized with small energy splitting
$\eta$. The resonance between left and right metastable state can be either
accidental (and thus characterized by finite detuning $\varepsilon\sim\eta$)
or induced by potential symmetry. _Right_ : potential with a single metastable
minimum separated from large region with dense quasi-continuous spectrum. Each
bound state inside the local minimum is coupled to the quasi-continuous
spectrum of states, inducing a small level width $\Gamma$, which is assumed to
be much larger than the level spacing $\delta$ of the continuum.

Because the original problem is non-interacting, the problem of calculating
$H_{\text{eff}}$ can be reduced to a similar problem for just a single
fermion. Indeed, at each phase of an individual Trotter step (7), the
fermionic system evolves according to a Hamiltonian that is quadratic in
fermionic operators:

$A_{k}=\sum_{ij}a_{i}^{\dagger}A_{k}^{ij}a_{j}$ (11)

where the sum goes over all sites $i,j$, and the matrix elements $A_{k}^{ij}$
are obtained from Eqs. (4)-5 by substitution (2). Due to the Wick’s theorem,
the composition (7) is still generated by a quadratic Hamiltonian:

$U_{\text{appr}}\left(\delta t\right)=\exp\left\\{-i\delta
t\sum_{ij}a_{i}^{\dagger}\mathcal{H}_{\text{eff}}^{ij}a_{j}\right\\}$ (12)

where the effective single-particle Hamiltonian $\mathcal{H}_{\text{eff}}$ is
defined via its evolution operator $\mathcal{U}=e^{-i\delta
t\mathcal{H}_{\text{eff}}}$ by the same Suzuki-Trotter formula (7), but with
each $A_{k}$ replaced by a single-particle matrix according to the
correspondence

$H_{k,k+1}\mapsto\mathcal{H}_{k,k+1}^{ij}=-\frac{J}{2}\left(\delta_{ik}\delta_{j,k+1}+\delta_{i+1,k}\delta_{jk}\right),$
(13)
$H_{k}\mapsto\mathcal{H}_{k}^{ij}=\left(h_{k}-\frac{1}{2}\right)\delta_{ik}\delta_{jk}.$
(14)

The behavior of a given spin state under the evolution with the approximate
unitary is thus entirely described by the single-particles effects and the
Fermi statistics of the Jordan-Wigner fermions.

###  II.3 Single-particle tunneling problem

In what follows, we are interested in the single spin flip subspace
corresponding to just a single fermion, which means that our problem is
described by a single-particle Hamiltonian of fermions:

$\mathcal{H}=\sum_{k=1}^{L}h_{k}\,\left|k\left\rangle\right\langle
k\right|-\sum_{k=1}^{L-1}\frac{J}{2}\left(\left|k\left\rangle\right\langle
k+1\right|+\left|k+1\left\rangle\right\langle k\right|\right)$ (15)

where we have dropped the additional $1/2$ after $h_{k}$ as it produces a
constant energy shift. Unless stated otherwise, we consider smooth potential
profiles of the form $h_{n}=h\left(an\right)$, where $h\left(x\right)$ is a
smooth function, and $a\ll 1$ plays the role of small lattice constant. We are
then interested in describing the problem of tunneling outside the local
minimum of the potential. In terms of the original spins, this corresponds to
creating a single-spin excitation with spin density localized in a local
minimum of $h\left(x\right)$ and observing how the spin density tunnels
through the potential barrier. The two potential profiles of particular
interest are sketched on 2 and discussed below. In both cases, we assume that
the initial state $\left|\psi_{\text{init}}\right\rangle$ of the excitation is
close to a certain localized metastable state $\left|\psi_{N}\right\rangle$ of
a local minimum of the potential.

##### 1\. Rabi oscillations.

Consider the double-well potential on 2, left. If the barrier is infinite,
each well hosts a series of discrete energy levels characterized by typical
level spacing $\Delta$ (energy-dependent for an anharmonic well). Suppose that
a pair of those levels from different wells has a small energy difference
$\varepsilon=E_{\text{right}}-E_{\text{left}}\ll\Delta$, either due to
reflection symmetry of the potential (in which case $\varepsilon=0$) or by
accident. Finite transparency of the potential barrier then leads to
hydridization of those levels according to the following Hamiltonian of the
resonant subspace:

$H_{\text{res.}}=\begin{pmatrix}-\varepsilon/2&\eta\\\
\eta&\varepsilon/2\end{pmatrix}$ (16)

where $\eta\ll\Delta$ is the tunneling amplitude between the two resonant
states, exponentially small for a large barrier.

This problem is readily implementable on existing quantum devices [1].
However, the approximate evolution operator (7) distorts the evolution in two
ways: _i)_ the energies of the discrete states in each of the wells and the
tunneling amplitude are distorted, leading to a change in the resonant
subspace Hamiltonian (16), and _ii)_ the wave functions of the discrete states
are distorted as well, so an arbitrary initial state $\left|\psi\right\rangle$
localized in either of the wells now has a different overlap with the target
resonance subspace. Our task is to to describe these effects.

##### 2\. Decay into continuous spectrum.

Consider now a local minimum separated from a large region of small potential
by a potential barrier, 2, right. The well still hosts a set of discrete
metastable states with typical level spacing $\Delta$, and the region beyond
the barrier is characterized by a nearly continuous spectrum with level
spacing $\delta\ll\Delta$. Finite width of the barrier arranges a coupling of
the two regions by a finite transmission amplitude $\eta$, which we assume to
be much smaller than the discrete level spacing $\Delta$. While
$\delta\gg\eta$ corresponds to previously considered case of discrete
resonance, the opposite situation $\delta\ll\eta$ causes the decay of the
metastable state due to coupling with the continuous spectrum. As a result,
each discrete metastable state in the well is now characterized by an
exponential decay of probability
$\left|\psi_{N}\left(t\right)\right|^{2}\propto\exp\left\\{-\Gamma_{N}t\right\\}$
with small characteristic decay rate $\Gamma_{N}\sim\eta^{2}/\Delta\ll\Delta$,
and after $t\sim\Gamma_{N}^{-1}$ such state decays into states of continuous
spectrum with energies close to $E_{N}$ according to a Lorentzian profile of
width $\sim\Gamma_{N}$.

Attempting to simulate this process using the Suzuki-Trotter approach distorts
this evolution by _i)_ introducing corrections to the values of
$E_{N},\Gamma_{N}$ and _ii)_ by altering the wave functions inside the well,
which changes the distribution of energies in the initial state. We are
interested in describing these changes. Although this problem is not practical
for existing hardware, as one needs exponentially large system sizes to render
the continuous part of the spectrum sufficiently dense to satisfy
$\delta\ll\Gamma$, it still presents interest from the physics point of view.

###  II.4 Semiclassical approximation for single-particle tunneling

In this section, we review some of the well-known results of the semiclassical
approximation in application to the problems of II.3, with the intention to
refer to these results in what follows. The semiclassical approach is a
conventional tool for describing the events of tunneling through a large
barrier. It consists of two parts [18]: _i)_ constructing semiclassical wave
functions in the region where the semiclassical approximation is applicable
and _ii_) asymptotic matching of two regions of applicability with a transfer
matrix obtained from exact solution of the neighborhood of the problematic
region. In our problem, there is a complication in the form of discretized
spatial dimension, but this does not reduce the power of the semiclassical
approach [19]. A detailed overview of the procedure for the particular system
in question is given in Discrete-Semiclassical-Approximaion.

The WKB approximation amounts to the following asymptotic series w.r.t the
small lattice constant $a\ll 1$ (being a closed analog of the Planck’s
constant) for the stationary wave function:

$\psi_{n}=\frac{e^{iS\left(an\right)/a}}{\sqrt{v\left(an\right)}}\left[1+O\left(a\right)\right],$
(17)

where $v,S$ are assumed to be smooth functions of $x$. The classical action
$S\left(x\right)$ obeys the classical Hamilton-Jacobi equation:

$\mathcal{H}\left(x,p\right)=E,\,\,\,p=S^{\prime}\left(x\right),$ (18)

which also defines the classical momentum $p\left(x\right)$. Here, the
classical counterpart of the Hamiltonian is defined through the Wigner
transform:

$\mathcal{H}\left(x=an,p\right)=\sum_{m}\mathcal{H}^{nm}e^{ip\left(m-n\right)}=-J\cos
p+h\left(x\right)$ (19)

where the sum in the r.h.s of the middle expression goes over all sites, and
the last expression corresponds to the particular case of Hamiltonian (15).
Finally, $v\left(x\right)$ in Eq. (17) is the classical group velocity:

$v\left(x\right)=\left.\frac{\partial\mathcal{H}}{\partial
p}\right|_{p=p\left(x\right)}=J\sin p\left(x\right)$ (20)

The approximation (17) is only applicable in the region where the relative
change of the de Broglie wave length at the spatial scale of the wave length
itself is small, i.e., the wave length is well defined locally:

$\frac{a\,h^{\prime}\left(x\right)}{p^{2}\left(x\right)\,v\left(x\right)}\ll
1.$ (21)

This condition is violated if the potential is not smooth (i.e.,
$h^{\prime}\left(x\right)\sim 1/a\gg 1$) or when either the classical momentum
$p$ or the group velocity $v$ vanishes. Such points are called the classical
turning points. In addition to the standard turning points $p,v=0$, there are
also anomalous turning points corresponding to $v=0,p=\pm\pi$, which is the
consequence of the discreteness of the problem. The existence of this points
is the semiclassical explanation of the Bloch’s oscillations [20].

In the regions where the semiclassical approximation breaks down, the
corresponding stationary Schrodinger equation has to be solved exactly. The
asymptotic behavior of the result is then matched to the semiclassical
expressions (17), and the latter are required to be applicable in a
sufficiently large region, i.e. $\Delta x\gg\lambda\left(x\right)$, which also
implies $S/a\sim\Delta x/\lambda\gg 1$. The procedure yields a connection
between the coefficients of two linearly independent semiclassical solutions
on either side of the non-classical region. The example of this procedure for
a simple turning point is presented in Matching-and-quantization-conditions.

The semicalssical procedure allows one to compute the answers to the two
problems described in II.3, as explained in Matching-and-quantization-
conditions or otherwise known from the literature [19, 18, 21]. The discrete
energy levels in a potential well are described by the celebrated Bohr
quantization rule:

$\frac{S_{12}}{a}=\frac{1}{a}\intop_{x_{1}}^{x_{2}}p\left(E_{N};x\right)dx=\pi\left(N+\frac{1}{2}\right)+\theta_{2}-\theta_{1},$
(22)

where $N\in\mathbb{N}\cup\left\\{0\right\\}$ is the level quantum number,
$p\left(E,x\right)$ is the classical momentum at energy $E$ found from Eq.
(18), $x_{1,2}$ are the positions of the classical turning points defined by
vanishing group velocity $v\left(x_{i}\right)=0$, and $\theta_{1,2}$ are the
additional phases introduced by the anomalous turning points [19]:

$\theta_{i}=\begin{cases}0,&p\left(x_{i}\right)=0,\\\ -\frac{\pi}{2}+\pi
x_{i}/a,&p\left(x_{i}\right)=\pm\pi.\end{cases}$ (23)

The phases $\theta_{i}$ in Eqs. (22) and (23) actually depend on the exact
position $n=x/a\in\mathbb{R}$ of the anomalous turning point _between_ the two
sites, which is essential for correctly reproducing the quantization
conditions in problems with discrete coordinate, as pointed out in Ref [19]
and explained in Matching-and-quantization-conditions.

Eq. (22) also produces the well-known expression for the single-particle level
spacing $\Delta=E_{N+1}-E_{N}$:

$T_{12}\left(E\right)\Delta\left(E\right)=2\pi a$ (24)

where $T_{12}=2\,\partial S_{12}/\partial E$ is the period of classical motion
at energy $E$.

Eqs. (22-24) are applicable if the region of applicability of the
semiclassical approximation inside the well is sufficiently large, i.e.
$S_{12}/a\gg 1$. This condition is typically violated for the lowest energy
levels of the well, so one is required to solve the eigenproblem exactly in
this case. However, even for such cases the parametric dependence of the
semiclassical result is correct up to a numeric coefficient of the order of
unity, and it turns out to be exact if the potential minimum is described by a
quadratic approximation.

The tunnel splitting of two semiclassical levels with close energies is
described by

$\eta=\frac{\sqrt{\Delta_{L}\Delta_{R}}}{2\pi}\exp\left\\{-\frac{1}{a}S_{B}\left(\frac{E_{L}+E_{R}}{2}\right)\right\\}$
(25)

where $\Delta_{L,R}$ are the local level spacings at energies $E_{L,R}$ in the
corresponding wells, and

$S_{B}\left(E\right)=\intop_{x_{L}}^{x_{R}}\left|p\left(E;x\right)\right|dx$
(26)

is the under-barrier action calculated from the imaginary momentum of the
forbidden region between the classical turning points $x_{L},x_{R}$ such that
$v\left(x_{L,R}\right)=0$. Finally, the same method gives the celebrated
Gamow’s formula [21, par. 50] for the decay rate:

$\Gamma_{N}=\frac{\Delta_{N}}{2\pi}\exp\left\\{-\frac{2}{a}S_{B}\left(E_{N}\right)\right\\}$
(27)

with $\Delta_{N}$ being the level spacing of the metastable states inside the
well, and $S_{B}$ given by the same Eq. (26) at the energy $E_{N}$ of the
corresponding metastable state. If the states inside the well are not
described by the semiclassical approximation (e.g., several first low-energy
states of the potential minimum), both expressions (25) and (27) typically
remain valid up to a numerical factor of the order unity, with the case of a
quadratic minimum of the potential described exactly.

##  III Influence of trotterization on single-particle tunneling

###  III.1 Matrix elements of the effective Hamiltonian

By construction, each operator applied within one Suzuki-Trotter step (7) has
nonzero matrix elements only at distances at most two, i.e., $\left\langle
i\left|\exp\left\\{-iA_{k}\delta t/2\right\\}\right|j\right\rangle=0$ for
$\left|i-j\right|>2$, as it is composed of commuting two-qubit gates.
Moreover, due to the smallness of the Trotter step $\delta t$, nondiagonal
matrix elements ought to be small as $O\left(J\delta t\right)$, with $J$ being
the characteristic scale of the kinetic energy. For our case, the single-
particle evolution operator $\mathcal{U}=e^{-i\delta
t\mathcal{H}_{\text{eff}}}$ is composed of either 5 such operator, two of
which are diagonal (as they represent the potential term), or just 3
operators. As a result, it can only have nonzero matrix elements at distances
not exceeding 3, and these matrix elements bear the corresponding order of
smallness in powers of $J\delta t$:

$\mathcal{U}^{ij}\sim\begin{cases}\left(J\delta
t/4\right)^{\left|i-j\right|},&\left|i-j\right|\leq 3,\\\
0,&\left|i-j\right|>3.\end{cases}$ (28)

It then follows from direct expansion of Eq. (8) in powers of $\delta t$ that
large-distance matrix elements of the effective Hamiltonian $H_{\text{eff}}$
decay exponentially with distance:

$\mathcal{H}_{\text{eff}}^{ij}\sim J\cdot\left(J\delta
t/4\right)^{\left|i-j\right|-1},\,\,\,\left|i-j\right|>3.$ (29)

Evidently, it is crucial that

$J\delta t/4<1,$ (30)

in order for the resulting Hamiltonian to be truly local, which is but one of
the applicability criteria (10).

Small-distance matrix elements of $\mathcal{H}_{\text{eff}}$ can be computed
in the leading order in powers of $\delta t$ by directly expanding Eqs. (7)
and (8). The exact expression are presented in general-form_for-2nd-order-
correction, while here we present a simple recursive formula for the leading
order defect $\delta
t^{2}\cdot\mathcal{D}=\mathcal{H}_{\text{eff}}-\mathcal{H}$ for the
Hamiltonian of the form $\mathcal{H}=A_{1}+...+A_{n}$:

$\displaystyle\delta\mathcal{\mathcal{D}}_{k+1}-\delta\mathcal{\mathcal{D}}_{k}=$
$\displaystyle-\frac{1}{12}\left\\{\left[\mathcal{H}_{k},\left[\mathcal{H}_{k},A_{k+1}\right]\right]-\frac{1}{2}\left[A_{k+1},\left[A_{k+1},\mathcal{H}_{k}\right]\right]\right\\}$
(31)

where $A_{k}$ are applied according to Eq. (7),
$\mathcal{H}_{k}=\sum_{j=1}^{k}A_{j}$ is the target Hamiltonian after adding
$k$ operators, and the recursion starts from $\mathcal{\mathcal{D}}_{0}=0$.

Straightforward application of the semiclassical description of II.4 is
hindered by the oscillating terms in the effective Hamiltonian. Indeed,
applying the semiclassical description requires that the Wigner’s transform
(19) of the Hamiltonian is a smooth function of $x=an$. However, the
correction $\delta t^{2}\mathcal{D}=\mathcal{H}_{\text{eff}}-\mathcal{H}$ has
the following spatial structure of the matrix elements:

$\mathcal{D}^{ij}=U^{ij}+\left(-1\right)^{i}W^{ij}+\left(-1\right)^{j}\overline{W^{ij}},$
(32)

with both operators $U,W$ having smooth Wigner’s transforms (19). The
Trotterization procedure thus introduces an oscillating term with matrix
elements proportional to $\left(-1\right)^{n}$ that seemingly renders the
semiclassical approximation inapplicable. Upon closer look, the problematic
term causes weak scattering of waves with momentum $p$ to those with momentum
$p\pm\pi$, thus requiring a change of energy by $\sim 2J$. Therefore, all such
processes are virtual and amount to an additive correction to the Hamiltonian
$\mathcal{H}_{\text{eff}}$ with a quantity of the order $\delta
t^{4}W^{2}/J^{2}$. The Wigner’s transform of this correction is yet another
smooth function and constitutes $O\left(\delta t^{4}\right)$ contribution,
which is beyond the precision of the analysis. As a result, we can safely
discard this term in our analysis, as far as semiclassical quantization
condition is concerned. A more detailed version of this argument is presented
in Correction-to-semiclassical-Hamiltonian.

There is additional simplification due to the same assumption of smoothness
that validates the semiclassical approximation. Consider, for instance, the
$\left\\{K_{\text{even}},K_{\text{odd}},P\right\\}$ ordering of the operators
in the Suzuki-Trotter approximation (7). Eq. (31) then suggests that the
corrections contains terms from noncommutativity of $K_{\text{even}}$ and
$K_{\text{odd}}$ as well as terms from noncommutativity of the full kinetic
energy $K=K_{\text{even}}+K_{\text{odd}}$ with the potential term $P$.
However, if the potential is homogeneous, it certainly commutes with $K$, so
the second contribution vanishes. Inhomogeneity of the potential then produces
contributions that are proportional to $a^{2}\left(\nabla
h\right)^{2},a^{2}\Delta h$, so they can safely be neglected as subleading
w.r.t the main semiclassical approximation. A more rigorous analysis for all
possible groupings of operators is presented in Correction-to-semiclassical-
Hamiltonian.

###  III.2 Trotterization as a periodic perturbation at large frequency
$\Omega=2\pi/\delta t$ and prethermalization

The Suzuki-Trotter algorithm (7) is described by the following piecewise-
constant time-dependent Hamiltonian:

$H_{\text{appr.}}\left(t\right)=n\times\begin{cases}A_{n},&\frac{1}{2n}<t/\delta
t\leq\frac{1}{n},\\\ ...,\\\ A_{k},&\frac{n-k}{2n}<t/\delta
t\leq\frac{n-k+1}{2n},\\\ ...,\\\ A_{1},&\frac{n-1}{2n}<t/\delta
t\leq\frac{n+1}{2n},\\\ ...,\\\ A_{n},&\frac{2n-1}{2n}<t/\delta t\leq
1,\end{cases}$ (33)

which is periodically repeated every Trotter step of length $\delta t$, with
the overall factor of $n$ needed simply to identify the “system time” $t$
describing the evolution of the target model (1) with the actual physical time
on the device. Being averaged over one Trotter step, it produces the target
Hamiltonian, i.e.

$\intop_{0}^{\delta t}dt\,H_{\text{appr.}}\left(t\right)=\delta t\,H.$ (34)

The difference $H_{\text{appr.}}\left(t\right)-H$ can be regarded as high-
frequency perturbation, as it can be expanded in cosine series with the main
frequency $\Omega=2\pi/\delta t$ (due to symmetry $t\leftrightarrow\delta t-t$
in $H_{\text{appr.}}\left(t\right)$, sine components are absent).

Because $H_{\text{appr.}}\left(t\right)$ is periodic with period $\delta t$,
the Floquet theorem guarantees the existence of the effective time-independent
Hamiltonian $H_{\text{eff}}$ that describes the evolution of the system at
stroboscopic times $t_{n}=n\delta t,\,\,n\in\mathbb{N}$:
$\psi\left(t_{n}\right)=\exp\left\\{-it_{n}H_{\text{eff}}\right\\}\psi\left(0\right)$.
One can then develop perturbative expansion [22, Sec. 3] to express the
effective Hamiltonian $H_{\text{eff}}$ via a series in powers of $1/\Omega$.
In particular, Eq. (31) describes the leading term of this expansion of order
$\Omega^{-2}$ (with $\Omega^{-1}$ term vanishing due to symmetry
$t\leftrightarrow\delta t-t$ of the circuit.

The question then arises whether the high-frequency expansion is applicable in
the thermodynamic limit. Indeed, one expects a system under periodic drive to
absorb energy, unless the latter is in many-body localized phase. However,
heating is clearly inconsistent with the exact conservation of any quasi-local
Hamiltonian $H_{\text{eff}}$ with finite energy density at stroboscopic times.
This implies that the power series of $H_{\text{eff}}$ must be asymptotic,
i.e. diverge beyond certain order $n_{*}$. However, as it was shown in Ref.
[15], truncating the power series at optimally chosen order $n_{*}$ allows to
correctly describe stroboscopic evolution up to time $\tau_{*}$ that is
exponentially large in frequency, viz., $\ln\tau_{*}\sim\Omega/\Lambda$, where
$\Lambda$ is the typical energy scale of the matrix elements of the
Hamiltonian, e.g. $\Lambda=\max\left\\{J,h\right\\}$ for Hamiltonian (1). We
therefore expect our analysis to be qualitatively valid for sufficiently small
Trotter step $\delta t$ if the simulation is carried out for times $T$ shorter
than the prethermalization time $\tau_{*}$. In particular, the latter has to
be longer than the characteristic tunneling times of interest.

However, for the system considered in the present work the issue of heating is
absent to begin with. The Hamiltonian (1) is mapped to noninteracting
fermions, so the exact dynamics of the many-body wave function is fully
determined by the single-particle evolution operator $\mathcal{U}$, see Eq.
(12). Because $\mathcal{U}$ now describes just one particle with bounded
spectrum, no heating is present, and the associated high-frequency expansion
of the effective single-particle Hamiltonian $\mathcal{H}_{\text{eff}}$ is
thus convergent, allowing one to infer its properties from the perturbation
theory in powers of $\Omega^{-1}$. As a result, our analysis is quantitatively
correct for arbitrary long times for the system (1) in question.

###  III.3 Correction to the semiclassical description and discrete energy
levels

Due to the smoothness of the potential, the procedure of calculating the
correction to the semiclassical description of the problem thus looks as
follows: one applies Eq. (31) only once to compute the contribution to
$H_{\text{eff}}$ arising from noncommutativity of the components
$K_{\text{even}},\,K_{\text{odd}}$ of the kinetic energy. In the resulting
operator, only the part that has smooth dependence on the coordinate is
preserved, and the associated correction to the Wigner’s transform (19) of the
effective Hamiltonian is computed. The result reads:

$\delta\mathcal{H}_{\text{eff}}\left(x,p\right)=J\frac{\left(J\delta
t\right)^{2}}{24}\cos p\,\sin^{2}p.$ (35)

Crucially, this function vanishes at $p=0,\pm\pi$, so the positions of all
classical turning points according to Eq. (18) remain intact. Note also that
this correction is not sensitive to the ordering of the operators in Eq. (7)
(see also 1), as illustrated in detail in (B.3).

We are now in position to describe the correction to the semiclassical
treatment arising from the presence of perturbation. Similarly to the
semiclassical approach itself, it consists of three parts: _i)_ the change of
the action in the semiclassical regions, both allowed and forbidden ones,
_ii)_ the modification of the matching conditions for the regions where the
semiclassical approximation is not applicable, and _iii)_ the shift in the
energy value induced by the perturbation. The first effect is found from the
perturbed version of the classical equation of motion (18) and reads

$\delta
S\left(E;x\right)=-\intop^{x}dx^{\prime}\,\frac{\delta\mathcal{H}_{\text{eff}}\left(x^{\prime},p\left(E;x^{\prime}\right)\right)}{v\left(p\left(E;x^{\prime}\right)\right)},$
(36)

where $p\left(E;x^{\prime}\right)$ is computed from the unperturbed equation
of motion (18). The deformation of the matching conditions at the classical
turning points is actually absent because it only relies on the fact that the
turning point is isolated and spawned by a potential with small but finite
gradient, as discussed in [App] in more detail. The energy shift of a
semiclassical bound state is described by a deformed version of the
quantization condition (22):

$S_{12}\left(E+\delta E\right)-S_{12}\left(E\right)+\delta
S_{12}\left(E\right)=0.$ (37)

where $\delta S_{12}$ is computed as the integral (36) between the two turning
points, similarly to $S_{12}$ in Eq. (22). While the correction $\delta
S_{12}/a$ can be much larger than $\pi$ if the allowed region is large enough,
both $S_{12}$ and its energy derivative are always large in the classically
allowed region, hence the equation admits perturbative solution:

$\delta E=-\frac{2\,\delta S_{12}\left(E\right)}{T_{12}\left(E\right)},$ (38)

where $T_{12}\left(E\right)=2\,\partial S_{12}/\partial E$ is the classical
period of motion in the allowed region. Remarkably, the answer remains correct
up to a numerical prefactor of order unity even for the states that are not
described by the semicalssical expression [21, par. 48]. Using Eq. (35), this
energy shift can be estimated as

$\left|\delta E\right|\apprle\frac{\left(J\delta
t\right)^{2}}{24}\,\frac{4S_{12}}{T_{12}}\leq\frac{\left(J\delta
t\right)^{2}}{24}\,\frac{2n_{\text{cl}}}{\pi}\Delta\left(E\right)$ (39)

with $\Delta$ being the semiclassical level spacing according to Eq. (24),
$S_{12}$ being the unperturbed classical action of the corresponding energy
level, and $n_{\text{cl}}=\left|x_{2}-x_{1}\right|/a$ is the total number of
sites in the classically allowed region, which serves as an upper bound for
$S_{12}/a$, as evident from Eq. (26).

###  III.4 Dynamics of Rabi oscillations are detuned in general case

The key distortion of the Rabi oscillations comes from the shift of the energy
levels within each of the wells. Within the semiclassical approximation, the
detuning parameter $\varepsilon$ in the resonant subspace Hamiltonian (16) is
found as

$\varepsilon_{\text{eff}}=\varepsilon+\left(-\frac{\delta
S_{L}}{T_{L}}+\frac{\delta S_{R}}{T_{R}}\right),$ (40)

where the index $L,R$ describes the characteristics of left and right well,
respectively, according to Eq.(38). If the two energy shifts are different,
the additional detuning in $\varepsilon_{\text{eff }}$ can be estimated as
$\left(J\delta t\right)^{2}n_{\text{cl}}\Delta$ according to Eq. (39), with
$n_{\text{cl}}$ being the total number of sites in the allowed region. On the
other hand, if $\varepsilon_{\text{eff}}$ becomes significantly larger than
the tunneling amplitude $\eta$, the hybridization of the two levels disappears
and no tunneling is observed. Thus, another applicability condition of the
Suzuki-Trotter approach emerges in the form

$\left(J\delta t\right)^{2}\apprle\eta n_{\text{cl}}/\Delta,$ (41)

where $\eta$ is the tunneling amplitude of Eq. (16). Given that the latter is
itself an exponentially small quantity, the simple detuning of the resonant
levels turns out to be the main limitation on the size of Trotter step $\delta
t$ (e.g., it is much more restrictive than the naive ciriteria (10)). In
particular, it allows one to completely neglect all other effects described
below as long as one is tasked with observing the tunneling in the first
place, as those effects are exponentially small if Eq. (41) is satisfied. The
criterion (41) can be viewed as comparison of the energy scale of
Trotterization $\delta t^{-1}$ to the ones defined by tunneling, with the
latter being exponentially smaller than the energy scales $P,J$ of kinetic and
potential terms of the Hamiltonian.

###  III.5 Rabi oscillations with no detuning: enhancement of tunneling

The change in the tunneling amplitude $\eta$ ought to be addressed as well.
According to Eq. (25), the tunnel splitting $\eta_{N}$ of the resonance
subspace is determined by the action in the forbidden region
$S_{B}\left(E_{N}\right)$ at energy $E_{N}$ of the resonant subspace, as well
as the derivative $T_{N}=2\partial S/\partial E$ of the action in the allowed
region $S_{L,R}$$\left(E_{N}\right)$ at the same energy. The correction to
$\eta_{N}$ thus consists of two effects: direct change of the the action $S$
according to (36) and the energy shift (38) of the resonance subspace. The new
value of the tunneling amplitude is given by

$\frac{\eta_{N,\text{eff}}}{\eta_{N}}\approx\exp\left\\{-\frac{1}{a}\left(\delta
S_{B}\left(E_{N}\right)+\frac{\partial S_{B}\left(E_{N}\right)}{\partial
E}\delta E_{N}\right)\right\\},$ (42)

where the change of the subexponential prefactor of Eq. (25) was neglected in
comparison to the change of the exponent due to the smoothness of the
potential expressed as $a\ll 1$, as explained in Change-in-classical-period.
The first term in the exponent of Eq. (42) reflects the change in the
underbarrier action due to perturbation, Eqs. (35)-(36):

$\delta S_{B}=-\frac{\left(J\delta
t\right)^{2}}{24}\,\intop_{x_{L}}^{X_{R}}dx\,\frac{\sinh
2\left|p\left(x\right)\right|}{2},$ (43)

with the purely imaginary momentum $p\left(x\right)$ found from the equation
of motion (18) in the forbidden region. The second term of Eq. (42) stands for
the change of the underbarrier action (26) due to the induced shift in the
energy of the resonance subspace:

$\delta E_{N}=-\frac{\delta S_{L}}{T_{L}}-\frac{\delta S_{R}}{T_{R}}$ (44)

with the same notations as in Eq. (40).

It is now evident that in the presence of detuning of an accidental resonance
the criterion (41) puts a strong limitation on $\delta t$, rendering the
change of $\eta$ negligible. The situation changes drastically if the effects
of detuning are irrelevant, so the main distortion of the evolution is due to
the change in the tunneling amplitude $\eta$, which can then constitute
multiple orders of magnitude, as it will be shown below. Consider the case of
two _symmetric_ wells (2, left), in which case each discrete state in one of
the wells is weakly coupled to its counterpart in the other well. The
introduced perturbation then causes symmetric shifts in the wells, provided
the symmetry of the potential is respected by the perturbation. As a result,
all of the perturbed discrete states remain perfectly hybridized
($\varepsilon\equiv 0$) with the corresponding partner in the other well, but
the tunneling amplitude $\eta$ of each resonant subspace is distorted by
perturbation according to Eq. (42).

Importantly, the corrected value $\eta_{\text{eff}}$ is always larger than
original $\eta$ because the correction to the tunneling action (43) is
negative for the majority of states. Indeed, the first term in the exponent of
Eq. (43) is always negative, as evident from Eq. (43), while the second term
is negative for states with absolute value of classical momentum either always
above or always below $\pi/2$, as can been verified from Eqs. 38 and (26) for
$\delta E_{N}$ and $S_{B}\left(E_{N}\right)$, respectively. This fact
expresses the suppression of the potential barrier by the Trotterization due
to next-to-nearest neighbor hopping.

Furthermore, if $\delta t$ is not limited by the condition (41) due to
detuning, large values of the exponent in Eq. (42) become possible, implying
the change of $\eta_{\text{eff}}$ by several orders of magnitude, as the
effect of next-to-nearest hopping accumulates along the length of the barrier.
Indeed, the first term in the exponent of Eq. (42) is always negative and can
roughly be estimated as

$\frac{\delta S_{B}}{a}\sim-\left(J\delta
t\right)^{2}n_{\text{barr}}\times\begin{cases}\sqrt{P/J},&P/J\ll 1,\\\
\left(P/J\right)^{2},&P/J\gg 1,\end{cases}$ (45)

where $n_{\text{barr}}$ is the number of sites in the forbidden region
corresponding to the barrier, $P=\max_{x\in\text{barr}}\left|h(x)-E\right|-J$
is the energy scale of the potential barrier. While it is difficult to provide
a similar estimation for the second term in the exponent of Eq. (42), it only
enhances the suppression of the potential barrier as it is also negative in
most cases. As a result, if $\left|\delta S_{B}\right|/a\apprge 1$, which
might happen even if $J\delta t\ll 1$ for sufficiently large heights $P/J$ and
lengths $n_{\text{barr.}}$ of the barrier, the effective tunneling amplitude
$\eta_{\text{eff}}$ might become exponentially larger than that of the
unperturbed Hamiltonian.

Of course, the resulting value of $\eta_{\text{eff}}$ still has to be small
compared to all other energy scales of the problem in order for the
semiclassical approximation to be applicable. In the spirit of Eq. (45), the
under-barrier action $S_{B}$ without perturbation can be estimated as

$\frac{S_{B}}{a}\sim n_{\text{barr}}\times\begin{cases}\sqrt{P/J},&P/J\ll
1,\\\ P/J,&P/J\gg 1.\end{cases}$ (46)

Requiring now that $\eta_{\text{eff}}\ll J,P,\Delta,...$ while still assuming
that the barrier is semiclassical, i.e., $n_{\text{barr}}\gg 1$ and
$S_{B}/a\gg 1$, thus leads to the following condition on $\delta t$:

$\frac{S_{B}+\delta S_{B}}{a}\gg 1\Rightarrow\left(J\delta
t\right)^{2}<\begin{cases}1,&P/J\ll 1,\\\ \left(P/J\right)^{-1},&P/J\gg
1,\end{cases}$ (47)

which is literally equivalent to the naive criteria (10). As a result, within
its entire range of applicability, the Suzuki-Trotter approximation is capable
of strongly enhancing tunneling as compared to the unperturbed evolution
provided that the detuning effects described in III.4 are irrelevant.

The case of high barriers $P/J\gg 1$ is additionally interesting in the
following regard: the estimations (45) and (46) show that the additional
tunneling action (43) due to the Suzuki-Trotter approximation grows faster
with the height $P/J$ of the potential barrier than the unperturbed tunneling
action (26), while the ensuing criteria of applicability Eq. (47) shows that
the true parameter of the perturbation theory is $J\delta t\cdot P\delta t$.
Should it be large, the perturbation $\delta\mathcal{H}$ actually governs the
motion of the particle under the barrier, as evident from Eq. (47), implying
that higher order in $\delta t$ corresponding to large-distance hopping also
become important. In certain cases, these terms can be analyzed analytically,
as explained in III.8.

###  III.6 Enhancement of the decay rate

Figure 3: _Left:_ visualization of the portion of the eigenlevels (horizontal
dashed lines) and eigenfunctions (curves with colored filling) for the problem
with $L=200$ sites and potential profile $h_{n}/J=5/4\,\cos 4\pi\,(n-1)/(L-1)$
(thick blue line) realizing the double-well experiment (2, left). Only the
left well is shown, corresponding to $25\leq n\leq 75$. The tunnel splitting
between symmetric and antisymmetric combination of the two wells is too small
to be visible. The apparent forbidden zone at $E/J\geq 0.75$ corresponds to
appearance of anomalous turning points as the quantity $E-h\left(x\right)$
exceeds the band width $2J$. All energies are shifted by $J$. _Right_ : color
plot of the absolute value of the overlap $\left|\left\langle
M|N\right\rangle\right|$ between a given local level $\left|N\right\rangle$ in
the left well of the exact Hamiltonian and the local level
$\left|M\right\rangle$ of the effective Hamiltonian of the left circuit on 1
with $\delta tJ=0.2$ for the system used on the left figure. For $N=M$ the
quantity $\sqrt{1-\left|\left\langle
N|N_{\text{eff}}\right\rangle\right|^{2}}$is plotted. All overlaps are
normalized by $\left(J\delta t\right)^{2}$. Diagonal white dashed lines
correspond to (from bottom to top) _i)_ $E_{M}=E_{N}$, _ii_)
$E_{M}=2J-\left(E_{N}+h\right)$, _iii)_ $E_{M}=2J+\left(E_{N}+h\right)$, so
significant overlap close to these lines corresponds to _i)_ smearing of the
target state over close-lying states, _ii)_ overlap with states that have
typical momenta $\pi-p$, where $p$ is the typical momentum of the original
state, and _iii)_ overlap with states with energy difference $2J$,
corresponding to a different band. Computing the same quantity by 1st order
perturbation theory $\left\langle N|M\right\rangle\approx\left\langle
N|\delta\mathcal{H}|M\right\rangle/\left(E_{N}-E_{M}\right)$ and Eq. (49) for
$N=M$ gives indistinguishable results.

The case of decay of a metastable state is governed by the same physics as the
Rabi oscillations, safe for the different structure of single-particle levels
at either side of the potential barrier (see 2, right). According to Eq. (27),
the decay rate of a given metastable state $\left|N\right\rangle$ inside the
well changes to

$\frac{\Gamma_{N,\text{eff}}}{\Gamma_{N}}\approx\exp\left\\{-\frac{2}{a}\left(\delta
S_{B}\left(E_{N}\right)+\frac{\partial S_{B}\left(E_{N}\right)}{\partial
E}\delta E_{N}\right)\right\\},$ (48)

with the only difference with Eq. (42) being the extra factor of 2 in the
exponent. Similarly to the previous case, the change of prefactor of Eq. (27)
is neglected as a subleading correction, although it can trivially be taken
into account, see Change-in-classical-period.

What is different with the case of Rabi oscillations is the effect of
detuning: because all states inside the well are metastable, detuning does not
prevent the decay process because it cannot detune the resonance with the
dense spectrum on the other side of the well. In other words, the situation is
similar to that of the previous III.5: _i)_ in most cases, trotterization
enhances tunneling because the exponent of Eq. (48) is negative, and _ii)_
this enhancement can constitute several orders of magnitude, but _iii)_ for
sufficiently high barriers the perturbative treatment of III.3 is
inapplicable, necessitating a more accurate description of the underbarrier
motion with large-distance jumps taken into account.

###  III.7 Probability defect

If the particle was originally prepared with the wave function
$\left|N\right\rangle$ of a resonant state inside either of the wells, the
perturbation due to Suzuki-Trotter approximation changes the stationary wave
function $\left|N_{\text{eff}}\right\rangle$, so the latter acquires a finite
overlap with other localized states $\left|M\right\rangle,\,\,M\neq N$ inside
the well. As far as tunneling processes are concerned, of particular interest
is the overlap with high-energy states, as the latter have much smaller
tunneling exponents or even no tunnel barriers at all, thus introducing
completely different character of evolution to the system. Motivated by this,
we estimate the total probability defect $\delta P=1-\left|\left\langle
N_{\text{eff}}|N\right\rangle\right|^{2}$ due to the distortion of the
semiclassical wave function in the allowed region. The detailed calculation is
presented in Probability-defect, and the result reads

$\delta P_{N}\approx\sum_{N^{\prime}\neq N}\frac{\left|\left\langle
N^{\prime}\left|\delta\mathcal{H}\right|N\right\rangle\right|^{2}}{\left(E_{N}-E_{N^{\prime}}\right)^{2}}\approx
C_{N}\left(J\delta t\right)^{4},$ (49)

where $C_{N}$ is a positive constant bounded from above by $n_{\text{cl}}^{2}$
(with $n_{\text{cl}}$ being the number of sites in the classically allowed
region at energy $E_{N}$) but typically of the order of 1, as shown in
Probability-defect. The resulting probability defect is guaranteed to be small
if

$\left(J\delta t\right)^{2}n_{\text{cl}}\ll 1.$ (50)

This result can be understood from energy conservation. Indeed, the typical
shift of the energy of a given state inside the well is given by Eq. (38), and
if this shift approaches the typical level spacing $\Delta$ (which happens
exactly when the criterion (50) is violated) then the initial state has no
choice but to get distributed across $\delta N\sim\delta
E_{N}/\Delta\sim\left(J\delta t\right)^{2}n_{\text{cl}}\sim 1$ neighboring
eigenstates of the distorted system.

It is important, however, that the distribution of the overlaps $\left\langle
N^{\prime}|N_{\text{eff}}\right\rangle$ might be quite different from the one
expected from naive energy conservation law. In fact, two additional effects
arrange a significant overlap with higher-energy states: _i)_ the perturbation
arranges overlap with states at energy difference $\approx J$ if those are
present, as both have comparable classical momenta and _ii)_ the presence of
oscillating term, Eq. (32), in the effective Hamiltonian (neglected up until
this point) arranges a finite overlap between states with momenta $\sim p$ and
$\sim\pm\pi+p$. Those two phenomena are present on Fig. 3, which illustrates
the matrix elements of the perturbation. While not directly relevant to the
single-particle problem, this phenomena might turn out to be important for
many-body problem by means of so-called orthogonality catastrophe [23].

###  III.8 Large Trotter steps within the semiclassical picture

We now discuss the effect of large-distance matrix elements of the full
effective Hamiltonian $\mathcal{H}_{\text{eff}}$, which become important if
$J\delta t$ is not very small, with special emphasis on the forbidden region,
as hinted by the estimations of III.5. The peculiarity of this region can be
seen in the Wigner’s transform (19) of the effective Hamiltonian given the
estimation (29) of the large-distance matrix elements. Indeed, the Wigner’s
transform of the full Hamiltonian can be represented as

$\displaystyle\delta\mathcal{H}\left(x=an,p\right)$
$\displaystyle=\left\langle
n\left|\delta\mathcal{H}\right|n\right\rangle+\sum_{\left|k\right|=1}^{3}\left[\left\langle
n\left|\delta\mathcal{H}\right|n+k\right\rangle e^{ipk}\right]$
$\displaystyle+\sum_{\left|k\right|>3}\left\langle
n\left|\delta\mathcal{H}\right|n+k\right\rangle e^{ipk},$ (51)

where the summation over $k$ includes both positive and negative values.
According to III.1, the first two terms are of the order of $\left(J\delta
t\right)^{2}$ and are discussed in detail in Sections III.1-III.6, while the
last one is $O\left(\left[J\delta t\right]^{3}\right)$ and the sum terms decay
with $k$ exponentially as $\left(J\delta t/4\right)^{k}$, according to Eq.
(29). It then follows that the latter term is not important as far as the
motion with real momentum is concerned. However, the higher-order term
diverges at a imaginary momentum corresponding to

$p_{c}\approx\pm i\ln\frac{4}{J\delta t}.$ (52)

Clearly, the presence of this singularity has important implications for the
tunneling processes, as those are represented by the motion with imaginary
momenta. If the latter approach the threshold value (58), one expects the
semiclassical behavior to be significantly altered. Estimating the maximum
momentum in the forbidden region by its unperturbed dynamics, Eqs. (18-19), we
obtain that the tunneling trajectory does not approach this singularity if

$\left(h(x)-E\right)\delta t\ll 1.$ (53)

In this way, the second part of the criteria of applicability (10) of the
Suzuki-Trotter approximation is reinstated. If condition (53) is violated,
large-distance hopping becomes more favorable than tunneling through the
barrier. Moreover, as evident from Eq. (53), the low-energy subspace is the
first to be affected, with the most prominent effect occurring exactly at the
point with the highest potential barrier. This also elucidates the observed
difference in the scaling of the underbarrier action (46) and the correction
(45) with the potential height $P$: Suzuki-Trotter error is dominant in the
regions with sufficiently high potential barriers.

####  III.8.1 Floquet theorem and practical implications

Figure 4: Visualization of the portion of the eigenlevels (horizontal _black_
dashed lines) and eigenfunctions (curves with colored filling) for the Suzuki-
Trotter simulation of the problem with the left circuit of 1 and trotter step
$J\delta t=1.49707035$, $L=50$ sites and potential profile $h_{n}/J=5/4\,\cos
4\pi\,(n-1)/(L-1)$ (thick blue line) realizing the double-well experiment (2,
left). The horizontal _red_ dashed lines represent the position of the energy
levels of the original Hamiltonian (15). The tunnel splitting between
symmetric and antisymmetric combination of the two wells is too small to be
visible. All energies are shifted by $J$. The precise choice of $\delta t$ is
required to tune the mode inside the forbidden region into resonance with one
of the two states in the two wells, so two instances (green and orange) of
hybridization of the antisymmetric state with the bound state inside the
barrier state are observed. The energy splitting between symmetric (blue) and
any of the the lowest antisymmetric states (green, orange) is $\sim 10^{9}$
bigger than the one in the unperturbed case, while for higher pairs of states
(the first one on the plot is purple and red) the amplification of splitting
is given by $9.3\times 10^{4}$, $4.1\times 10^{3}$. Note that the bound state
deep inside the potential barrier corresponds to momenta close to $\pi$ as its
wave function demonstrates osculations of the form
$\psi_{n}\propto\left(-1\right)^{n}.$

As explained above, when either $J\delta t$ or $P\delta t$ is of order unity
one expects certain qualitative differences between the evolution under the
original Hamiltonian (15) and that under the corresponding Suzuki-Trotter
approximation. Indeed, as the total bandwidth of the original Hamiltonian
$\Lambda\approx 2J+P$ approaches the value of $\Omega=2\pi/\delta t$,
inelastic processes involving changes of energy by a multiple of $\Omega$
cease to be purely virtual, thus causing a complete reconstruction of the
eigenstates of the evolution. Mathematically, this is expressed by the Floquet
theorem: because the energy spectrum of a system evolving under periodic drive
with period $\delta t$ is defined modulo $\Omega=2\pi/\delta t$, high-energy
states of the original Hamiltonian get mixed with low-energy ones due to the
emerging periodicity in energy. The latter effect causes a strong distortion
of the evolution in the low-energy subspace of interest.

One particular example of such distortion is shown on 4: the system is find-
tuned close to $J\delta t\sim h\delta t\sim\pi/2$ to bring the high-energy
states localized _inside_ the barrier into resonance with bound states of a
potential well. The resulting tunneling amplitude is exponentially enhanced
due to the presence of an intermediate state.

One important practical conclusion also follows from the Floquet theorem:
Suzuki-Trotter approximation is a poor choice if one is interested in
simulating quadratic dispersion law $\varepsilon_{p}\approx Jp^{2}/2$
corresponding to $J\gg E+J\sim P$. Indeed, the criteria of applicability (10)
then prescribe one to use Trotter steps of size $\delta t\ll J^{-1}$, thus
giving no benefit from small energy scales $E+J\sim P\ll J$. Conversely, as
soon as one reaches $J\delta t\sim 1$ the target low-energy states are
rendered strongly distorted by the Suzuki-Trotter approximation, as hinted by
Eq. (49) for the total probability defect and promptly visible on 4. The
quantitative description of this distortion is given in the next subsection.

####  III.8.2 Semiclassical description

The situation of large Trotter steps $\left(2J+P\right)\delta t\sim 1$ can be
understood semiclassically. Indeed, the semiclassical approximation
essentially relies on the fact that the classical momentum $p$ changes slowly
with spatial coordinate due to smoothness of the potential profile. On the
other hand, the latter also manifests itself in the smallness of corrections
(31) to the effective Hamiltonian due to the noncommutativity of the kinetic
and potential terms. The main correction thus originates from noncommutativity
of the kinetic terms on even and odd sites _in all orders_ in powers of
$J\delta t$, as the latter is not small anymore. However, the procedure of
collecting all those corrections exactly corresponds to computing the Suzuki-
Trotter evolution operator (7) in absence of the potential term, which is
straightforward in the momentum space. The semiclassical description of the
Suzuki-Trotter evolution in a smooth potential thus amounts to diagonalizing
the effective Hamiltonian in absence of the potential term and using the
resulting wave functions to perform the semiclassical ansatz. The detailed
procedure is presented in Full-effective-Hamiltonian, and here we just quote
the result.

The semiclassical Hamiltonian is given by

$\mathcal{H}_{\text{eff}}\left(x,p\right)=-\frac{2}{\delta
t}\arcsin\left\\{\sin\frac{J\delta t}{2}\,\cos p\right\\}+h(x),$ (54)

where we assume that $J\delta t<2\pi$. Naturally, expanding this expression in
powers of $\delta t$ restores the unperturbed Hamiltonian (19) and the
correction (35). On the other hand, the form of the kinetic energy coincides
with that of Ref [1] for $J\delta t=\pi/2$.

$\psi\left(x=na\right)=\frac{e^{iS(x)/a}}{\sqrt{v_{p}}}\left(\cos\theta_{p}-ie^{i\pi
x/a}\,\sin\theta_{p}\right),$ (55)

where $p=S^{\prime}(x)$, the angle $\theta_{p}\in\left[0,\pi\right]$ is
determined by the problem without potential:

$\sin 2\theta_{p}=\frac{\sin^{2}\frac{J\delta t}{4}\,\sin
2p}{\sqrt{1-\left(\cos p\sin\frac{J\delta t}{2}\right)^{2}}},$ (56)

and the normalization $v$ is the group velocity given by

$v_{p}=\frac{\partial T}{\partial p}\,\cos
4\theta_{p}+\frac{1}{2}T\,\frac{\partial\cos 4\theta_{p}}{\partial p},$ (57)

with $T(p)$ being the kinetic energy of the semiclassical Hamiltonian 54.

Crucially, the amplitude of wave function (55) now contains both a smooth
contribution and strongly oscillating contribution proportional to $(-1)^{n}$,
which is the manifestation of the fact that the true translation period of the
system under the Suzuki-Trotter evolution operator is doubled. It is for the
same reason that the group velocity (57) differs from its original value
$J\sin p$. For small $J\delta t$, the difference of the new wave function (55)
and its counterpart (17) in the original system is $O\left(\left(J\delta
t\right)^{2}\right)$, in agreement with considerations of III.1, while the
correction to the group velocity $v$ is $O\left(\left(J\delta
t\right)^{4}\right)$, which is beyond the precision of the previous
subsections and hence also absent therein. On the contrary, for $J\delta
t\sim\pi/2$ the magnitudes of the two terms in the prefactor of the wave
function become comparable even for small momenta
$\left|p\right|\apprge\left|\pi/2-J\delta t\right|$, indicating the presence
of strong oscillations even in the wave functions at low energies. At $J\delta
t\sim\pi$, $\theta_{p}$ approaches $\pi/2$, rendering the oscillating
component the dominant one.

Figure 5: Phase portraits of the effective classical Hamiltonian (54) with the
same parameters as on 4 for three lowest (solid lines) and three highest
(dashed lines) eigenenergies of the unperturbed problem. Above the main figure
presented is the corresponding potential profile. The hybridization observed
on 4 appears due to the exact resonance between the energies of the bound
states inside the three regions around $x/L=1/4,\,1/2,\,3/4$ (expressed as a
match in the areas enclosed by the corresponding classical trajectories, in
accordance with the quantization condition). This also explains huge
amplification of the tunnel splitting quoted on 4.

In what concerns the motion with imaginary momenta corresponding to tunneling,
the kinetic term of (54) has singularity at imaginary momenta

$p_{c}=\pm i\text{arccosh}\,\frac{1}{\sin J\delta t/2}$ (58)

which reduces to Eq. (52) in the limit of small $\delta t$. The associated
equation of motion corresponding to the classical Hamiltonian (54) reads

$\sin\frac{J\delta t}{2}\,\cos p=-\sin\frac{\left(E-h(x)\right)\delta t}{2}.$
(59)

The violation of condition (53) by approaching to the singularity (58) in the
momentum space then corresponds to the appearance of an allowed region inside
the barrier, as shown on 5. This region hosts bound states, which are the same
states that would otherwise be observed at _high_ energies $E\sim 2h+J$. In
particular, this implies that at certain configurations a resonance is
possible between such states and the states inside the well. If, in addition,
$J\delta t\sim\pi/2$, both wave functions exhibit noticeable oscillations and
thus have little overlap with the original bound state inside the well, as
also hinted by the perturbative estimation (49) for the total probability
defect.

##  IV Discussion

###  IV.1 Main qualitative results

We constructed an analytical description of the distortion introduced by the
Suzuki-Trotter algorithm (7) to the evolution of a single particle obeying the
time-independent hopping Hamiltonian (15) with position-dependent potential,
with emphasis on tunneling processes.

Our approach is based on constructing an approximation to the effective time-
independent Hamiltonian $H_{\text{eff}}$ that governs the evolution of the
system between discrete moments of time $t_{n}=n\delta t$, with $\delta t$
being the Trotter step. We do so by examining formal expansion of
$H_{\text{eff}}$ in powers of $\delta t$. While in general this procedure is
hindered by heating due to the periodic drive, the target Hamiltonian (15) can
be mapped to noninteracting fermions, which eliminates the issue of heating,
as explained in III.3. As a result, the evolution under the Suzuki-Trotter
algorithm is fully characterized by an effective time-independent single-
particle Hamiltonian $\mathcal{H}_{\text{eff}}$ that can further be found as a
series in powers of $\delta t$. Due to discrete nature of the evolution, the
eigenenergies of this effective Hamiltonian are defined up to a multiple of
$\Omega=2\pi/\delta t$. That is why the criteria of applicability based on
naive comparison of energy scales of the target Hamiltonian with $\Omega$
accurately capture the magnitude of the evolution distortion. Indeed, the fact
that $\Omega$ is the largest energy scale in the problem simply implies that
all transitions induced by the Suzuki-Trotter error are virtual.

The situation of time-independent Hamiltonian and the emergent periodic
evolution is to be contrasted to the case of aperiodic time-dependent
Hamiltonian, for which low-frequency tails in the continuous spectrum of the
perturbation might lead to significant modification of the tunneling dynamics
by cascades of inelastic transitions [24].

Under the assumption of smooth potential profile, we employed the
semiclassical approximation and shown that the main contribution to the error
comes from noncommutativity of kinetic terms at neighboring pairs of sites.
This contribution produces an additional term to the semiclassical
Hamiltonian:

$\delta\mathcal{H}_{\text{eff}}\left(x,p\right)=J\frac{\left(J\delta
t\right)^{2}}{24}\cos p\,\sin^{2}p,$ (60)

which, in turn, leads to quantitative changes in certain tunneling processes.

One of the main effects of Trotterization is the shift of energies caused by
the perturbation in the classically allowed region. The magnitude of this
shift is given by (38) and can roughly be estimated as

$\left|\delta E\right|\apprle\frac{\left(J\delta
t\right)^{2}}{24}\,\frac{2n_{\text{cl}}}{\pi}\Delta,$ (61)

where $J$ is the characteristic kinetic energy scale, $\Delta$ is the relevant
level spacing, $n_{\text{cl}}$ is the size of the classically allowed region
(that is, the region of localization of the metasatable state that
participates in tunneling). This shift has the potential to detune the
resonance and thus render the tunneling impossible in the first place. This
effect puts a limitation on the Trotter size $\delta t$:

$\eta\leq\delta E\Leftrightarrow\left(J\delta t\right)^{2}\apprle\eta
n_{\text{cl}}/\Delta,$ (62)

where $\eta$ is the tunneling amplitude, an exponentially small quantity,
which makes Eq. (62) the most essential limitation on the size of the Trotter
step $\delta t$.

Hand in hand with the detuning goes the defect in the overlap between the
initial state and the target subspace. If the initial state is prepared to
exactly match the wave function corresponding to a metastable state that
experiences tunneling, the perturbation of the classically allowed region
deforms this target state. The corresponding overlap defect can be estimated
as

$\delta
P_{N}=1-\left|\left\langle\psi_{\text{init}}|\psi_{N}\right\rangle\right|^{2}\sim
C_{N}\left(J\delta t\right)^{4},$ (63)

where $C_{N}\leq n_{\text{cl}}^{2}$ is state-dependent constant that can be
expressed semiclassically, and $n_{\text{cl}}$ is the size of the allowed
region for the $N$-th level. Notably, this probability defect gets distributed
across all kinds of states, including high-energy states, as demonstrated on
3, although the majority of the probability defect is still localized around
the target energy level.

The situation changes significantly if the resonance is not spoiled, e.g., due
to the symmetry of the wells (in which cases it has to be preserved by the
simulation protocol), or continuous spectrum with level spacing
$\delta\ll\eta$ outside the well. In this case, the main effect of the
perturbation is to alter the tunneling amplitude, and the decay rate,
respectively:

$\frac{\eta_{N,\text{eff}}}{\eta_{N}}=\sqrt{\frac{\Gamma_{N,\text{eff}}}{\Gamma_{N}}}\approx\exp\left\\{-\frac{1}{a}\left(\delta
S_{B}+\frac{\partial S_{B}}{\partial E}\delta E\right)\right\\},$ (64)

where $\delta S_{B}$ is the correction to the underbarrier action given by Eq.
(43), $S_{B}$ is the original tunneling action given by Eq. (26), and $\delta
E\sim\left(J\delta t\right)^{2}n_{\text{cl}}\Delta$ is the energy shift
induced by the perturbation. In most cases, the exponent of Eq. (64) is
positive, and the resulting effect can thus represent and exponentially large
_increase_ in the tunneling rate, even for small values of $\delta t$.

An increase in performance due to discretization has been reported for Quantum
Annealing Algorithms [25, 26] that supposedly rely on tunneling. While the
systems used to study quantum optimization algorithms typically feature
strongly disordered potential landscapes, rendering our approach inapplicable
even in principle, one might still expect that the discussed large-distance
hops induced by the Suzuki-Trotter approximation do increase tunneling rates.

We also demonstrated that for sufficiently large values of $\delta t$ the
unitary evolution under the Suzuki-Trotter approximation is starkly different
from the original Hamiltonian evolution. As it is shown in III.8, large-
distance hops dominate the behavior of the system, with the most prominent
effect observed for the states with the largest potential barriers $P\delta
t\sim 1$, where $P\sim\max_{x}h(x)-E$ is the height of the potential barrier.
This also provides an evidence against the use of large-step Suzuki-Trotter
approximation for examining low-energy physics. The values of $\delta t$ that
do not cause such behavior are given by the naive comparison of energy scales:

$J\delta t,P\delta t\ll 1,$ (65)

where $J,P$ is the typical magnitude of the kinetic and potential terms,
respectively.

One practical conclusion follows straightforwardly from the present analysis:
the Suzuki-Trotter algorithm is not suitable for simulating Hamiltonians in
the range of energies that corresponds to quadratic dispersion law
$\varepsilon_{p}\approx Jp^{2}/2,\,\,\,p\ll 1$. Indeed, the latter is achieved
for small energies $E+J\ll J$, while the Suzuki-Trotter error is sensitive to
large energy scales $\sim J$ and influences equally strongly the evolution at
all energies (as suggested both by Eq. (63) and III.8), so the natural
division of scales of the original problem is not respected.

###  IV.2 Relation with the rigorous error bounds

It is instructive to compare our results with the available rigorous bounds
and see whether they are are saturated. The rigorous upper bound on the
spectral norm distance
$\left\|O\right\|=\sqrt{\frac{1}{2^{n}}\text{Tr}O^{\dagger}O}$ (with $n$ being
the number of qubits) between the exact $U$ and approximate $U_{\text{appr}}$
evolution operators during one Trotter step is given by [6, Sec. F.2.1, Eq.
(86)]:

$\left\|U\left(\delta t\right)-U_{\text{appr}}\left(\delta
t\right)\right\|\leq C\frac{x^{3}}{n}+\frac{x^{4}}{12}e^{x}$ (66)

where $x=3n\Lambda\delta t$, $C<1/4$ is certain numerical constant,
$\Lambda=\max\left\\{h,J/2\right\\}$ is the maximum magnitude of local terms
in the Hamiltonian, and $n$ is the system size. This implies a rigorous upper
bound on eigenvalue difference:

$\left|E_{N}-E_{N,\text{eff}}\right|\leq 3C\Lambda\left(3n\Lambda\delta
t\right)^{2},$ (67)

where we neglected the subleading second term in Eq. (66). In comparison to
the estimate of Eq. (39), the upper bound is tight as a function of $\delta t$
and $\Lambda=\text{max}\left\\{P,J/2\right\\}$, while the overall energy scale
$n_{\text{cl}}\Delta$ is replaced by its model-blind upper bound $n\Lambda$ up
to an overall numerical coefficient.

However, the rigorous bound (66) predictably fails to accurately reflect the
structure of error at large times. After $r=t/\delta t$ Trotter steps, the
upper bound on the error in the evolution operator reads [6, Sec. F.2.1, Eq.
(86)]

$\left\|U\left(t\right)-\left[U_{\text{appr}}\left(\delta
t\right)\right]^{t/\delta t}\right\|\leq\frac{t}{\delta
t}\left(C\frac{x^{3}}{n}+\frac{x^{4}}{12}e^{x}\right),$ (68)

which implies that at time scales of the order of inverse detuning (67):

$\Lambda t_{\text{bound}}\sim\frac{\Lambda}{\delta E}\sim\left(3n\Lambda\delta
t\right)^{-2},$ (69)

the difference between the exact and approximate evolution operators can be
expected to attain order of unity, implying a noticeable divergence of the
simulation from the target. However, this statement is strongly dependent on
the physical observable in question and the required precision. As an example,
consider the expectation value of the number of particles in the left well
$\mathcal{N}_{\text{left}}=2\sum_{i=1}^{L/2}S_{i}$ for the double-well
experiment described in II.3 and measurement precision $\epsilon\leq 1$ that
is still larger than the total probability defect, Eq (63). Each pair of
resonant states in left and right well is then characterized by its respective
effective Hamiltonian (16) in the basis of states localized in left and right
well. The initial state is prepared as one of the exact left states
$\left|\psi_{N,\text{left}}\right\rangle$. In the simulation by means of
Suzuki-Trotter algorithm, the time dynamics of the expectation value
$\left\langle\mathcal{N}_{\text{left}}\right\rangle$ is approximately
described as

$\left\langle\mathcal{N}_{\text{left}}\left(t\right)\right\rangle\approx\left(\frac{\varepsilon_{N,\text{eff}}/2}{\Omega_{N,\text{eff}}}\right)^{2}+\left\\{1-\left(\frac{\varepsilon_{N,\text{eff}}/2}{\Omega_{N,\text{eff}}}\right)^{2}\right\\}\cos^{2}\Omega_{N,\text{eff}}t,$
(70)

where
$\Omega_{N,\text{eff}}=\sqrt{\eta_{N,\text{eff}}^{2}+\varepsilon_{N,\text{eff}}^{2}/4}$
, and both the exponentially small contribution of states in the right well to
$\left\langle\mathcal{N}_{\text{left}}\right\rangle$ and the total probability
defect, Eq. (63), were neglected, as the allowed error $\epsilon$is assumed to
be larger than those effects. The dynamics of
$\left\langle\mathcal{N}_{\text{left}}\right\rangle$ is thus starkly different
depending on whether the detuning effect due to Suzuki-Trotter error is
present or not, although in both cases the actual error reveals itself only at
times much larger than $t_{\text{bound}}$. In the former case, the tunneling
is destroyed as $\varepsilon_{N,\text{eff}}\gg\eta_{N}$, so the difference
between now nearly constant
$\left\langle\mathcal{N}_{\text{left}}\right\rangle$ with its unperturbed
counterpart exhibiting Rabi oscillations reaches order of $\epsilon$ at
exponentially large times $t\sim\eta_{N}^{-1}\sqrt{\epsilon}\gg
t_{\text{bound}}$. In the case of no detuning, Rabi oscillations are still
observed in the perturbed evolution, but the frequency of oscillations is
changed according to Eq (64), so the difference of order $\epsilon$ is reached
at
$t\sim\sqrt{\left|\eta_{N,\text{eff}}^{-2}-\eta_{N}^{-2}\right|}\sqrt{\epsilon}\gg
t_{\text{bound}}$. Certainly, this analysis looses the short-time error
occurring due to both detuning (if present) and probability defect $\delta
P_{N}$, which is precisely the content of the rigorous bound (66). However, it
clearly illustrates that rigorous upper bounds overstate the actual magnitude
of error at large times if tunneling is involved.

###  IV.3 Immediate experimental realization

In this section, we propose a particular circuit that demonstrates the
presented physics and its extensions with circuit depth of $1.2\cdot 10^{5}$,
with all the details presented in Experimental-design While Ref. [1] performed
spectral analysis of noisy circuits with depths up to $1.4\cdot 10^{3}$ for a
similar experiment, we expect depths $\sim 10^{5}$ to be possible in the near
future. The latter anticipation is motivated by two reasons: _i)_ the problem
is manifestly single-particle, which renders error mitigation by post-
selection efficient [1], and _ii)_ the coherent times of superconducting
devices are already in the range of $3\times 10^{2}\,\mu s$ [4, 27], more than
10 times larger than that of the device in Ref. [1], and further improvement
is likely to occur. Moreover, due to optimized design choices, the proposed
experiment can tolerate relative static coherent error of order $5\times
10^{-3}$, which is well within the capabilities of the device used in Ref.
[1]. As a result, we expect this experiment to be feasible in the near
feature.

The idea behind the proposed design is to observe the Rabi oscillations in the
double-well potential. One prepares the initial state in one of the wells, and
then observes the occupancy number in either of the wells. To better resolve
the oscillations, one can analyze the Fourier transform of the occupancy
number, as done in Ref. [1]. Due to the limitations of the current devices,
several optimizations have to be made: _i)_ the energy of the initial state
has to be close to the top of the potential barrier, to avoid too large
tunneling times; _ii)_ the typical momentum of the initial state has to be
close to $\pi/2$, in order to reduce the perturbative energy shift (38) that
may otherwise lift the target state above the barrier; _iii)_ the initial
state has to be close to one of the resonant states to maximize the visibility
of the Rabi oscillations. The detailed discussion of these requirements and
the resulting experimental design for a system of $L=50$ qubits are presented
in Experimental-design.

The plot of the squared Fourier harmonic of the occupancy of the left well
$\left|n_{\omega}\right|^{2}$ for various values of the Trotter step $\delta
t$ is presented on 6. The Rabi oscillations are clearly resolved, and the
corresponding period is well described by the semiclassical approximation. In
particular, one can clearly see fivefold acceleration of tunneling at large
Trotter step (the value of $2\pi J/\omega$ drops from $5\times 10^{3}$ to
$1\times 10^{3}$) while the energy shift of the corresponding state is
practically absent (it is smaller than the tunnel splitting at large Trotter
steps).

Figure 6: Density map of the squared absolute value of the Fourier harmonic
$\left|n_{\omega}\right|^{2}$ of the occupancy number in the left well (170)
as a function of period $T=2\pi/\omega$ (with logarithmic scale) and Trotter
step $\delta t$. The details of the simulation are presented in Experimental-
design. Red curves correspond to semiclassical description of the Rabi
oscillations, as detailed in Experimental-design. The remaining discrepancy
between the semiclassical description and direct simulation is due to
corrections beyond the semiclassical description (the potential is not
sufficiently smooth).

###  IV.4 Outlook

Recent literature demonstrates a revival of the debate on the exact
interpretation of what happens in the classically forbidden region during a
tunneling event, with particular emphasize on the question of the tunneling
time [28]. The discussion is complicated by nontrivial interpretation of the
measurement outcomes in many experimental designs. Based on the result of our
analysis, we believe that the existing superconducting qubit systems provide a
promising platform for an in-depth study of the issue, especially given the
unique opportunity to perform simultaneous measurements at every spatial
point. However, a practical challenge of characterizing tunneling
experimentally is to design a clock to measure the time spent under the
barrier, with a notable implementation by a controlled periodic modulation of
the barrier height [29]. The latter is straighforwardly realized in existing
qubit platforms, hence the experiment described in the previous subsection can
be tailored to contribute to the debate on the tunneling time. However, the
structure of error due to the Suzuki-Trotter approximation applied to the
resulting _time-dependent_ Hamiltonian should also be carefully analyzed, as
described below.

The present theoretical analysis can be generalized to higher-order product
formulas that are expected to yield $O\left(\left(J\delta
t\right)^{2k}\right)$ naive accuracy of $H_{\text{eff}}-H$ in powers of
$\delta t$ [11]. Those corrections are clearly more nonlocal and are thus more
sensitive to the exact shape of the potential (expansion in gradients of the
potential converges slower). However, we expect the qualitative content of the
analysis to remain the same as that of $k=1$: effects in the classically
allowed region are perturbative, while the effects on tunneling can be
described by the leading-order correction to the classical tunneling action,
as long as the barrier is not too high. It also follows that the tunneling
events still impose the limitation $P\delta t\ll 1$ regardless of the order
$k$ of the product formula due to the Floquet theorem reasoning (in terms of
III.8, the singularities in the kinetic part of the effective semiclassical
Hamiltonian are still located at $p\sim\pm i\ln\frac{4}{J\delta t}$).

There are many ways to extend the proposed model to more physically nontrivial
situations. One particularly important question is the interplay between
fermionic statistics and Trotterization, which appears as soon as many-body
states are considered. One immediate example is the decay of large spin
density (or, equivalently, fermionic particle density) confined to a local
minimum of the potential such as on 2, right.

Connected to the many-body physics is the question of interaction: our model
enjoys a mapping to non-interacting fermions, whereas all known applications
of the existing quantum hardware necessarily imply true many-body dynamics
that is capable of taking advantage of large Hilbert space [12, 13, 14].
Analyzing the qualitative effects of Trotterization in this class of problems
is thus of great practical importance. One particularly important effect
completely absent in a single-particle setting is that generic many-body
system subject to periodic perturbation (such as the one induced by the
Suzuki-Trotter approximation) will heat up even at high frequencies of the
perturbation [15].

An instance of non-integrable many-body evolution in _homogeneous_ potential
was studied numerically in Ref. [30]. In this work, considered is the Ising
spin chain, with $Z$ being the quantization direction, supplemented by a
certain local interaction term subject to external homogeneous field along $X$
and $Z$ directions, with the $Z$ component being small. By using a carefully
crafted state featuring a localized boundary between regions with opposite
spin polarizations, the authors of Ref. [30] quantitatively analyzed _i)_ the
Bloch’s oscillations of the boundary due to finite band width, and _ii)_
excitation production due to collision of two such boundaries. The question of
whether Trotterization has a qualitative effect on these results then becomes
practically important as soon as one wishes to reproduce and extend these
results using quantum hardware. However, addressing many-body physics of this
type is not possible with existing noisy quantum devices [30] and requires
instead implementing quantum error correction protocols [31, 32].

Several immediate developments of the proposed model concern the introduction
of potential disorder with typical energy scale $W$ to the single-particle
Hamiltonian (15), as has been done in Ref. [1] for small systems. One the one
hand, arbitrarily small disorder induces Anderson localization in one
dimension at length scales comparable to the mean free path [33]. While our
method is not straightforwardly applicable to local potential disorder, we
still expect that the key effect of Suzuki-Trotter approximation for small
disorder amounts to smearing the potential disorder, causing only a
renormalization of the mean free path and, consequently, of the localization
length. On the other hand, sufficiently strong disorder drives the system deep
into the localized phase, which is obviously beyond the approximation of
smooth potential employed in the present work. Such situations are much more
important for the applications [12], and our semiclassical treatment is
clearly inapplicable for the case of strong disorder, so we leave this for
future studies.

Extending the presented analysis to the case of time-dependent Hamiltonian
represents another interesting challenge. While rigorous upper bounds on the
Suzuki-Trotter error for time-dependent Hamiltonians are also available [9],
the question of their saturation is much less trivial, as it depends on the
exact dynamics arranged by the target time-dependent Hamiltonian. In
particular, explicit time dependence renders the Floquet theorem inapplicable,
enabling a wide range of additional physical effects. One particular example
is the dynamic inelastic enhancement of tunneling due to a periodic drive of
the Hamiltonian with frequencies much larger than the the tunneling rate
$\eta$ but smaller than the height of the barrier potential barriers [24].

Finally, an important issue is that of external time-dependent noise, both as
a part of the Hamiltonian and as an inescapable consequence of coupling
between the physical device implementing the qubit system and the environment.
Among the most substantial effects for tunneling is the suppression of
tunneling due to coupling to external bath [34, Sec. 3.5]. It is then of
practical importance to find out whether sufficiently strong large-distance
hopping or other effects induced by Trotterization are capable of influencing
this effect.

###### Acknowledgements.

The authors would like to thank Mikhail V. Feigel’man for numerous fruitful
discussions. A.V.K. is also grateful for the support by Laboratoire
d’excellence LANEF in Grenoble (ANR-10-LABX-51-01).

## References

  * Neill  _et al._ [2021] C. Neill, T. McCourt, X. Mi, Z. Jiang, M. Y. Niu, W. Mruczkiewicz, I. Aleiner, F. Arute, K. Arya, J. Atalaya, _et al._ , Nature 594, 508 (2021).
  * Georgescu  _et al._ [2014] I. M. Georgescu, S. Ashhab, and F. Nori, Rev. Mod. Phys. 86, 153 (2014).
  * Arute  _et al._ [2019] F. Arute, K. Arya, R. Babbush, D. Bacon, J. C. Bardin, R. Barends, R. Biswas, S. Boixo, F. G. S. L. Brandao, D. A. Buell, _et al._ , Nature 574, 505–510 (2019).
  * Place  _et al._ [2021] A. P. M. Place, L. V. H. Rodgers, P. Mundada, B. M. Smitham, M. Fitzpatrick, Z. Leng, A. Premkumar, J. Bryon, A. Vrajitoarea, S. Sussman, _et al._ , Nature communications 12, 1779 (2021).
  * Satzinger  _et al._ [2021] K. J. Satzinger, Y.-J. Liu, A. Smith, C. Knapp, M. Newman, C. Jones, Z. Chen, C. Quintana, X. Mi, A. Dunsworth, _et al._ , Science 374, 1237–1241 (2021).
  * Childs  _et al._ [2018] A. Childs, D. Maslov, T. Nam, N. J. Ross, and Y. Su, Proceedings of the National Academy of Sciences 115, 9456–9461 (2018).
  * Berry  _et al._ [2007] D. W. Berry, G. Ahokas, R. Cleve, and B. C. Sanders, Communications in Mathematical Physics 270, 359–371 (2007).
  * Haah  _et al._ [2018] J. Haah, M. B. Hastings, R. Kothari, and G. H. Low, SIAM Journal on Computing Special Section FOCS, 250 (2018).
  * Childs and Su [2019] A. M. Childs and Y. Su, Physical review letters 123, 050503 (2019).
  * Childs  _et al._ [2019] A. M. Childs, A. Ostrander, and Y. Su, Quantum 3, 182 (2019).
  * Suzuki [1991] M. Suzuki, Journal of Mathematical Physics 32, 400–407 (1991).
  * Smelyanskiy  _et al._ [2020] V. N. Smelyanskiy, K. Kechedzhi, S. Boixo, S. V. Isakov, H. Neven, and B. Altshuler, Physical Review X 10, 011017 (2020).
  * Albash and Lidar [2018] T. Albash and D. A. Lidar, Rev. Mod. Phys. 90, 015002 (2018).
  * Bharti  _et al._ [2022] K. Bharti, A. Cervera-Lierta, T. Kyaw, T. Haug, S. Alperin-Lea, A. Anand, M. Degroote, H. Heimonen, J. Kottmann, T. Menke, _et al._ , Rev. Mod. Phys. 94, 015004 (2022).
  * Abanin  _et al._ [2017] D. A. Abanin, W. De Roeck, W. W. Ho, and F. Huveneers, Physical Review B 95, 014112 (2017).
  * Lieb  _et al._ [1961] E. Lieb, T. Schultz, and D. Mattis, Annals of Physics 16, 407–466 (1961).
  * Yoshida [1990] H. Yoshida, Physics letters A 150, 262–268 (1990).
  * Berry and Mount [1972] M. V. Berry and K. E. Mount, Reports on Progress in Physics 35, 315 (1972).
  * Braun [1993] P. A. Braun, Reviews of modern physics 65, 115 (1993).
  * Bloch [1929] F. Bloch, Zeitschrift für physik 52, 555–600 (1929).
  * Landau and Lifshitz [2013] L. D. Landau and E. M. Lifshitz,  _Quantum mechanics: non-relativistic theory_ , Vol. 3 (Elsevier, 2013).
  * Bukov  _et al._ [2015] M. Bukov, L. D’Alessio, and A. Polkovnikov, Advances in Physics 64, 139–226 (2015).
  * Anderson [1967] P. W. Anderson, Physical Review Letters 18, 1049 (1967).
  * Ivlev and Melnikov [1992] B. I. Ivlev and V. I. Melnikov, in  _Quantum tunnelling in condensed media._ , edited by Y. Kagan and A. J. Leggett (Elsevier, 1992) Chap. 5, p. 265–312.
  * Heim  _et al._ [2015] B. Heim, T. F. Rønnow, S. V. Isakov, and M. Troyer, Science 348, 215–217 (2015).
  * Steiger  _et al._ [2015] D. S. Steiger, T. F. Rønnow, and M. Troyer, Physical review letters 115, 230501 (2015).
  * Sivak  _et al._ [2023] V. V. Sivak, A. Eickbusch, B. Royer, S. Singh, I. Tsioutsios, S. Ganjam, A. Miano, B. L. Brock, A. Z. Ding, L. Frunzio, _et al._ , Nature 616, 50 (2023).
  * Sainadh  _et al._ [2020] U. S. Sainadh, R. T. Sang, and I. V. Litvinyuk, Journal of Physics: Photonics 2, 042002 (2020).
  * Büttiker [1990] M. Büttiker, Traversal, Reflection and Dwell Time for Quantum Tunneling, in  _Electronic Properties of Multilayers and Low-Dimensional Semiconductor Structures_, edited by J. M. Chamberlain, L. Eaves, and J.-C. Portal (Springer US, Boston, MA, 1990) pp. 297–315.
  * Milsted  _et al._ [2022] A. Milsted, J. Liu, J. Preskill, and G. Vidal, PRX Quantum 3, 020316 (2022).
  * Google Quantum AI [2023] Google Quantum AI, Nature 614, 676 (2023).
  * Kitaev [2003] A. Kitaev, Annals of Physics 303, 2–30 (2003).
  * Abrahams  _et al._ [1979] E. Abrahams, P. W. Anderson, D. C. Licciardello, and T. V. Ramakrishnan, Phys. Rev. Lett. 42, 673 (1979).
  * Kamenev [2023] A. Kamenev,  _Field theory of non-equilibrium systems_ (Cambridge University Press, 2023).
  * Maslov and Fedoriuk [2001] V. P. Maslov and M. V. Fedoriuk,  _Semi-classical approximation in quantum mechanics_ , Vol. 7 (Springer Science & Business Media, 2001).

##  Appendix A Discrete Semiclassical Approximation

In this Appendix, the semiclassical approximation for 1D systems is reviewed.
The derivation of the semiclassical ansatz is in line with Ref. [35], and the
matching procedure is done similarly to Ref. [18].

###  A.1 Semiclassical wave function

We start from the exact Shroedinger equation:

$\sum_{m}\mathcal{H}_{nm}\psi_{m}=E\psi_{n},$ (71)

where $\mathcal{H}_{nm}$ are matrix elements of the Hamiltonian in certain
basis, such as the site position basis for the Hamiltonian (15) of the present
paper. We assume that in this basis the matrix elements $\mathcal{H}_{nm}$
obey two conditions: _i)_  _locality_ : $\mathcal{H}_{nm}$ decay quickly
(e.g., exponentially) with the difference $\left|n-m\right|$, and _ii)_
_smoothness_ : $\mathcal{H}_{nm}$ depend smoothly on $a\left(n+m\right)$ with
some characteristic scale $a\ll 1$. For the case of hopping Hamiltonians of
the form (15) this implies only short-distance hopping, with both the hopping
constants and the potential terms depending on site number $n$ only via a
combination $an$. We then search the wave function in the form of the
following asymptotic series:

$\psi_{n}=\exp\left\\{\frac{iS\left(x\right)}{a}\right\\}\left[\varphi_{0}\left(x\right)+a\,\varphi_{1}\left(x\right)+...\right],\,\,\,\,x=an$
(72)

where $S\left(x\right)$ and $\varphi_{i}\left(x\right)$ are slow functions of
their arguments. The action of the Hamiltonian on such wave functions can be
expanded in powers of small $a$:

$\displaystyle\sum_{m}\mathcal{H}_{nm}\psi_{m}$
$\displaystyle=\sum_{m}\mathcal{H}_{nm}\exp\left\\{\frac{iS\left(am\right)}{a}\right\\}\left[\varphi_{0}\left(am\right)+a\,\varphi_{1}\left(am\right)+...\right]$
$\displaystyle\approx\exp\left\\{\frac{iS\left(x\right)}{a}\right\\}\mathcal{H}\left(x,p\right)\left[\varphi_{0}\left(x\right)+a\,\varphi_{1}\left(x\right)+...\right]$
$\displaystyle+\left(-ia\right)\,\exp\left\\{\frac{iS\left(x\right)}{a}\right\\}\,\left\\{\partial_{p}\mathcal{H}\left(x,p\right)\,\varphi_{0}^{\prime}\left(x\right)+\frac{1}{2}\,S^{\prime\prime}\left(x\right)\,\partial_{p}^{2}\mathcal{H}\left(x,p\right)\,\varphi_{0}\left(x\right)\right\\}+O\left(a^{2}\right),$
(73)

where $x=an$, $p=dS/dx$ is the semiclassical momentum, and
$\mathcal{H}\left(n,p\right)$ is the Wigner’s transform of the Hamiltonian:

$\mathcal{H}\left(x=an,p\right)=\sum_{m}\mathcal{H}_{nm}\exp\left\\{ip\left(n-m\right)\right\\}.$
(74)

According to the locality and smoothness conditions, this function depends
smoothly on both of its arguments, with $n$-dependence expressed via the
combination $an$, $a\ll 1$. The two leading orders in powers of $a$ of Eq.
(73) and the Schrodinger equation (71) imply the classical equation of motion

$\mathcal{H}\left(x,p\right)=E,$ (75)

as well as an equation on $\varphi_{0}$ called the transport equations (as it
ensures the conservation of probability density):

$\left\\{\partial_{p}\mathcal{H}\left(x,p\right)\,\frac{d}{dx}+\frac{1}{2}\,S^{\prime\prime}\left(x\right)\,\partial_{p}^{2}\mathcal{H}\left(x,p\right)\right\\}\,\varphi_{0}\left(x\right)=0,$
(76)

with similar equations also available for higher orders
$\varphi_{i},\,\,\,i\geq 1$, if one were to retain higher orders in the
expansion (LABEL:action-of-operator-on-semiclassical-wave-function). The
transport equation can be solved explicitly, yielding

$\varphi_{0}\left(x\right)=\frac{\text{const}}{\sqrt{\partial_{p}\mathcal{H}\left(x,p\right)}}\exp\left\\{\frac{1}{2}\intop^{x}dx\,\frac{\partial_{x}\partial_{p}\mathcal{H}\left(x,p\right)}{\partial_{p}\mathcal{H}\left(x,p\right)}\right\\}$
(77)

where the position-dependent value of the classical wave length
$p\left(x\right)$ found from Eq. (75) is implied. In particular, if the
Hamiltonian has the form $\mathcal{H}_{\text{kin}}+\mathcal{H}_{\text{pot}}$,
with kinetic term being position-independent and $H_{\text{pot}}$ being
diagonal, the cross term $\partial_{x}\partial_{p}\mathcal{H}$ vanishes, and
one arrives to the standard form of the WKB anstaz for the wave function:

$\psi_{n}=\frac{\text{const}}{\sqrt{\left|v\left(x\right)\right|}}\exp\left\\{\frac{i}{a}\intop^{x}p\left(x\right)\>dx\right\\},\,\,x=an$
(78)

where $v\left(x\right)=\partial_{p}\mathcal{H}\left(x,p\left(x\right)\right)$
is the classical velocity, and $p\left(x\right)$ is the classical wave length
found from the classical equation of motion, Eq. (75), with the latter having
solutions with two different signs of $p$, corresponding to particles
propagating in different directions according to the sign of the semiclassical
velocity.

The ansatz (72) is an asymptotic series whose applicability is controlled by
the two initial assumptions of locality and smoothness of the Hamiltonian. _A
posteriori_ , they are expressed in the requirement for the de Broglie wave
length $\lambda\left(x\right)=2\pi a/p\left(x\right)$ to change slowly with
$x=an$ at the scale of the wave length itself:

$\frac{d\lambda}{dx}\lambda\ll\lambda\Leftrightarrow
a\,\frac{1}{p^{2}}\frac{\partial_{x}\mathcal{H}}{\partial_{p}\mathcal{H}}\ll
1,$ (79)

While this condition is satisfied in general case due to smallness of $a$, it
is clearly violated if either $p$ or $v=\partial_{p}\mathcal{H}$ approach
zero, which correspond to classical turning points. In this region, the
semiclassical expansion (72) is not applicable, and one has to employ the
matching procedure.

###  A.2 Matching conditions and quantization rules for the cosine spectrum

For the particular case of Hamiltonian (15), the classical Hamiltonian is
given by $-J\cos p+h\left(x\right)$, and the positions of the classical
turning points are described by $\sin p=0$, or $E-h\left(x\right)=\pm J$. If
$dp/dx\neq 0$ at such point (where $p(x)$ is found from the classical equation
of motion), the $x$ axis is split by the turning point into classically
allowed and classically forbidden region. At sufficiently large distances from
the tuning point, the semiclassical solution becomes applicable, but one still
has to connect the amplitudes of the semiclassical wave functions at either
sides of the turning point. To do that, one needs to provide a solution for
the vicinity of the tuning point, possibly exploiting. These two solutions are
then matched together by their asymptotic forms. In order for this procedure
to work, there has to be an overlap in the regions of applicability of the two
approximations with spatial sizes that fit many de Broglie wave length, so the
phase change
$\frac{1}{a}\intop^{x}p\left(x\right)dx\propto\intop^{x}dx/\lambda(x)$ is
large as one traverses the overlap region.

One particularly important case where the vicinity of the critical point can
be analyzed is the case of smooth potential that can be expanded as
$h\left(x\right)\approx h_{0}+\alpha\left(x-x_{0}\right)$ in a sufficiently
large neighborhood of the classical turning point $x_{0}$. One then
distinguishes two cases: standard turning point $p=0$, corresponding to
$h_{0}=E+J$, and anomalous tuning point, for which one has $h_{0}=E-J$. Both
these cases are analyzed below.

####  A.2.1 Standard turning point

Let $\alpha>0$ be the potential gradient, and $E$ is such that
$\alpha\,a\delta=E+J$ for some $\delta\in\left[0,1\right]$ (one can always
shift the numeration of sites in order to force the exact position of the
tuning point to lie in the $[0,1]$ section). The exact Shroedinger equation
(71) then reads

$-\frac{J}{2}\left[\psi_{n+1}+\psi_{n-1}\right]+\alpha\,an\,\psi_{n}=\left(\alpha\,a\delta-J\right)\psi_{n}.$
(80)

Although this equation is exactly solvable in terms of Bessel functions, we
will resort to an controllable approximation that exploits the smoothness of
the potential and formally corresponds to asymptotic expansion of the exact
solution in powers of the gradient of the potential $\alpha a\ll 1$. Near
$n=0$, we expect to find a solution that is a smooth function of $x=an$, as
the relevant classical momenta are small due to the proximity to the standard
turning point where $p\left(x_{0}=a\delta\right)=0$. One can then replace the
discrete difference $\psi_{n+1}+\psi_{n-1}-2\psi_{n}$ by its continuous
approximation $a^{2}d^{2}\psi/dx^{2}+O\left(a^{4}\right)$, which is equivalent
to replacing the true $p$-dependence of the classical Hamiltonian (19) with
its quadratic approximation near $p=0$. The criterion of applicability for
such a replacement is given by the smallness of the momentum that validates
the quadratic expansion:

$\left|p\right|^{2}\approx\left|\frac{E+J-\alpha a\,x}{J}\right|\ll
1\Leftrightarrow\frac{\alpha}{J}\left|x-a\delta\right|\ll 1,$ (81)

where the momentum is estimated from the energy conservation. Performing the
substitution, one obtains the following differential equation:

$-\frac{Ja^{2}}{2}\,\psi^{\prime\prime}\left(x\right)+\alpha
x\,\psi\left(x\right)=\alpha\,a\delta\,\psi\left(x\right),$ (82)

which is exactly solvable in terms of the Airy functions:

$\psi\left(x\right)=C_{1}\text{Ai}\left(K\left(x-a\delta\right)\right)+C_{2}\text{Bi}\left(K\left(x-a\delta\right)\right),\,\,\,K=\left(\frac{2\alpha}{Ja^{2}}\right)^{1/3}\gg
1.$ (83)

Despite the condition (81), $K\left(x-a\delta\right)$ can attain large values,
so it is legitimate to use the asymptotic expansion of the Airy functions in
the region $K^{-1}\ll\left|x-a\delta\right|\ll\left(\alpha/J\right)^{-1}$,
rendering

$\psi_{n}\approx\begin{cases}z>0:&C_{1}\,\frac{1}{2\sqrt{\pi}}\sqrt[4]{\frac{1}{z}}\exp\left(-\frac{2}{3}z^{3/2}\right)+C_{2}\,\frac{1}{\sqrt{\pi}}\sqrt[4]{\frac{1}{z}}\exp\left(+\frac{2}{3}z^{3/2}\right),\\\
z<0:&C_{1}\,\frac{1}{\sqrt{\pi}}\sqrt[4]{\frac{1}{\left|z\right|}}\cos\left(-\frac{2}{3}\left|z\right|^{3/2}+\frac{\pi}{4}\right)+C_{2}\,\frac{1}{\sqrt{\pi}}\sqrt[4]{\frac{1}{\left|z\right|}}\cos\left(\frac{2}{3}\left|z\right|^{3/2}+\frac{\pi}{4}\right),\end{cases}\,\,z=K\left(x-a\delta\right).$
(84)

On the other hand, the semiclassical approximation can be applied in the
region described by its applicability condition (79), which amounts to
$K\left|x-a\delta\right|\gg 1$. The semiclassical wave functions then reads:

$\psi\left(x\right)\approx\begin{cases}\text{allowed
region}:&\frac{1}{\sqrt{\left|v\left(x\right)\right|}}\exp\left\\{\pm\frac{i}{a}\intop_{a\delta}^{x}p\left(x\right)dx\right\\},\\\
\text{forbiden
region}:&\frac{1}{\sqrt{\left|v\left(x\right)\right|}}\exp\left\\{\pm\frac{1}{a}\intop_{a\delta}^{x}\left|p\left(x\right)\right|dx\right\\},\end{cases}$
(85)

where $p\left(x\right)=\arccos\left(h\left(x\right)-E\right)/J$, and
$v\left(x\right)=J\sin p\left(x\right)$. Given the condition (81), one can
replace $\arccos 1-u\approx\sqrt{2u}$, which is equivalent to replacing the
kinetic energy with its quadratic approximation, so the semiclassical
approximation for the wave function is given by

$\psi_{n}\approx\begin{cases}z>0:&C_{\text{decaying}}\frac{1}{z^{1/4}}\exp\left\\{-\frac{2}{3}z^{3/2}\right\\}+C_{\text{growing}}\frac{1}{z^{1/4}}\exp\left\\{+\frac{2}{3}z^{3/2}\right\\},\\\
z<0:&C_{\text{towards}}\frac{1}{\left|z\right|^{1/4}}\exp\left\\{+\frac{2i}{3}\left|z\right|^{3/2}\right\\}+C_{\text{away}}\frac{1}{\left|z\right|^{1/4}}\exp\left\\{-\frac{2i}{3}\left|z\right|^{3/2}\right\\},\end{cases}\,\,\,\,z=K\left(n-a\delta\right),$
(86)

where the labeling of the $C$ coefficients represents the direction of
propagation w.r.t the turning point according to the sign of the group
velocity $v=J\sin p\approx Jp$ for $z>0$ and the behavior with the distance
from the tuning point for $z<0$. Both expressions (84) and (86) are applicable
in the same region, which contains $\frac{2}{3}\left(Ka\right)^{-1}$ de
Broglie wave lengths, so for the matching procedure to be valid it is
necessary to have $Ka\ll 1\Leftrightarrow\alpha a/J\ll 1$, which is yet
another expression of the smoothness of the potential. Matching the two
asymptotic expressions then yields

$\begin{cases}n>\delta:&C_{\text{decaying}}=C_{1}\,\frac{1}{2\sqrt{\pi}},\,\,\,C_{\text{growing}}=C_{2}\,\frac{1}{\sqrt{\pi}},\\\
n<\delta:&C_{\text{towards}}=C_{1}\,\frac{1}{\sqrt{\pi}}\frac{e^{i\pi/4}}{2}+C_{2}\,\frac{1}{\sqrt{\pi}}\frac{e^{-i\pi/4}}{2},\,\,\,C_{\text{away}}=C_{1}\,\frac{1}{\sqrt{\pi}}\frac{e^{-i\pi/4}}{2}+C_{2}\,\frac{1}{\sqrt{\pi}}\frac{e^{i\pi/4}}{2}.\end{cases}$
(87)

Excluding $C_{1},C_{2}$, one obtains the well-known [18] connection between
the coefficients of the semiclassical expansion on either side of the standard
turning point:

$\begin{cases}C_{\text{towards}}=C_{\text{decaying}}e^{i\pi/4}+C_{\text{growing}}\frac{e^{-i\pi/4}}{2},\\\
C_{\text{away}}=C_{\text{decaying}}e^{-i\pi/4}+C_{\text{growing}}\frac{e^{i\pi/4}}{2}.\end{cases}$
(88)

Although these relations are obtain for the case $\alpha>0$, they remain
identical for the case $\alpha<0$ upon correctly identifying growing and
decaying solutions in the forbidden region and left- and right-propagating
waves in the allowed region.

####  A.2.2 Anomalous turning point

Let $\alpha>0$ be the potential gradient, and $E$ is such that
$\alpha\,a\delta=E-J$ for some $\delta\in\left[0,1\right]$. This corresponds
to anomalous turning described by $p\left(x_{0}=a\delta\right)=\pm\pi$, with
$x>a\delta$ being the classically allowed region and $x<a\delta$ being the
classically forbidden one. To solve the exact Schroedinger equation (71) for
this case, we first perform the substitution
$\psi_{n}=\eta_{n}\left(-1\right)^{n}$, which maps the problem to the previous
case, Eq. (80), with reversed sign of $\alpha$:

$-\frac{J}{2}\left[\psi_{n+1}+\psi_{n-1}\right]+\alpha\,an\,\psi_{n}=\left(\alpha\,a\delta+J\right)\psi_{n}\Leftrightarrow-\frac{J}{2}\left[\eta_{n+1}+\eta_{n-1}\right]-\alpha\,an\,\eta_{n}=\left(-\alpha\,a\delta-J\right)\eta_{n}.$
(89)

Similarly to the previous subsection, the latter equation is solved by Airy
function, so in the region
$K^{-1}\ll\left|x-a\delta\right|\ll\left(\alpha/J\right)^{-1}$ the asymptotic
form of the solution is given by:

$\eta_{n}=\psi_{n}\left(-1\right)^{n}=\begin{cases}z<0:&C_{1}\,\frac{1}{2\sqrt{\pi}}\sqrt[4]{\frac{1}{\left|z\right|}}\exp\left(-\frac{2}{3}\left|z\right|^{3/2}\right)+C_{2}\,\frac{1}{\sqrt{\pi}}\sqrt[4]{\frac{1}{\left|z\right|}}\exp\left(\frac{2}{3}\left|z\right|^{3/2}\right),\\\
z>0:&C_{1}\,\frac{1}{\sqrt{\pi}}\sqrt[4]{\frac{1}{z}}\cos\left(-\frac{2}{3}z^{3/2}+\frac{\pi}{4}\right)+C_{2}\,\frac{1}{\sqrt{\pi}}\sqrt[4]{\frac{1}{z}}\cos\left(\frac{2}{3}z^{3/2}+\frac{\pi}{4}\right),\end{cases}\,\,\,z=-K\left(x-a\delta\right),$
(90)

where $K=\left(2\alpha/Ja^{2}\right)^{1/3}$. On the other hand, the
semiclassical momentum is close to $\pm\pi$, so semiclassical wave function is
described by

$\psi\left(x\right)\approx\begin{cases}\text{allowed
region}:&\frac{1}{\sqrt{\left|v\left(x\right)\right|}}\exp\left\\{\pm\frac{i}{a}\intop_{a\delta}^{x}\left[\pi-\delta
p\left(x\right)\right]dx\right\\},\\\ \text{forbiden
region}:&\frac{1}{\sqrt{\left|v\left(x\right)\right|}}\exp\left\\{\pm\frac{1}{a}\intop_{a\delta}^{x}\left[i\pi+\left|\delta
p\left(x\right)\right|\right]dx\right\\},\end{cases}$ (91)

where $\delta
p=\pi-\arccos\left\\{\left(h\left(x\right)-E+J\right)/J-1\right\\}\approx\sqrt{2\left(E-J-h\left(x\right)\right)/J}$
(positive quantity in the allowed region), while the group velocity reads
$v=J\sin\left(\pi-\delta p\right)\approx J\delta p$. The anomalous turning
point is characterized by the fact that the real contribution to the phase
$\pm i\pi\left(x-a\delta\right)/a=\pm i\pi\left(n-\delta\right)$ is present at
both sides of the turning point and causes oscillatory behavior of the wave
function. The semiclassical solution then reads

$\psi_{n}\approx\begin{cases}z<0:&C_{\text{decaying}}\frac{1}{\left|z\right|^{1/4}}\exp\left\\{+i\pi\frac{x-a\delta}{a}-\frac{2}{3}\left|z\right|^{3/2}\right\\}+C_{\text{growing}}\frac{1}{\left|z\right|^{1/4}}\exp\left\\{-i\pi\frac{x-a\delta}{a}+\frac{2}{3}\left|z\right|^{3/2}\right\\},\\\
z>0:&C_{\text{towards}}\frac{1}{\left|z\right|^{1/4}}\exp\left\\{-i\pi\frac{x-a\delta}{a}+\frac{2i}{3}z^{3/2}\right\\}+C_{\text{away}}\frac{1}{\left|z\right|^{1/4}}\exp\left\\{+i\pi\frac{x-a\delta}{a}-\frac{2i}{3}z^{3/2}\right\\},\end{cases}\,\,\,z=-K\left(x-\delta\right),$
(92)

where the labeling of the constants is again chosen to reflect the
characteristic behavior of the wave functions (e.g., the group velocity is
positive for right-propagating wave, which corresponds to moving away from the
turning point). Performing the matching procedure in the region
$K^{-1}\ll\left|x-a\delta\right|\ll\left(\alpha/J\right)^{-1}$ then yields

$\begin{cases}C_{\text{towards}}=C_{\text{decaying}}\exp\left\\{-2i\pi\delta\right\\}e^{-i\frac{\pi}{4}}+C_{\text{growing}}\frac{1}{2}e^{i\frac{\pi}{4}},\\\
C_{\text{away}}=C_{\text{decaying}}e^{i\frac{\pi}{4}}+C_{\text{growing}}\exp\left\\{+2i\pi\delta\right\\}\frac{1}{2}e^{-i\frac{\pi}{4}}.\end{cases}$
(93)

where we have taken into account that $\exp\left\\{\pm i\pi
x/a\right\\}=\exp\left\\{\pm i\pi n\right\\}=\left(-1\right)^{n}.$

One important feature of this matching condition in comparison to the standard
one, Eq. (88), is that it includes and additional phase factor featuring the
exact position $\delta\in\left[0,1\right]$ of the turning point, which has
been a point of controversy for many authors [19]. While one might argue that
the potential is not even well-defined at non-integer coordinate values of
$n=x/a$, it becomes well-defined after one requires it to be a smooth function
of the coordinate. Due to the fact that the wave function has momentum
$p=\pm\pi$ at the turning point, this exact coordinate becomes essential. It
is important for accurate description of the quantization conditions.

####  A.2.3 Quantization conditions

For a given energy $E$, consider a potential that arranges a bounded
classically allowed region $x_{1}<x<x_{2}$, where $x_{1,2}$ are turning
points. If one does not consider tunneling outside this region, the matching
conditions allow one to derive the quantization rules for the bound states.
The latter is defined as a purely decaying wave in the forbidden region close
to both of the tuning points $x_{i},i=1,2$:

$C_{\text{grow}}^{\left(i\right)}=0\Rightarrow\frac{C_{\text{towards}}^{\left(i\right)}}{C_{\text{away}}^{\left(i\right)}}=\begin{cases}e^{i\pi/2},&p\left(x_{i}\right)=0,\\\
e^{-i\pi/2}\exp\left\\{-2i\pi\delta^{\left(i\right)}\right\\},&p\left(x_{i}\right)=\pm\pi.\end{cases}$
(94)

where $\delta^{\left(i\right)}\in\left[0,1\right]$ is the exact position of
the anomalous turning point. Obviously, one has
$\left|C_{\text{towards}}^{\left(i\right)}/C_{\text{away}}^{\left(i\right)}\right|=1$,
so the phase
$\theta_{i}=\text{arg}C_{\text{towards}}^{\left(i\right)}/C_{\text{away}}^{\left(i\right)}$
defines the shape of the semiclassical wave function in the allowed region up
to an overall phase as seen defined by the corresponding turning point:

$\psi^{\left(1\right)}\left(x\right)=\frac{C_{\text{away}}^{\left(1\right)}}{\sqrt{\left|v\left(x\right)\right|}}\exp\left\\{\frac{i}{a}\intop_{x_{1}}^{x}p\left(x\right)dx\right\\}+\frac{C_{\text{towards}}^{\left(1\right)}}{\sqrt{\left|v\left(x\right)\right|}}\exp\left\\{-\frac{i}{a}\intop_{x_{1}}^{x}p\left(x\right)dx\right\\}=\frac{2C_{\text{away}}^{\left(1\right)}e^{\theta_{1}/2}}{\sqrt{\left|v\left(x\right)\right|}}\cos\left\\{\frac{1}{a}\intop_{x_{1}}^{x}p\left(x\right)dx-\frac{\theta_{1}}{2}\right\\},$
(95)
$\psi^{\left(2\right)}\left(x\right)=\frac{C_{\text{towards}}^{\left(2\right)}}{\sqrt{\left|v\left(x\right)\right|}}\exp\left\\{\frac{i}{a}\intop_{x_{2}}^{x}p\left(x\right)dx\right\\}+\frac{C_{\text{away}}^{\left(2\right)}}{\sqrt{\left|v\left(x\right)\right|}}\exp\left\\{-\frac{i}{a}\intop_{x_{2}}^{x}p\left(x\right)dx\right\\}=\frac{2C_{\text{away}}^{\left(2\right)}e^{i\theta_{2}/2}}{\sqrt{\left|v\left(x\right)\right|}}\cos\left\\{\frac{1}{a}\intop_{x_{2}}^{x}p\left(x\right)dx+\frac{\theta_{2}}{2}\right\\}$
(96)

where we have taken into account $x_{1}<x<x_{2}$ , and $p$ is chosen to be a
positive root of the classical equation of motion (75) in the vicinity of the
tuning point $x_{i}$ (so the corresponding group velocity is also positive).
The two expressions for $i=1,2$ should actually yield the same wave function
(up to an overall sign), which is equivalent to requiring that the difference
of phases of the cosines is $\pi N$, yielding

$\frac{1}{a}\intop_{x_{1}}^{x_{2}}p\left(x\right)dx-\frac{\theta_{1}}{2}-\frac{\theta_{2}}{2}=\pi
N,$ (97)

which is equivalent to Eqs. (22) and (23) upon substitution of
$\theta_{i}=\text{arg}C_{\text{towards}}^{\left(i\right)}/C_{\text{away}}^{\left(i\right)}$
from Eq. (94).

Let us illustrate this result for an exactly solvable case of linear potential
$h_{n}=\alpha an$ with $a\ll 1$. The exact eigenvalues of the corresponding
Schroedinger equation (71) are given by $E_{N}=\frac{J}{\alpha}+\alpha
aN,\,\,N\in\mathbb{Z}$, and the associated eigenfunctions read

$\left\langle n|N\right\rangle=J_{n-N}\left(\frac{J}{\alpha}\right),$ (98)

where $J_{k}\left(x\right)$ is the Bessel function of order $k$ (which
promptly reproduces the special cases of Eqs. (84) and (90) with $C_{2}=0$
near the corresponding tuning points for $a\ll 1$). From the semiclassical
point of view, at energy $E$ the classically allowed region of Hamiltonian
$-J\cos p+\alpha x$ is given by $x_{1}<x<x_{2}$ with $x_{1,2}=\left(E\mp
J\right)/\alpha$, and $x_{1}$ is the anomalous turning point with
$p\left(x_{1}\right)=\pm\pi$, while $x_{2}$ is the standard turning point with
$p\left(x_{2}\right)=0$. Applying Eqs. (94) and (97) then renders
$\theta_{2}=+\pi/2$, and $\theta_{1}=-\pi/2-2\pi\left\\{x_{1}/a\right\\}$,
with $\left\\{n\right\\}$ corresponding to the fractional part of $n$, so

$\frac{1}{a}\intop_{x_{1}}^{x_{2}}\arccos\left\\{\frac{\alpha
x-E}{J}\right\\}dx+\pi\left\\{x_{1}/a\right\\}=\pi N,$ (99)

The integral of momentum is $\pi J/\alpha$, so one obtains an equation on the
value of $E$:

$\frac{J}{\alpha a}+\left\\{\frac{E}{\alpha a}-\frac{J}{\alpha
a}\right\\}\in\mathbb{Z},$ (100)

which is solved by the correct value $E=\alpha aN$. Crucially, the correct
quantization condition can be traced down to the presence of the additional
dependence of the matching conditions (94) on $\delta$ .

##  Appendix B Correction to the semiclassical description due to
Trotterization

This Appendix covers several technical details of calculations leading to the
results of (III) by using the framework of semiclassical approximation for the
matrix Hamiltonians as presented in Ref. [35, Ch. 11]., with the latter
naturally arising due to the different periodicity of the target Hamiltonian
(15) and the Suzuki-Trotter approximation (7).

###  B.1 Expression for the 2nd-order correction

The matrix elements of the difference between the single-particle effective
Hamiltonian (8) of the Suzuki-Trotter evolution and the target Hamiltonian
(15) are presented on B1. Clearly, they are dominated by the close-neighbor
hopping, which corresponds to leading-order perturbation theory in powers of
the Trotter step $\delta t$. The Suzuki-Trotter evolution operator composed of
$k$ terms is described by the following recursion relation:

$\mathcal{U}_{k}\left(\delta t\right)=e^{-iA_{k}\delta
t/2}\mathcal{U}_{k-1}\left(\delta t\right)e^{-iA_{k}\delta
t/2},\,\,\,\mathcal{U}_{1}=e^{-iA_{1}t}.$ (101)

The effective Hamiltonian at each step has the following expansion in powers
of $t$:

$\mathcal{H}_{\text{eff}}^{\left(k\right)}=-\frac{1}{i\delta
t}\text{Ln}\mathcal{U}_{k}\left(\delta t\right)=\mathcal{H}_{\text{
}}^{\left(k\right)}+\delta t^{2}\mathcal{D}^{\left(k\right)}+O\left(\delta
t^{4}\right),\,\,\,\mathcal{H}_{\text{
}}^{\left(k\right)}=\sum_{j=1}^{k}A_{j},$ (102)

and we are interested in the value of $\mathcal{D}^{\left(k\right)}$. To do
this, we apply the second order of Baker–Campbell–Hausdorff formula:

$\text{Ln}\left\\{e^{-itA}e^{-itB}\right\\}=-it\left[A+B+\frac{t}{2i}\left[A,B\right]+\frac{t^{2}}{12}\left(\left[A,\left[A,B\right]\right]+\left[B,\left[B,A\right]\right]\right)+O\left(t^{3}\right)\right].$
(103)

Used twice on Eq. (101) and then substituted in Eq. (102), this yields Eq.
(31) of the main text:

$\mathcal{D}^{\left(k+1\right)}=\mathcal{D}^{\left(k\right)}-\frac{1}{12}\left\\{\left[\mathcal{H}_{\text{
}}^{\left(k\right)},\left[\mathcal{H}_{\text{
}}^{\left(k\right)},A_{k+1}\right]\right]-\frac{1}{2}\left[A_{k+1},\left[A_{k+1},\mathcal{H}_{\text{
}}^{\left(k\right)}\right]\right]\right\\}.$ (104)

If the target Hamiltonian can be decomposed as $\mathcal{H}=A_{1}+A_{2}$,
where both $A_{1,2}$ can be exponenitated exactly, the leading-order
correction to the effective Hamiltonian simply reads

$\mathcal{H}_{\text{eff}}=\mathcal{H}_{\text{ }}+\delta
t^{2}\mathcal{D}^{\left(2\right)}+O\left(\delta
t^{4}\right),\,\,\,\mathcal{D}^{\left(2\right)}=-\frac{1}{12}\left\\{\left[A_{1},\left[A_{1},A_{2}\right]\right]-\frac{1}{2}\left[A_{2},\left[A_{2},A_{1}\right]\right]\right\\}.$
(105)

Note that this expression is not symmetric to the order of $A_{i}$. Similarly,
for the three-term decomposition $\mathcal{H}=A_{1}+A_{2}+A_{3}$ one has

$\mathcal{H}_{\text{eff}}=\mathcal{H}_{\text{ }}+\delta
t^{2}\mathcal{D}^{\left(3\right)}+O\left(\delta t^{4}\right),$ (106)
$\mathcal{D}^{\left(3\right)}=-\frac{1}{12}\left\\{\left[A_{1},\left[A_{1},A_{2}\right]\right]-\frac{1}{2}\left[A_{2},\left[A_{2},A_{1}\right]\right]+\left[A_{1}+A_{2},\left[A_{1}+A_{2},A_{3}\right]\right]-\frac{1}{2}\left[A_{3},\left[A_{3},A_{1}+A_{2}\right]\right]\right\\}.$
(107)

Clearly, the expression is not invariant w.r.t permutation of indices. We will
later show, however, that the difference is inessential as far as the
semiclassical approximation is concerned.

Figure B1: Visualization of the spatial matrix elements of the Suzuki-Trotter
error $\left\langle n\left|\delta\mathcal{H}\right|m\right\rangle$ normalized
by $J\,\left(J\delta t\right)^{2}/24$ with the same parameters as 4 but
$J\delta t=0.3$. _Left:_ color plot of the matrix elements _. Right:_ matrix
element as a function of distance between sites $\left|m-n\right|$ for fixed
$n=L/2$, with the blue line corresponding to $\left(J\delta
t/4\right)^{\left|m-n\right|-3}$, according to estimation (29). On both plots,
red dashed line denotes $\left|m-n\right|=4$.

###  B.2 Treatment of alternating perturbation

As we can see from, e.g., B1, even in the case of a smooth potential, the
effective Hamiltonian $\mathcal{H}_{\text{eff}}$ features a “staggered”
structure:

$\delta\mathcal{H}^{ij}:=\mathcal{H}_{\text{eff}}^{ij}-\mathcal{H}^{ij}=U^{ij}+\left(-1\right)^{i}W^{ij}+\left(-1\right)^{j}\overline{W^{ij}},$
(108)

where both $U$ and $W$ depend smoothly on $m+n$ and decay quickly with $m-n$.
In this case, the spatial dependence of the wave function is not smooth
anymore. It essentially implies that even in the limit of homogeneous system
the wave function has two components on even and odd sites respectively, thus
reflecting the doubling of the translation period of the system. The
corresponding Brillouin zone is halved to $p\in\left[-\pi/2,\pi/2\right]$.

Let’s apply the operator of the form (108) to a semiclassical wave function. A
calculation identical to that of A.1 yields:

$\displaystyle\sum_{m}\mathcal{H}_{\text{eff}}^{nm}\psi_{m}$
$\displaystyle\approx\exp\left\\{\frac{iS\left(x\right)}{a}\right\\}\left[\mathcal{H}\left(x,p\right)\varphi-\left(ia\right)\left\\{\partial_{p}\mathcal{H}\left(x,p\right)\,\varphi^{\prime}+\frac{1}{2}\,p^{\prime}\,\partial_{p}^{2}\mathcal{H}\left(x,p\right)\,\varphi\right\\}+O\left(a^{2}\right)\right]$
$\displaystyle\exp\left\\{\frac{iS\left(x\right)}{a}\right\\}\left[U\left(x,p\right)\varphi-\left(ia\right)\left\\{\partial_{p}U\left(x,p\right)\,\varphi^{\prime}+\frac{1}{2}\,p^{\prime}\,\partial_{p}^{2}U\left(x,p\right)\,\varphi\right\\}+O\left(a^{2}\right)\right]$
$\displaystyle+\left(-1\right)^{n}\exp\left\\{\frac{iS\left(x\right)}{a}\right\\}\left[W\left(x,p\right)\varphi-\left(ia\right)\left\\{\partial_{p}W\left(x,p\right)\,\varphi^{\prime}+\frac{1}{2}\,p^{\prime}\,\partial_{p}^{2}W\left(x,p\right)\,\varphi\right\\}+O\left(a^{2}\right)\right],$
(109)

where $p=S^{\prime}\left(x\right)$, $\mathcal{H}\left(x,p\right)$ is given by
Eq. (112), and

$U\left(x=an,p\right)=\sum_{m}U^{nm}e^{ip\left(m-n\right)},\,\,\,W\left(x=an,p\right)=\sum_{m}\left[W^{nm}e^{ip\left(m-n\right)}+\overline{W^{mn}}e^{i\left(p+\pi\right)\left(m-n\right)}\right].$
(110)

We thus need to consider the two-component wave function:

$\psi_{n}=\exp\left\\{\frac{iS\left(an\right)}{a}\right\\}\left(\varphi\left(an\right)+\left(-1\right)^{n}\eta\left(an\right)\right).$
(111)

Application of the full effective Hamiltonian $\mathcal{H}_{\text{eff}}$ to
the such wave function then yields

$\sum_{m}\mathcal{H}_{\text{eff}}^{nm}\psi_{m}\approx\exp\left\\{\frac{iS\left(an\right)}{a}\right\\}\,\left[\hat{\mathcal{H}}_{\text{eff}}\left(x,p\right)-\left(ia\right)\left\\{\partial_{p}\hat{\mathcal{H}}_{\text{eff}}\left(x,p\right)\,\frac{d}{dx}+\frac{1}{2}\,p^{\prime}\,\partial_{p}^{2}\hat{\mathcal{H}}_{\text{eff}}\left(x,p\right)+i\frac{\partial\,\hat{\mathcal{H}}_{\text{eff}}\left(x,p\right)}{\partial
a}\right\\}+O\left(a^{2}\right)\right]\overrightarrow{v},$ (112)

where $p\left(x\right)=S^{\prime}\left(x\right)$, and the matrix
$\hat{\mathcal{H}}_{\text{eff}}$ and vector $\overrightarrow{v}$ are given by

$\hat{\mathcal{H}}_{\text{eff}}\left(x,p\right)=\begin{pmatrix}\mathcal{H}\left(x,p\right)+U\left(x,p\right)&W\left(x,p+\pi\right)\\\
W\left(x,p\right)&\mathcal{H}\left(x,p+\pi\right)+U\left(x,p+\pi\right)\end{pmatrix},\,\,\,\overrightarrow{v}\left(x\right)=\begin{pmatrix}\varphi\left(x\right)\\\
\eta\left(x\right)\end{pmatrix}.$ (113)

In Eq. (112), we have also included the explicit dependence of
$\hat{\mathcal{H}}$ on $a$, which will come handy in the following sections.

Following Ref. [35], we then seek the solution to the Schroedinger equation
expansion in powers of $a$: