# Quantum Simulation of Lindbladian Dynamics via Repeated Interactions

Matthew Pocrnic  Department of Physics, 60 Saint George St., University of
Toronto, Toronto, Ontario, M5S 1A7, Canada  Dvira Segal  Department of
Chemistry and Centre for Quantum Information and Quantum Control, University
of Toronto, 80 Saint George St., Toronto, Ontario, M5S 3H6, Canada  Department
of Physics, 60 Saint George St., University of Toronto, Toronto, Ontario, M5S
1A7, Canada  Nathan Wiebe  Department of Computer Science, University of
Toronto, Toronto ON, Canada  Pacific Northwest National Laboratory, Richland
Wa, USA  Canadian Institute for Advanced Research, Toronto ON, Canada

###### Abstract

The Lindblad equation generalizes the Schrödinger equation to quantum systems
that undergo dissipative dynamics. The quantum simulation of Lindbladian
dynamics is therefore non-unitary, preventing a naive application of state-of-
the-art quantum algorithms. Here, we make use of an approximate correspondence
between Lindbladian dynamics and evolution based on Repeated Interaction (RI)
CPTP maps to write down a Hamiltonian formulation of the Lindblad dynamics and
derive a rigorous error bound on the master equation. Specifically, we show
that the number of interactions needed to simulate the Liouvillian
$e^{t\mathcal{L}}$ within error $\epsilon$ scales in a weak coupling limit as
$\nu\in O(t^{2}\|\mathcal{L}\|_{1\rightarrow 1}^{2}/\epsilon)$. This is
significant because explicit error bounds in the Lindbladian approximation to
the dynamics are not explicitly bounded in existing quantum algorithms for
open system simulations.We then provide quantum algorithms to simulate these
maps using an iterative Qubitization approach and Trotter-Suzuki formulas and
specifically show that for iterative qubitization the number of operations
needed to simulate the dynamics (for a fixed value of $\nu$) scales in a weak
coupling limit as $O(\nu(t\alpha_{0}+\log(1/\epsilon)/\log\log(1/\epsilon)))$
where $\alpha_{0}$ is the coefficient $1$-norm for the system and bath
Hamiltonians. This scaling would appear to be optimal if the complexity of
$\nu$ is not considered, which underscores the importance of considering the
error in the Liouvillian that we reveal in this work.

##  I Introduction

The simulation of quantum systems is arguably the most natural application of
quantum computers [1]. The simulation problem can be summarized as follows:
For a given Hamiltonian $H$ and evolution time $t$, construct a unitary $U$
that approximates the dynamics generated by $H$ up to an error $\epsilon$,
with respect to some appropriate distance metric. In the case of time-
independent $H$ we wish to construct a quantum algorithm to build $U$, and
approximate the solution to the Schrödinger equation
$\partial_{t}\ket{\psi}=-iH\ket{\psi}$ such that
$\norm{U-e^{-iHt}}\leq\epsilon$. This problem has deservedly received much
attention, and the field has achieved many remarkable results [2, 3, 4, 5, 6].
However, significantly less attention has been given to non-unitary quantum
dynamics, or open quantum systems. These are systems that interact and
exchange energy with a larger environment. The field of open quantum systems
has put forward a variety of methods for studying such systems to varying
levels of approximation and accuracy [7]. In particular, the formalism of
quantum master equations (QMEs) has been proven useful to describe dissipative
quantum dynamics. A QME is differential equation for the system’s density
matrix dictating its evolving in time, incorporating the effects of both
coherent evolution and dissipative processes influenced by the environment
[8]. The approach we implement in this work is the Lindblad quantum master
equation [8, 9]. From the perspective of quantum computing, this is perhaps
the most natural master equation to implement given that it generates dynamics
described by completely positive trace preserving (CPTP) maps. This equation
also leads to Markovian dynamics, meaning that correlations between the bath
and the system do not build up. Physically, this means that the bath dynamics
occur at a much faster time scale than the system, so any correlations in the
bath due to the coupled system propagate to infinity before they influence the
system once again. Conceptually, one would therefore expect to be able to
simulate Markovian dynamics in a smaller Hilbert space, and therefore, with
less resources than non-Markovian quantum master equations.

In this work, we present quantum algorithms based on an existing approximate
correspondence between the Lindblad equation and repeated interaction (RI)
maps, also known as collision models. Repeated interactions model is a scheme
whereby the environment repeatedly interact with the system. The environment
comprises many independent subsystems with each interacting with the system
one by one. These bath’s subsystems do not interact with one another, and
never interact with the system again, naturally leading to a Markov process.
While this might seem contrived, open quantum systems such as the Maser can be
naturally described in this fashion: In this case, an electromagnetic field in
a cavity is repeatedly perturbed by particles passing through, whereby the
reduced dynamics of interest is that of the quantum field [10, 11]. Here, we
provide quantum algorithms to implement these repeated interactions to
effectively model Lindbladian dynamics. Our algorithms yield low order
integrators in general, which we expect to be nearly optimal given that the
error between the Lindblad dynamics and those generated by repeated
interactions is second order in the interaction time in general. Low order
formulas here are not as problematic as they may seem given the physics
captured by the Lindblad equation: This QME can be derived from microscopic
principles by integrating the exact equation of motion generated by an
interaction Hamiltonian $H_{I}$ that governs the entire system bath dynamics
[8]. This integral equation is then expanded to second order in the
interaction time, and a trace is performed over the environment Hilbert space.
This therefore gives an accurate description of the exact dynamics only up to
second order, which implies that higher order integrators are not necessarily
capturing additional physics.

##  II Contributions and Related Work

When it comes to related works, state of the art Lindbladian simulation is
presented in Refs. [12, 13]. The former utilizes Kraus maps to approximate the
Lindbladian to second order, and then constructs a variant of linear
combinations of unitaries (LCU) for channels to implement the Kraus maps in a
purified space. The approach then utilizes a compression encoding scheme
involving the Chernoff bound to achieve complexity
$O(t\text{polylog}(t/\epsilon))$. The latter approach [13] utilizes an
interesting application of Duhamel’s principle to build high-order integral
formulas, which are approximated as a sum, derived from Gaussian Quadrature
methods. Linear combinations of unitaries are then constructed once again for
implementation, again achieving time complexity
$O(t\text{polylog}(t/\epsilon))$, with slightly better dependence on the
number of dissipators in the Lindbladian. Other notable algorithmic
contributions in this field can be found in Refs. [14, 15, 16]. Alternative
implementations that also conduct hardware experiments can be found in Refs.
[17, 18]. Most recently [19], work has been done to formulate integrators of
arbitrary order based on repeated Hamiltonian simulations, which is
reminiscent of the theme of this work. The authors in Ref. [19] propose a
strategy that is also based on maps of a Stinespring form with $\sqrt{\tau}$
evolution of an effective Hamiltonian. However, this approach is derived from
a method known as unravelling, in which the Lindblad equation is transformed
into a stochastic Schrödinger equation. Then, the goal is to solve the
dynamics for the expectation value of the final state. This approach is
similar to an earlier correspondence proposed between the repeated
interactions scheme that we investigate here, and quantum Langevin equations
[20]. The intuition between this relation is that coupling of the system to a
random ancilla can lead to the generation of random quantum forces. Ref. [19]
also generalizes results to time dependent Lindbladians, a topic that has
seldom been studied in the context of quantum simulations.

In terms of our contributions, we were largely inspired by Ref. [21], which to
the best of our knowledge, was the first to propose quantum simulation of the
RI scheme. Our work builds upon yet differs significantly from this analysis.
Here, the implementation of the RI maps is abstracted from the physical
picture considered in Ref. [21], whereby we analyze the complexity of
converging to arbitrary time-independent Lindbladian dynamics using an RI map
with minimal assumptions. In doing so, we provide a rigorous error bound on
the difference between these two CPTP maps. This is one of the main
contributions of our work, which to the best of our knowledge has not been
done before in this context. Bounds given in Ref. [21] solve a problem with a
physical model in mind. There, a k-local Lindbladian is assumed to be
simulated with a physical multipartite interaction model in which subsystems
interact locally with a system governed by $H_{S}$, supported on some lattice.
This necessarily involves some complicated indices such as the set of all
ancilla-system site pairings, which can make the bounds difficult to
interpret. As well, for the purpose of simulation, the analysis in Ref. [21]
uses Trotter-Suzuki formulas to implement each propagator from the outskirts.
Here, we specifically bound the convergence of an RI map to Lindbladian
dynamics, where we use the number of interactions with ancillary subsystems to
control the error. We then use modern quantum algorithms based on block
encodings to implement these general maps as efficiently as possible. For
means of simplicity and comparison, we further provide a simulation of our RI
maps with Trotter-Suzuki formulas. We also note that multipartite and local
dynamics can be simulated with our model by restricting the support of the
operators in the unitary decomposition of the interaction Hamiltonian. We
summarize our contributions below:

  1. 1.

We rigorously bound the error between RI maps and Lindbladian dynamics in
terms of the number of environment subsystems $m$ and the number of
applications $\nu$ of an RI map composed of $m$ interactions in Theorem 2.

  2. 2.

In Corollary 2.1, we then immediately obtain a bound on the number of ancilla
bath’s subsystems, given that the Trace is performed over the entire
environment register at the end of the simulation.

  3. 3.

We describe a strategy to achieve the same simulation accuracy and time
complexity using a constant number of qubits, given that traces are performed
intermittently.

  4. 4.

For either of the ancilla strategy above, we implement a simulation of
Lindbladian dynamics using unitary RI maps, with a probability of success
$\delta=1$ and an error $\epsilon$, utilizing an iterative Qubitization
strategy to implement each discrete time map. We give a bound on the time
complexity in terms of queries to block encoding oracles in Theorem 3.

  5. 5.

To compare with Ref. [21], we implement our RI maps with Trotter-Suzuki
formulas, and bound the step size $r$ to achieve error $\epsilon$ for formulas
of first and arbitrary order, the former of which we utilize commutator
structure. This is presented in Theorems 5 and 6.

##  III Preliminaries

In this Section, we introduce the relevant notations and assumptions for the
work. We also provide a brief introduction to the Lindblad equation and
repeated interaction maps, and some of their relevant properties.

###  III.1 Notation and Assumptions

We use lowercase letters to denote variables, $a$, uppercase letters to
represent operators, $A$, and reserve calligraphic letters for super-
operators, $\mathcal{A}$. In this paper, we consider only finite-dimensional
Hilbert spaces, so all operators $A$ are matrices, and all super-operators can
be written as some functions of matrices. An important part of analyzing
quantum algorithms is bounding the error $\epsilon$, which allows one to
calculate the complexity of the algorithm. Since the algorithms here are
written in terms of super-operators, we implement an induced norm to calculate
$\epsilon$. The familiar Schatten p-norms for matrices are defined in the
following way: $\norm{A}_{p}=(\sum_{k}s_{k}^{p})^{\frac{1}{p}}$, where $s_{k}$
represent the singular values of the matrix $A$. The Schatten norms then
induce a norm on super-operators, that of focus here will be the induced
$1\rightarrow 1$ norm, which has the following definition:
$\norm{\mathcal{A}}_{1\rightarrow
1}=\max_{A:\norm{A}_{1}=1}\norm{\mathcal{A}[A]}_{1}$. This norm obeys the
triangle inequality and submultiplicativity, which we will make use of.
Throughout the work, we write $\rho$ to indicate a density matrix, which is a
positive semi-definite operator with $\Tr\rho=1$. Some other recurring
notations are the use of $\mathcal{L}$ as a Lindblad generator, $t$ for the
total simulation time, $m$ as the number of subsystems or ancilla used in the
repeated interactions scheme, $\nu$ as the number of applications of an RI map
with $m$ interactions, and $\tau=t/\nu$ as the interaction time interval. We
use $\circ$ to indicate the composition of super-operators, all of which are
quantum channels in this work. When the composition is indexed over many
channels, we use $\bigcirc$ to be distinct from a product.

###  III.2 Lindblad Equation

The Lindblad equation, often seen as a generalization of the Schrödinger
equation for Markovian open quantum systems, takes the following form
$(\hbar\equiv 1)$:

$\frac{d\rho_{S}}{dt}=-i[H_{0},\rho_{S}]+\sum_{j=1}^{l}L_{j}\rho_{S}L_{j}^{\dagger}-\frac{1}{2}\\{L_{j}^{\dagger}L_{j},\rho_{S}\\}\equiv\mathcal{L}(\rho_{S}),$
(1)

where the $L_{j}$ operators are derived from the interaction between the
system and environment. Therefore, this equation describes the dynamics of the
reduced state $\rho_{S}$ of a quantum system, and does not track the dynamics
of the environment. If $L_{j}=0~{}\forall~{}j$ then we recover the von Neumann
equation, or the Schödinger equation for density matrices. Here, the
Hamiltonian is a Hermitian operator, however, the jump operators $L_{j}$ are
not necessarily so. We also have $H,~{}L_{j}\in\mathbb{C}^{d\times d}$, where
$d$ is the dimension of the system.

###  III.3 Repeated Interactions Model

Repeated interactions allow us to study the reduced dynamics of a quantum
system under the influence of an environment, essentially through repeated
applications of Stinespring’s dilation theorem. The conceptual picture of the
RI scheme is that of a quantum system which, for short times, repeatedly
interacts with an environment comprised of small discrete subsystems or
ancilla, that are disjoint. This means that the subsystems comprising the
environment do not interact with each other. In addition, once interacting
with the quantum system of interest, the subsystems travel off to infinity and
never see the system again. It is intuitive that such a construction leads to
a Markov process, and it turns out that under certain conditions on the
coupling energies, this scheme approximates an evolution generated by the
Lindblad QME discussed in the previous Section [11, 10]. Studies regarding a
more general RI form including interacting ancilla and recurring interactions
have also been studied, resulting in non-Markovian dynamics, as expected [22].
However, this generalization is beyond the scope of our work. We now present
the formalism for the Markovian repeated interactions scheme described above.
The system Hamiltonian $H_{0}$ lives in the Hilbert space $\mathscr{H}_{0}$,
and the $n$th subsystem of the environment lives in the smaller Hilbert space
$\mathscr{E}_{n}$. We represent the entirety of the environment space as the
tensor product of all said subsystems:

$\mathscr{H}_{E}=\bigotimes_{n=1}^{m}\mathscr{E}_{n},$ (2)

and the total Hilbert space of the problem is

$\mathscr{H}=\mathscr{H}_{0}\otimes\mathscr{H}_{E}.$ (3)

We write the total Hamiltonian $H$ in the following way,

$\displaystyle H=H_{0}+\sum_{n=1}^{m}(H_{E_{n}}+\lambda H_{I_{n}}).$ (4)

The operators act as identities outside the support of their respected
subspace, such that it is implied that $H_{0}\equiv
H_{0}\otimes\openone_{\mathscr{H}_{E}}$, and so on. The interaction
Hamiltonian takes a tensor product form: $H_{n}=V^{\dagger}\otimes
a_{n}+h.c.$, which is required for the repeated interaction formalism;
$\lambda$ is the system-bath coupling parameter. Here, $V$ is some general
system operator, and $a_{n}$ is an annihilation operator on the $n$th
subsystem. For example, $a_{n}$ could be a bosonic annihilation operator or
the spin lowering operator. This total Hamiltonian $H$ generates the repeated
interactions map $\mathcal{R}_{n}$ according to the following definition:

######  Definition 1 (Repeated Interactions Map).

Given a system density operator $\rho_{S}$ of $\dim=d$ that evolves under the
free Hamiltonian $H_{0}$, and a discrete environment comprised of $m$
subsystems $\rho_{E_{n}}$ of $\dim=2$ that evolve freely under $H_{E_{n}}$ and
interact with $\rho_{S}$ via an interaction Hamiltonian $H_{I_{n}}$ with
coupling $\lambda$, then the finite time reduced dynamics of the system
$\rho_{S}(t_{0})\rightarrow\rho_{S}(t_{0}+\tau)$ is generated by the following
repeated interactions map $\mathcal{R}_{n}$:

$\mathcal{R}_{n}(\rho_{S})=\Tr_{\mathscr{E}_{n}}\left(e^{-i\tau(H_{0}+H_{E_{n}}+\lambda
H_{I_{n}})}\rho_{S}\otimes\rho_{E_{n}}e^{i\tau(H_{0}+H_{E_{n}}+\lambda
H_{I_{n}})}\right),$ (5)

with interaction time $\tau=t/m$. The entire reduced dynamics for a time $t$
yielding $\rho(t_{0}+t)$ is produced by $m$ applications of the map.

Note that the map here is applied only to $\rho_{S}$, and part of the action
of the map is to tensor on the subsystem $\rho_{E_{n}}$. The repeated
interactions map is necessarily applied to the joint system-environment space,
but due to the trace it yields the reduced dynamics of the system as required.
These maps are Markovian by construction, and given the simplest case where
each of the $n$ subsystems are identical, we also have the semi-group
property:

$\mathcal{R}^{\circ m}=\mathcal{R}^{\circ(m-q)}\circ\mathcal{R}^{\circ q},$
(6)

for $q<m$, which shows that any composition of RI maps is also an RI map.
These maps do not form a complete group as their inverses are not guaranteed
to exist [10]. In this work, to maintain generality, we do not assume the bath
qubits to be identical, and subscript the bath’s subsystems by $n$. Another
important property is that the discretization of the dynamics in this way
yields no Trotter error. This is because $[H_{0},H_{E_{n}}]=0$, as is clear
from above, and the interaction $H_{I_{n}}$ is only turned on for the
interaction time $\tau$, similar to a delta kick Hamiltonian. For more
properties and applications of RI maps, consult the following reviews [10,
11]. The most important feature discussed in both, is the correspondence to
Lindblad dynamics. Below, we give a similar derivation to what can be found in
[10, 11]. However, we provide a more rigorous treatment than what is done in
Refs. [10, 21], and show a simpler derivation than that contained in Ref.
[11].

##  IV Lindbladian Derivation

Let us now show that in a specified limit, the repeated interactions scheme
converges to Lindbladian dynamics in the limit of infinitesimally small
interaction time. This is provided more so as a demonstration and is not a
novel result. Similar work can be found in Refs. [11, 10, 21]. However, our
derivation tries to achieve the goal of being both pedagogical and rigorous.
We will work in the regime where $\lambda^{2}\tau=1$, and will study the limit
$\tau\rightarrow 0$. Given the nature of this limit, we can think of the
repeated interactions as strongly coupling the system to the bath’s subsystem
for very short times. This also provides a somewhat interesting physical
picture, given that the Lindblad equation is derived for weak coupling.

To derive the Lindblad equation, we write the map in Liouvillian super-
operator form and expand for small interaction times,

$\displaystyle\mathcal{R}_{n}(\rho_{S})=\Tr_{\mathscr{E}_{n}}\left(\exp(-i\tau[(H_{0}+H_{E_{n}}+\lambda
H_{I_{n}}),~{}\cdot~{}])(\rho_{S}\otimes\rho_{E_{n}})\right),$ (7)

or we define a super operator $\mathcal{H}(\rho)=-i[(H_{0}+H_{E_{n}}+\lambda
H_{I_{n}}),\rho]$ and equivalently write

$\displaystyle\mathcal{R}_{n}(\rho_{S})=\Tr_{\mathscr{E}_{n}}\left(e^{\tau\mathcal{H}}(\rho_{S}\otimes\rho_{E_{n}})\right).$
(8)

For each application, this map effectively performs a dilation by tensoring on
an environmental subsystem, followed by a unitary evolution and a reduction.
The map is CPTP by Stinespring’s dilation theorem. For a single application of
the map lets expand to second order in the interaction time

$\displaystyle\mathcal{R}_{n}(\rho_{S})$
$\displaystyle=\Tr_{\mathscr{E}_{n}}\left(e^{\tau\mathcal{H}}(\rho_{S}\otimes\rho_{E_{n}})\right)$
$\displaystyle=\Tr_{\mathscr{E}_{n}}\Biggr{(}\rho_{S}\otimes\rho_{E_{n}}-i\tau[(H_{0}+H_{E_{n}}+\lambda
H_{I_{n}}),\rho_{S}\otimes\rho_{E_{n}}]$
$\displaystyle\quad-\frac{\tau^{2}}{2}[(H_{0}+H_{E_{n}}+\lambda
H_{I_{n}}),[(H_{0}+H_{E_{n}}+\lambda
H_{I_{n}}),\rho_{S}\otimes\rho_{E_{n}}]]\Biggr{)}$
$\displaystyle=\rho_{S}-i\tau[H_{0},\rho_{S}]-\Tr_{\mathscr{E}_{n}}\left(\frac{\tau^{2}}{2}[(H_{0}+H_{E_{n}}+\lambda
H_{I_{n}}),[(H_{0}+H_{E_{n}}+\lambda
H_{I_{n}}),\rho_{S}\otimes\rho_{E_{n}}]]\right),$ (9)

In the last step we use the fact that
$\Tr_{E}[H_{E},\rho_{S}\otimes\rho_{E_{n}}]=0$. As well, we can either
redefine the interaction Hamiltonian, or choose the initial environment state
such that $\Tr_{E}[H_{I},\rho_{S}\otimes\rho_{E_{n}}]=0$, as is a standard
choice in the derivation of Lindblad master equations. This leads to

$\displaystyle\mathcal{R}_{n}(\rho_{S})$
$\displaystyle=\rho_{S}-i\tau[H_{0},\rho_{S}]-\frac{\tau^{2}}{2}\Tr_{\mathscr{E}_{n}}\Biggl{(}\lambda^{2}\bigl{[}H_{I_{n}},\bigl{[}H_{I_{n}},\rho_{S}\otimes\rho_{E_{n}}\bigr{]}\bigr{]}+\lambda\bigl{[}(H_{0}+H_{E_{n}}),\bigl{[}H_{I_{n}},\rho_{S}\otimes\rho_{E_{n}}\bigr{]}\bigr{]}$
$\displaystyle\quad+\lambda\bigl{[}H_{I_{n}},\bigl{[}(H_{0}+H_{E_{n}}),\rho_{S}\otimes\rho_{E_{n}}\bigr{]}\bigr{]}+\bigl{[}(H_{0}+H_{E_{n}}),\bigl{[}(H_{0}+H_{E_{n}}),\rho_{S}\otimes\rho_{E_{n}}\bigr{]}\bigr{]}\Biggr{)}$
(10)

By taking $\lambda\rightarrow\frac{1}{\sqrt{\tau}}$ we get

$\displaystyle\widetilde{\mathcal{R}_{n}}(\rho_{S})=\rho_{S}-i\tau[H_{0},\rho_{S}]-\frac{1}{2}\Tr_{\mathscr{E}_{n}}\Biggl{(}\tau[H_{I_{n}},[H_{I_{n}},\rho_{S}\otimes\rho_{E_{n}}]]+\tau^{3/2}[(H_{0}+H_{E_{n}}),[H_{I_{n}},\rho_{S}\otimes\rho_{E_{n}}]]$
$\displaystyle+\tau^{3/2}[H_{I_{n}},[(H_{0}+H_{E_{n}}),\rho_{S}\otimes\rho_{E_{n}}]]+\tau^{2}[(H_{0}+H_{E_{n}}),[(H_{0}+H_{E_{n}}),\rho_{S}\otimes\rho_{E_{n}}]]\Biggr{)},$
(11)

where we use a tilde to indicate that this limit has been enforced. This
replacement of variables is thoroughly discussed in Ref. [11]. Rearranging and
taking the limit as $\tau\rightarrow 0$ leads to

$\displaystyle\lim_{\tau\rightarrow
0}\left(\frac{\widetilde{\mathcal{R}_{n}}(\rho_{S})-\rho_{S}}{\tau}\right)$
$\displaystyle=-i[H_{0},\rho_{S}]-\frac{1}{2}\Tr_{\mathscr{E}_{n}}\left([H_{I_{n}},[H_{I_{n}},\rho_{S}\otimes\rho_{E_{n}}]]\right)$
$\displaystyle=-i[H_{0},\rho_{S}]+\Tr_{\mathscr{E}_{n}}\left(H_{I_{n}}\rho_{S}\otimes\rho_{E_{n}}H_{I_{n}}-\frac{1}{2}\\{H_{I_{n}}H_{I_{n}},\rho_{S}\otimes\rho_{E_{n}}\\}\right),$
(12)

which is in Lindblad form. To then eliminate the trace and write this solely
in terms of the system Hamiltonian, we can use the full form of
$H_{I_{n}}=V\otimes a^{\dagger}+h.c.$. Following [11], if we chose the
environment to be in some thermal state at inverse temperature $\beta$, we
obtain the following

$\displaystyle\lim_{\tau\rightarrow
0}\left(\frac{\mathcal{R}_{n}(\rho_{S})-\rho_{S}}{\tau}\right)=-i[H_{0},\rho_{S}]+\sum_{j=1}^{2}\left(L_{j}\rho_{S}L_{j}^{\dagger}-\frac{1}{2}\\{L_{j}^{\dagger}L_{j},\rho_{S}\\}\right),$
(13)

where $L_{j}=V\sqrt{z}$ for $j=1$, $L_{j}=V^{\dagger}\sqrt{z-1}$ for $j=2$,
and $z=\nicefrac{{\Tr_{\mathscr{E}_{n}}\left(a_{n}^{\dagger}\exp(-\beta
a_{n}^{\dagger}a_{n})a_{n}\right)}}{{Z(\beta)}}$. Here, the $a^{\dagger}$ and
$a$ operators are the creation and annihilation operators, and $Z$ is the
partition function. Later we will consider an implementation where the
environment is comprised of single qubit thermal states, and as a result, we
do not have to worry about exchange statistics. For this reason, $a^{\dagger}$
and $a$ are sometimes referred to as qubit creation and annihilation
operators. With these definitions, the full RI map, the $j$ index will grow to
have $2m$ terms for the entire repeated interaction map, where $m$ is the
number of subsystems. For the simple case of a single interaction with the RI
map, we have therefore shown how the Lindblad equation can be derived.

##  V Error Bounding in the Limiting RI case

In this Section we prove an error bound between $\nu$ applications of a RI map
composed of $m$ interactions, and Lindbladian dynamics. We work here only in
the limiting case $\lambda\rightarrow\frac{1}{\sqrt{\tau}}$ of the RI map
defined in Equation 1. We can attain equivalent results by using some
constant: $\lambda^{2}\tau\rightarrow\gamma$. However, we allow this freedom
of the Lindblad jump coefficients to be taken into account through the
spectral norm of the Interaction Hamiltonian. We work in this limit, which can
be regards as a replacement of variables, to effectively control the error of
the approximation. Working with a general $\lambda$ yields error in the first
non-trivial order, which is unfavourable. To see this argument worked out, we
refer readers to Appendix A. As in the previous Section, we use a tilde to
remind the reader that this limiting case of the RI map is being enforced, and
write a single application of the map in the following way:
$\widetilde{\mathcal{R}_{n}}(\rho_{S})$. This conditioned is understood to be
enforced in the Taylor series expansion of the super-operator, as was also
considered in the above derivation. For generality, we leave interactions and
the free Hamiltonians of the environment ancilla subscripted with $n$. In this
way, the bath can be prepared as an ensemble of single qubit thermal states
with different temperatures $\beta$, if the user wishes. Our next objective is
to show convergence of the RI map to a general Lindblad QME. The RI map we use
to converge to this dynamics, with an arbitrary accuracy, takes the following
form

$\widetilde{\mathcal{R}}=\Bigr{(}\bigcirc_{n=1}^{m}\widetilde{\mathcal{R}_{n}}(\nicefrac{{t}}{{\nu}})\Bigr{)}^{\circ\nu},$
(14)

where $\bigcirc$ indicates an indexed multi-composition and
$\widetilde{\mathcal{R}_{n}}$ is given by Eq. (9) upon substitution of
$\lambda^{2}\tau=1$. Intuitively, an immediate issue arises, in that the
general Lindblad picture consists of the system simultaneously coupled to many
degrees of freedom, undergoing evolution coupled all subsystems for a time
$t$, while the RI map simulates the free evolution for a time $t$, and evolves
for a time $\nicefrac{{t}}{{m}}$ under each of the $m$ ancilla. To deal with
this, we simply re-normalize the free Hamiltonian of the system
$H_{0}\rightarrow\frac{1}{m}H_{0}$, which allows us to converge to the
following general Lindbladian,

$\mathcal{L}(\rho_{S})=-i[H_{0},\rho]-\frac{1}{2}\sum_{k=1}^{m}\Tr_{\mathscr{H}_{E}}[H_{I_{k}},[H_{I_{k}},\rho_{S}\otimes\rho_{E_{k}}]].$
(15)

We can use the trace to extract the more natural jump operators later, as was
done in the derivation of Equation (13). With this in mind, we provide an
error bound in the following Theorem:

######  Theorem 2 (RI-Lindblad Correspondence Error).

Given an RI map from Definition 1 and a Lindblad Equation 15, then the induced
Schatten 1-norm of their difference has the following upper bound:

$\displaystyle\norm{e^{t\mathcal{L}}-\Bigr{(}\bigcirc_{n=1}^{m}\widetilde{\mathcal{R}_{n}}(\nicefrac{{t}}{{\nu}})\Bigr{)}^{\circ\nu}}_{1\rightarrow
1}$ $\displaystyle\in
O\Biggr{(}\frac{t^{2}}{\nu}\Biggr{(}\norm{\mathcal{L}}_{1\rightarrow
1}^{2}+m\max_{n}\lambda_{n}^{4}\Biggr{)}\Biggr{)},$ (16)

where $\nu$ is the number of applications of the RI map, $t$ the total
simulation time, and $m$ the number of subsystems in the bath.

###### Proof.

Our goal is to show that in terms of the time scaling we have error that goes
like $O(\tau^{2})$, in accordance with other analyses in the literature [11,
21].

We begin by defining the following super-operator to represent the dynamics
generated by the RI map:

$\widetilde{\mathcal{H}_{n}}(\rho)=\left[\frac{1}{m}H_{0}+H_{E_{n}}+\frac{1}{\sqrt{\tau}}H_{I_{n}},\rho\right].$
(17)

This superoperator already accounts for the condition
$\lambda^{2}\tau\rightarrow 1$ and renormalizes the free Hamiltonian by $m$ to
capture the correct dynamics.

$\displaystyle\norm{e^{t\mathcal{L}}-\Bigr{(}\bigcirc_{n=1}^{m}\widetilde{\mathcal{R}_{n}}(\nicefrac{{t}}{{\nu}})\Bigr{)}^{\circ\nu}}_{1\rightarrow
1}=\norm{(e^{\frac{t}{\nu}\mathcal{L}})^{\circ\nu}-\Bigr{(}\bigcirc_{n=1}^{m}\widetilde{\mathcal{R}_{n}}(\nicefrac{{t}}{{\nu}})\Bigr{)}^{\circ\nu}}_{1\rightarrow
1}\leq\nu\norm{e^{\frac{t}{\nu}\mathcal{L}}-\bigcirc_{n=1}^{m}\widetilde{\mathcal{R}_{n}}(\nicefrac{{t}}{{\nu}})}_{1\rightarrow
1}$ $\displaystyle=\nu\max_{\norm{\rho_{S}}_{1}\leq
1}\Biggr{|}\Biggr{|}\Biggr{(}\openone-i\frac{t}{\nu}\left[\frac{1}{m}H_{0},\cdot\right]-\frac{t}{2\nu}\sum_{k=1}^{m}\Tr_{\mathscr{H}_{E}}[H_{I_{k}},[H_{I_{k}},\cdot\otimes\rho_{E_{k}}]]+\sum_{j=2}^{\infty}\frac{\mathcal{L}^{j}t^{j}}{\nu^{j}j!}\Biggl{)}[\rho_{S}]$
$\displaystyle-\Tr_{\mathscr{H}_{E}}\prod_{n=1}^{m}\Biggr{(}\openone-i\frac{t}{\nu}\left[\frac{1}{m}H_{0}+H_{E_{n}}+\sqrt{\frac{\nu}{t}}H_{I_{n}},\cdot\otimes\rho_{E_{n}}\right]-\frac{t}{2\nu}[H_{I_{n}},[H_{I_{n}},\cdot\otimes\rho_{E_{n}}]]$
$\displaystyle-\frac{t^{3/2}}{2\nu^{3/2}}\left(\left[\left(\frac{1}{m}H_{0}+H_{E_{n}}\right),\left[H_{I_{n}},\cdot\otimes\rho_{E_{n}}\right]\right]+\left[H_{I_{n}},\left[\left(\frac{1}{m}H_{0}+H_{E_{n}}\right),\cdot\otimes\rho_{E_{n}}\right]\right]\right)$
$\displaystyle-\frac{t^{2}}{2\nu^{2}}\left[\left(\frac{1}{m}H_{0}+H_{E_{n}}\right),\left[\left(\frac{1}{m}H_{0}+H_{E_{n}}\right),\cdot\otimes\rho_{E_{n}}\right]\right]+\sum_{l=3}^{\infty}\frac{\widetilde{\mathcal{H}_{n}}^{l}t^{l}}{\nu^{l}l!}\Biggl{)}[\rho_{S}]\Biggl{|}\Biggl{|}_{1}.$
(18)

To ensure convergence of the Taylor series we make necessary smallness
assumptions on the Hamiltonian and the interaction time:
$\\{\norm{H_{0}}t/\nu,\norm{H_{E}}t/\nu\\}\in o(1)$ and
$\norm{H_{I}}\sqrt{t/\nu}\in o(1)$. Then, by expanding the product above and
tracking only terms $\in O\left(t^{2}/\nu^{2}\right)$ we obtain

$\displaystyle\norm{e^{t\mathcal{L}}-\Bigr{(}\bigcirc_{n=1}^{m}\widetilde{\mathcal{R}_{n}}(\nicefrac{{t}}{{\nu}})\Bigr{)}^{\circ\nu}}_{1\rightarrow
1}$ $\displaystyle\in O\Biggr{(}\nu\Biggr{(}m\max_{n,\norm{\rho_{S}}_{1}\leq
1}\bigr{|}\bigr{|}\Tr_{\mathscr{E}_{n}}\frac{t^{2}}{2\nu^{2}}[\frac{1}{m}H_{0},[\frac{1}{m}H_{0},\cdot\otimes\rho_{E_{n}}]]\bigl{|}\bigl{|}$
$\displaystyle+\norm{\sum_{j=2}^{\infty}\frac{\mathcal{L}^{j}t^{j}}{\nu^{j}j!}}_{1\rightarrow
1}+m\max_{n}\norm{\sum_{l=3}^{\infty}\frac{\widetilde{\mathcal{H}_{n}}^{l}t^{l}}{\nu^{l}l!}}_{1\rightarrow
1}\Biggl{)}\Biggr{)}$ (19)

This bound was achieved due to the following simplifications: First, we
collected from the first line of Eq. (18) the product of each order of the map
with the identity, as these give leading orders in all terms. Then, note that
the partial trace commutes with the free Hamiltonian $H_{0}$. We also recall
that the Lindbladian is derived with the property
$\Tr_{\mathscr{E}_{n}}[H_{I},\rho_{S}\otimes\rho_{E_{n}}]=0$, meaning that the
$\sqrt{\frac{\nu}{t}}\Tr_{\mathscr{E}_{n}}[H_{I_{n}},\cdot\otimes\rho_{E_{n}}]$,
and all terms they multiply necessarily vanish. Consequentially, terms that
scale like $O(\frac{t^{1/2}}{\nu^{1/2}})$ and $O(\frac{t^{3/2}}{\nu^{3/2}})$
vanish. Note that the sum in the final term in Equation 19 contains a single
term that goes like $O(\frac{t^{3/2}}{\nu^{3/2}})$ due to the inverse square
root that multiplies the interaction Hamiltonian. This term involves a triple
commutator like so: $\tau^{3/2}[H_{I_{n}},[H_{I_{n}},[H_{I_{n}},\rho]]]$,
however, given the restrictions that impose
$\Tr_{\mathscr{E}_{n}}[H_{I},\rho]=0$, it can be argued that the partial trace
over $\mathscr{E}_{n}$ for any odd order commutator is 0. To make this
argument clear, if we are to choose that $\rho_{E_{n}}$ is a thermal state, as
in the case of our prior derivation of Equation 13, and $H_{I_{n}}$ also has
the form defined in (13), then any odd combination of $a$ and $a^{\dagger}$
makes this necessarily 0, which is a well understood property of creation-
annihilation (or raising lowering) operators.

We now proceed with the bound (19), and will revisit
$\widetilde{\mathcal{H}}_{n}^{3}$ at the end of the proof. We proceed by using
tail bounds for the two exponential series,

$\displaystyle\norm{e^{t\mathcal{L}}-\Bigr{(}\bigcirc_{n=1}^{m}\widetilde{\mathcal{R}_{n}}(\nicefrac{{t}}{{\nu}})\Bigr{)}^{\circ\nu}}_{1\rightarrow
1}$ $\displaystyle\in
O\Biggr{(}\nu\Biggr{(}\frac{t^{2}}{m\nu^{2}}\norm{H_{0}}^{2}+\frac{\norm{\mathcal{L}}_{1\rightarrow
1}^{2}t^{2}}{\nu^{2}}e^{\frac{t\norm{\mathcal{L}}_{1\rightarrow 1}}{\nu}}$
$\displaystyle\quad+m\max_{n}\Bigr{(}\frac{t^{3}\norm{\widetilde{\mathcal{H}_{n}}^{3}}_{1\rightarrow
1}}{\nu^{3}}+\frac{t^{4}\norm{\widetilde{\mathcal{H}_{n}}}_{1\rightarrow
1}^{4}}{\nu^{4}}e^{\frac{t\norm{\widetilde{\mathcal{H}_{n}}}_{1\rightarrow
1}}{\nu}}\Bigr{)}\Biggr{)}\Biggr{)}.$ (20)

By making the requirements that $\\{\frac{t\norm{\mathcal{L}}_{1\rightarrow
1}}{\nu},~{}\frac{t\norm{\widetilde{\mathcal{H}_{n}}}_{1\rightarrow
1}}{\nu}\\}\in[0,\ln{2}]~{}\forall~{}n$, we linearize the bound into the
following

$\displaystyle\norm{e^{t\mathcal{L}}-\Bigr{(}\bigcirc_{n=1}^{m}\widetilde{\mathcal{R}_{n}}(\nicefrac{{t}}{{\nu}})\Bigr{)}^{\circ\nu}}_{1\rightarrow
1}$ $\displaystyle\in
O\Biggr{(}\frac{t^{2}}{m\nu}\norm{H_{0}}^{2}+\frac{\norm{\mathcal{L}}_{1\rightarrow
1}^{2}t^{2}}{\nu}+m\max_{n}\Bigr{(}\frac{t^{3}\norm{\widetilde{\mathcal{H}_{n}}^{3}}_{1\rightarrow
1}}{\nu^{2}}+\frac{t^{4}\norm{\widetilde{\mathcal{H}_{n}}}_{1\rightarrow
1}^{4}}{\nu^{3}}\Bigr{)}\Biggr{)}.$ (21)

Now, as a result of prior arguments about odd order commutators with
$H_{I_{n}}$, we have the following

$\displaystyle\frac{t^{3}}{\nu^{2}}\norm{\widetilde{\mathcal{H}_{n}}^{3}}_{1\rightarrow
1}$ $\displaystyle\in
O\Bigr{(}\frac{t^{3}}{m^{3}\nu^{2}}\norm{H_{0}}^{3}+\frac{t^{2}}{\nu}\norm{H_{E_{n}}}\norm{H_{I_{n}}}^{2}+\frac{t^{2}}{m\nu}\norm{H_{0}}\norm{H_{I_{n}}}^{2}\Bigr{)}$
$\displaystyle\in
O\Bigr{(}\frac{t^{3}}{m^{3}\nu^{2}}\norm{H_{0}}^{3}+\lambda_{n}^{3}\frac{t^{2}}{\nu}\Bigr{)},$
(22)

where
$\lambda_{n}\equiv\max\\{\norm{H_{0}},\norm{H_{E_{n}}},\norm{H_{I_{n}}}\\}$.
Now the first term in the above equation is clearly upper bounded by
$\frac{t^{2}}{m\nu}\norm{H_{0}}^{2}$, so we now have

$\displaystyle\norm{e^{t\mathcal{L}}-\Bigr{(}\bigcirc_{n=1}^{m}\widetilde{\mathcal{R}_{n}}(\nicefrac{{t}}{{\nu}})\Bigr{)}^{\circ\nu}}_{1\rightarrow
1}$ $\displaystyle\in
O\Biggr{(}\frac{t^{2}}{m\nu}\norm{H_{0}}^{2}+\frac{\norm{\mathcal{L}}_{1\rightarrow
1}^{2}t^{2}}{\nu}+m\max_{n}\Bigr{(}\lambda_{n}^{3}\frac{t^{2}}{\nu}+\frac{t^{4}\norm{\widetilde{\mathcal{H}_{n}}}_{1\rightarrow
1}^{4}}{\nu^{3}}\Bigr{)}\Biggr{)}.$ (23)

Now we deal with the final term

$\displaystyle\frac{t^{4}\norm{\widetilde{\mathcal{H}_{n}}}_{1\rightarrow
1}^{4}}{\nu^{3}}$
$\displaystyle\leq\nu\Bigr{(}\frac{t}{\nu}\norm{H_{0}}+\frac{t}{\nu}\norm{H_{E_{n}}}+\sqrt{\frac{t}{\nu}}\norm{H_{I_{n}}}\Bigr{)}^{4}$
$\displaystyle\in
O\Biggr{(}\nu\Bigr{(}\sqrt{\frac{t}{\nu}}\norm{H_{0}}+\sqrt{\frac{t}{\nu}}\norm{H_{E_{n}}}+\sqrt{\frac{t}{\nu}}\norm{H_{I_{n}}}\Bigr{)}^{4}\Biggr{)}$
$\displaystyle\in O\Bigr{(}\frac{t^{2}}{\nu}\lambda_{n}^{4}\Bigr{)}.$ (24)

To get the second line we recalled that $\tau=t/\nu$ and that we are
interested in the asymptotic limit where $\tau\rightarrow 0$. As such,
$\tau<\sqrt{\tau}$. We are now left with

$\norm{e^{t\mathcal{L}}-\Bigr{(}\bigcirc_{n=1}^{m}\widetilde{\mathcal{R}_{n}}(\nicefrac{{t}}{{\nu}})\Bigr{)}^{\circ\nu}}_{1\rightarrow
1}\in
O\Biggr{(}\frac{t^{2}}{m\nu}\norm{H_{0}}^{2}+\frac{\norm{\mathcal{L}}_{1\rightarrow
1}^{2}t^{2}}{\nu}+m\max_{n}\Bigr{(}\lambda_{n}^{4}\frac{t^{2}}{\nu}\Bigr{)}\Biggr{)}.$
(25)

Noting that the first term in (25) is asymptotically upper bounded by
$m\max_{n}\lambda_{n}^{4}\frac{t^{2}}{\nu}$ completes the proof. ∎

Some useful corollaries immediately follow from the Theorem above:

######  Corollary 2.1 (Number of Interactions).

The minimum necessary number of iterations $\nu$ of an RI map composed of $m$
interactions, written as in Equation (14), required to simulate the time
evolution of a Lindbladian within error $\epsilon$ has the following upper-
bound:

$\nu\in
O\Biggr{(}\frac{t^{2}}{\epsilon}\Bigr{(}\norm{\mathcal{L}}_{1\rightarrow
1}^{2}+m\max_{n}\lambda_{n}^{4}\Bigr{)}\Biggr{)}.$ (26)

Here $\lambda_{n}$ is the strength of the interaction used in the $n^{\rm th}$
interaction step. This then implies that the total number of RI maps
$\kappa=m\nu$ that needs to be implemented to realize dynamics $\epsilon$
close to that generated by $e^{t\mathcal{L}}$ has the following upper bound:

$\kappa\in
O\Biggr{(}\frac{t^{2}}{\epsilon}\Bigr{(}m\norm{\mathcal{L}}_{1\rightarrow
1}^{2}+m^{2}\max_{n}\lambda_{n}^{4}\Bigr{)}\Biggr{)}.$ (27)

As a direct consequence, this also upper bounds the number of ancilla qubits
used in the simulation strategy outlined later in Figure 1.

Furthermore, one of the insights that can be gained through the analysis
contained in Theorem 2 is that in terms of the time scaling of the error, the
higher orders do not necessarily contribute to converging towards Lindblad
dynamics. Part of the reason is that analyzing the higher order terms is
extremely difficult, and potentially intractable using methods from this
analysis. In this way, the methods are both poorly sorted to track the higher
order contributions, as well as there is no obvious structure or symmetries to
indicate higher order convergence. This leads us to conclude that the RI maps
defined in this work likely do not lead to higher order integrators, which
given the prior discussion in the Introduction, is not of an immense concern.
However, since a large part of the role of $\nu$ is simply to control the
higher order terms, it is possible that a cost effective simulation may simply
involve not simulating them. We obtain the following Corollary:

######  Corollary 2.2 (Second order Correspondence Error).

If an RI map defined as in Equation 14 is implemented up to second order in
the interaction time $\tau$, then the induced Schatten $1\rightarrow 1$ norm
of the difference between the RI map and the dynamics generated by the
Lindblad equation $e^{t\mathcal{L}}$ has the following upper bound:

$O\Biggr{(}\frac{t^{2}}{m\nu}\norm{H_{0}}^{2}+\frac{\norm{\mathcal{L}}_{1\rightarrow
1}^{2}t^{2}}{\nu}\Biggr{)},$ (28)

which can be seen through noting that the $\lambda_{n}$ terms are derived only
from higher order expansions.

Given the superior error scaling over Theorem 2, it may be possible to compose
more efficient quantum simulations. However, this truncation leads to a non-
unitary map. Implementing this map could therefore be done in a manner that
combines strategies proposed in Ref. [12]. Here, it is argued that a
Hamiltonian embedding of the Lindbladian dynamics with repeated trace-out is
inefficient compared to the method of using Kraus maps master equations and
then implementing them in a purified space using a variant of LCU for
isometries. If we were to instead implement repeated interactions to second
order, this could be seen as an intermediate approach between these schemes,
of which can be made much more compact in the number of qubits in comparison
to the Hamiltonian approach in [12]. However, this then leads to an algorithm
where $\delta\neq 1$, and would require a total of $\kappa$ rounds of
oblivious amplitude amplification. Therefore, it is unclear how this strategy
would compare, a topic of a separate investigation. In this work, we continue
to focus on unitary implementations of the RI map using iterative Hamiltonian
simulations.

To summarize, this section has achieved the goal of bounding the total error
between the dynamics generated by a Lindblad equation, and the dynamics
realized by $\nu$ applications of an RI map composed of $m$ interactions with
an environment. As a corollary, we have derived a bound on the minimum number
of system bath interactions to use an RI map to simulate $e^{t\mathcal{L}}$ up
to precision $\epsilon$. We can now proceed with strategies from quantum
computing to analyze the complexity of implementing said repeated
interactions, in order to derive the overall complexity of a quantum algorithm
based on repeated interactions to simulate Lindbladian Markovian dynamics.

##  VI RI Implementation with Hamiltonian Simulation

In this Section, we write down a unitary form of the map from the previous
Section, and show that using the principle of deferred measurement, one can
make the entire process unitary prior to tracing out the entire environment.
We therefore apply two approaches. In the first one, we iteratively use
Hamiltonian simulation strategies to perform the trace, which allows one to
bound the number of necessary ancilla qubits by a constant. The other approach
preserves unitarity for the entire evolution time, but then requires a number
of ancilla that scales polynomially with the simulation time $t$. Here, we use
Qubitization and Trotter-Suzuki formulas to implement the RI maps. We utilize
these algorithms given that they can be performed without amplitude
amplification, as both of them have a probability of success $\delta=1$. In
addition, Qubitization is utilized to achieve optimal asymptotics, whereas
Trotter-Suzuki formulas are chosen for their implementation simplicity, as
well as to take advantage of the commutator structure of the problem. The
bounds that follow will necessarily depend on $\nu$, however, physically
speaking there may exist cases where the RI maps actually provide a more
accurate description of the Markovian dynamics than the Lindblad equation, as
is likely the case of the Maser discussed in Section I and expanded on in Ref.
[11]. In these cases, the complexity of $\nu$ given in Corollary 2.1 need not
apply, and may instead be interpreted as constant. For these reasons, we state
our bounds in terms of $\nu$ rather than hashing out the full time complexity,
which can simply be obtained by combining the theorems that follow with
Corollary 2.1.

###  VI.1 Unitary RI Map and Environment Preparation

The map $\widetilde{\mathcal{R}}$ derived in the previous section required a
Taylor expansion in the interaction time, and then a replacement of the
coupling parameter $\lambda\rightarrow 1/\sqrt{\tau}$. This procedure leaves
us with a map we can implement to a desired order. However, the map does not
have a clear functional form, which can make it difficult to implement with
algorithms that are not LCU. The sum form of the map also does not give
insight on which parts of the evolution are unitary. Here, we show an RI map
that can be written with only unitary operators that produces the same output.

$\Qcircuit@C=1em@R=1em{\lstick{\rho_{S}}&\qw/^{\log
d}\qw\multigate{1}{U_{1}}\qw\qw\qw\qw\gate{U_{2}}\qw\qw\cdots~{}\qw\gate{U_{m}}\qwa~{}~{}{\approx\exp(\mathcal{L}(\rho_{S}))}\\\
\lstick{\rho_{E_{1}}}\qw\qw\ghost{U_{1}}\qw\qw\qw\qw\qw\qwx\qw\qw\cdots\qw\qw\qwx\meter\\\
\lstick{\rho_{E_{2}}}\qw\qw\qw\qw\qw\qw\qw\gate{U_{2}}\qwx\qw\qw\cdots\qw\qw\qwx\meter\\\
\vdots~{}~{}~{}~{}\vdots~{}~{}~{}~{}\vdots~{}~{}~{}~{}\qwx\vdots\\\
\lstick{\rho_{E_{m}}}\qw\qw\qw\qw\qw\qw\qw\qw\qw\qw\cdots\qw\gate{U_{m}}\qwx\meter}$
Figure 1: This illustration presents the repeated interactions implementation
utilizing the principle of deferred measurement. In this way, the entire
process is unitary, availing the use of other quantum algorithms prior to the
final measurement, such as amplitude application. The cost of this feature,
however, is $\kappa$ qubits, which is upper bounded in Corollary 2.1.
$\Qcircuit@C=1em@R=1em{\lstick{\rho_{S}}&\qw/^{\log
d}\qw\multigate{1}{U_{1}}\qw\qw\qw\qw\multigate{1}{U_{2}}\qw\qw\cdots~{}\qw\qw\multigate{1}{U_{m}}\qwa~{}~{}{\approx\exp(\mathcal{L}(\rho_{S}))}\\\
\lstick{\rho_{E_{1}}}\qw\qw\ghost{U_{1}}\meter~{}~{}\lstick{\rho_{E_{2}}}\ghost{U_{2}}\meter~{}\cdots~{}~{}\lstick{\rho_{E_{m}}}\ghost{U_{m}}\meter}$
Figure 2: Illustrated above is the iterative way of applying the repeated
interactions map. Given that the complexity of implementing $U_{n}$ is far
greater than preparing the environmental states $\rho_{E_{n}}$, we can utilize
the same register to iteratively prepare our environment and then apply the
interaction unitary with our algorithm of choice, to achieve a simulation with
a constant number of qubits.

Suppose we have the unitary $U_{n}(\tau)=\text{exp}(-i(H_{0}+H_{E_{n}})\tau-
iH_{I_{n}}\sqrt{\tau})$. Lets consider the channel that implements this
unitary and an ancillary system to yield the evolution of the reduced
dynamics:
$\hat{\mathcal{R}}_{n}[\cdot]=\Tr_{\mathscr{E}_{n}}U_{n}(\cdot\otimes\rho_{E_{n}})U_{n}^{\dagger}$.
As before, we also define a super-operator, which in this case is time
dependent:
$\hat{\mathcal{H}}_{n}(\rho)=-i[\tau(H_{0}+H_{E_{n}})+\sqrt{\tau}H_{I_{n}},\rho]$.
As before, using the adjoint endomorphism of the Lie group, we can write the
evolution super-operator $e^{\hat{\mathcal{H}}_{n}}$. For the a RI map with
$m$ interactions/ancilla this yields the following:

$\bigcirc_{n=1}^{m}\hat{\mathcal{R}}_{n}(\rho_{S}\otimes\rho_{E_{n}})=\Tr_{\mathscr{H}_{E}}\prod_{n=1}^{m}e^{\hat{\mathcal{H}}_{n}}(\rho_{S}\otimes\rho_{E_{n}})$
(29)

Through writing this as a Taylor series of commutators, one can see that this
leads to the same error scaling as above (for all orders). This then
introduces a subtlety in the algorithmic analysis going forward that is not
present in standard Hamiltonian evolution. We have introduced different
evolution time scaling for the free and interaction Hamiltonians, ie. without
assumptions of the functional scaling of $H_{I}/H_{0}$ or carrying out a more
careful analysis with respect to $\tau$, one cannot absorb the spectral norm
of $H$ into $t$ and simply apply state of the art bounds. The analyses
performed in the following sections reflect this, and carefully track powers
of $t$ to modify error bounds from the literature.

Observing Equation 29, we observe that we can perform the RI implementation in
one of two ways. We can either perform intermittent trace-out measurements, or
performed the trace over the entire environment at the end. The latter is more
useful when non-unitary algorithms are used to implement each of the $U_{n}$s
such that amplitude amplification can be used at the end. This method is
illustrated in Figure 1. Given that we have $m$-ancilla qubits to generate the
interactions in this case, the number of qubits is naturally upper bounded by
Corollary 2.1. The other approach is a prepare trace-out refresh strategy, in
which the environment is measured (with the information discarded), and the
thermal state is once again prepared or reset. Given that part of the
definition of the RI map consists of an environment or ancilla qubit that is
tensored on to the system with each application, we must also consider the
cost of preparing such a state. Since we are using single qubit thermal
states, these can be prepared very efficiently, with complexity far favourable
to that of implementing $U_{n}$. In this way, they can be prepared on the fly,
without the need to pre-compute and store ancilla. To show this, we can use
Givens rotations strategies to prepare a state vector that is a purification
of the thermal state density matrix $\rho_{E_{n}}$ [23]. For an $\eta$-qubit
state:

$\ket{\psi}=\sum_{j=0}^{2^{\eta}-1}\frac{\sqrt{\alpha_{j}}\ket{j}\ket{j}}{\sqrt{\sum_{k=0}^{2^{\eta}-1}\alpha_{k}}},$
(30)

the preparation has complexity $O(\eta 2^{\eta}\log(1/\epsilon))$ in general
[24, 23]. However, given we are initializing a single qubit state, this
becomes $O\log(1/\epsilon)$. With this state we simply trace over the second
register:

$\displaystyle\rho_{E_{n}}$ $\displaystyle=\Tr_{2}(\outerproduct{\psi}{\psi})$
(31)
$\displaystyle=\frac{\alpha_{0}\outerproduct{0}{0}+\alpha_{1}\outerproduct{1}{1}}{\alpha_{0}+\alpha_{1}},$
(32)

where the $\alpha_{j}$s are chosen based on the Givens rotations. In this way,
we can efficiently prepare our environment for repeated interactions. For our
purposes, when we have a thermal state $\rho_{E_{n}}=\exp(-\beta
a^{\dagger}a)/\Tr(...)$, and we write $a=X+iY$, our state has the following
coefficients:

$\rho_{E_{n}}=\frac{\outerproduct{0}{0}+e^{-\beta}\outerproduct{1}{1}}{e^{-\beta}+1}.$
(33)

###  VI.2 Iterative Qubitization of the RI Map

Since we have been able to write down the RI map as a unitary process, or a
process that is at least iteratively unitary with intermediate measurement
(depending on where the traces are taken), we can apply state of the art
Hamiltonian simulation algorithms to this problem. For qubitization, we wish
to write the time evolution as a sum of sine and cosine matrix functions. We
can do this here with a redefinition of the interaction Hamiltonian
$H_{I_{n}}^{\tau}\equiv H_{I_{n}}/\sqrt{\tau}$. This leaves us with the
following:

$U_{n}=\text{exp}(-i\tau(H_{0}+H_{E_{n}}+H_{I_{n}}^{\tau}))$ (34)

Now, in order to implement this unitary, we construct a block encoding using
LCU, which requires that we provide a unitary decomposition of each of our
operators here, each of which decomposition has at most polynomially many
terms in the number of qubits. We assume the system Hamiltonian $H_{0}$ has
$l_{0}$ terms and can be written like so (recalling that it has been re-
scaled):

$H_{0}=\sum_{j=1}^{l_{0}}\frac{h_{j}}{m}H_{j},$ (35)

where all phase factors are absorbed into $H_{j}$. We also normalize each of
these terms such that $h_{j}\equiv\norm{H_{j}}$ and after normalization each
$\norm{H_{j}}=1$. This can always be done WLOG. Recall, the interaction
Hamiltonian was given the earlier definition: $H_{I}=V_{n}\otimes
a_{n}^{\dagger}+V_{n}^{\dagger}\otimes a_{n}$. Here we use the re-scaled
$H_{I_{n}}^{\tau}$, and also write the general system operator as a sum of
unitaries, and the creation-annihilation operators in terms of Pauli operators
like so:

$\displaystyle H_{I_{n}}^{\tau}$
$\displaystyle=\frac{1}{\sqrt{\tau}}\Biggr{(}\sum_{p=1}^{l_{I_{n}}}v_{p_{n}}V_{p_{n}}\otimes\frac{1}{2}(X-iY)+\sum_{q=1}^{l_{I_{n}}}v_{q_{n}}V_{q_{n}}^{\dagger}\otimes\frac{1}{2}(X+iY)\Biggr{)}$
(36)
$\displaystyle=\frac{1}{\sqrt{\tau}}\Biggr{(}\sum_{p=1}^{l_{I_{n}}}\frac{v_{p_{n}}}{2}\bigr{(}V_{p_{n}}\otimes
X+\hat{V}_{p_{n}}\otimes
Y\bigr{)}+\sum_{q=1}^{l_{I_{n}}}\frac{v_{p_{n}}}{2}\bigr{(}V_{q_{n}}^{\dagger}\otimes
X+\hat{V}_{q_{n}}^{\dagger}\otimes Y\bigr{)}\Biggr{)},$ (37)

where we have defined $\hat{V}=-iV$. Finally, we write the free bath
Hamiltonians, simply as

$H_{E_{n}}=\omega_{n}Z,$ (38)

which is the same as the counting operator $a^{\dagger}a$ up to an energy
shift. With these decompositions defined, we can discuss the construction of
the block encoding. This can be done more neatly by writing the total
Hamiltonian with a re-scaled interaction like so:

$A^{\tau}_{n}=H_{0}+H_{E_{n}}+H_{I_{n}}^{\tau}=\sum_{s=1}^{l_{0}+4l_{I_{n}}+1}\alpha_{s_{n}}A_{s_{n}}.$
(39)

Here, the indices work in the order that the respective Hamiltonian matrices
were defined above, i.e. $\\{A_{s_{n}}=V_{s_{n}}\otimes
X~{}:~{}s\in[l_{0}+1,l_{0}+l_{I_{n}}]\\}$. With this definition, the block
encoding machinery can be explained more clearly. Our construction is similar
to that of Ref. [5] We define the prepare operator $P$ as a unitary to load
the necessary coefficient. This works in the following way

$P_{n}\ket{0}^{\otimes
c}=\frac{1}{\sqrt{\alpha_{n}}}\sum_{\\{s_{n}\\}}\sqrt{\alpha_{s_{n}}}\ket{s},$
(40)

where $P_{n}$ is a unitary of dimension $c=l_{0}+4l_{I_{n}}+1$, and
$\alpha_{n}=\sum_{\\{s_{n}\\}}\alpha_{s_{n}}$ for normalization. We will refer
to this sum of coefficients as the alpha norm, which we be an important factor
in the final bound we present in the Section. Next, we build a select operator
$S_{n}$ which applies a certain Hamiltonian term conditioned on the
superposition state produced by $P_{n}$. The operator $S_{n}$ acts in the
following way:

$S_{n}\ket{s_{n}}\ket{\psi}=\ket{s_{n}}A_{s_{n}}\ket{\psi}.$ (41)

Here $\ket{s_{n}}$ represents our control state register, while $\ket{\psi}$
represents our input quantum state. In general, our input state is a mixed
quantum state $\rho$, but given the representation of the map in Equation 29,
this formulation is equivalent, and writing all of these operations as
channels makes the derivation unnecessarily complex. To remind the reader, the
$n$ subscript is allowing for the freedom to define different interactions
between the environment and system as well as the bath frequency $\omega_{n}$.
For this reason, a different choice of interaction may yield a different
unitary decomposition that requires a prepare operator $P$ that has a
different Hilbert space dimension. For this reason, we continue to subscript
our prepare and select unitaries.

With the above machinery defined, we now show how to utilize it to prepare a
block encoding:

$(P_{n}^{\dagger}\otimes\openone_{d})S(P_{n}\otimes\openone_{d})\ket{0}^{\otimes
c}\ket{\psi}=\ket{0}^{\otimes
c}\frac{A^{\tau}_{n}}{\alpha_{n}}\ket{\psi}+\ket{\phi^{\perp}},$ (42)

where $\ket{\phi^{\perp}}$ is a properly normalized state that is not of
interest, and $d=\dim\mathscr{H}_{0}$. With the above operation, we have
successfully produced the following block encoding of the matrix
$A_{n}^{\tau}$:

$(P_{n}^{\dagger}\otimes\openone_{d})S_{n}(P_{n}\otimes\openone_{d})=\scalebox{1.5}{
{\hbox{\begin{bmatrix}\frac{A_{n}^{\tau}}{\alpha_{n}} &\cdot\\\ \cdot&\cdot\\\
\end{bmatrix}}} },$ (43)

which can be accessed in the following way:

$(\bra{0}^{\otimes
c}\otimes\openone_{d})(P_{n}^{\dagger}\otimes\openone_{d})S_{n}(P_{n}\otimes\openone_{d})(\ket{0}^{\otimes
c}\otimes\openone_{d})=\frac{A_{n}^{\tau}}{\alpha_{n}}.$ (44)

The above example illustrates that the block encoding can be accessed upon
measurement of the $\ket{0}^{\otimes c}$. In implementations of LCU
strategies, such as in Ref. [5], it is required that one uses the amplitude
amplification algorithm to boost the probability of success $\delta$ of this
measurement. However, in our strategy we are interested in methods with
$\delta=1$. Therefore, we use Qubitization and QSVT to simulate our RI maps
which implement via queries to $P$ and $S$ defined above. The number of
queries $q_{PS}$ to the prepare and select oracles required to simulate
$\mathcal{L}$ with $\mathcal{R}$ up to precision $\epsilon$ is bounded in the
following Theorem.

######  Theorem 3 (Iterative Qubitization Complexity).

Given a RI map from Definition 1, an initial state $\rho_{S}$ and total time
$t$, the time evolution $\exp(\mathcal{L}t)$ of the corresponding Lindbladian
from Equation 15 can be implemented with the following number of queries $q$,
to both oracles $P$ from Equation 40 and $S$ from Equation 41:

$q_{PS}\in
O\Biggr{(}\max_{n}\biggr{(}t\nu(\alpha_{0}+m\omega_{n})+\alpha_{I_{n}}m\nu^{\frac{3}{2}}\sqrt{t}+\frac{m\nu\log\frac{1}{\epsilon}}{\log\log\frac{1}{\epsilon}}\biggr{)}\Biggr{)},$
(45)

with $\nu$ applications of an RI map composed as in Equation 14, and $m$
environment ancilla.

###### Proof.

The art of qubitization requires us to define a reflection operator about the
2-dimensional subspace in which the prepare state is embedded. This can be
achieved as follows:

$R_{n}=(\openone-2(P_{n}\otimes\openone_{d})\ket{0}\bra{0}(P_{n}^{\dagger}\otimes\openone_{d})).$
(46)

With this reflection, we construct the walk operator below, which is used in
the implementation of the dynamics of the $n$th subsystem

$W_{n}=R_{n}S_{n}.$ (47)

Qubitization bounds are formulated in terms of the number of queries to the
walk operator, however, for better interpretability, one can see that this
requires 1 query to select $S_{n}$ and 2 queries to prepare $P_{n}$, which is
equivalent in terms of asymptotics. The Quantum Singular Value Transformation
then provides the machinery for using this walk operator to implement the
operator exponential using a series of $X$ and $Z$ rotations in the
2-dimensional subspace. To see this discussion, consult Refs [25, 26]. The
aforementioned analysis leads to the following complexity upper-bound:

$q_{PS}\in O\Biggl{(}\alpha
t+\frac{\log\frac{1}{\epsilon}}{\log\log\frac{1}{\epsilon}}\Biggr{)},$ (48)

where $t$ is the total simulation time and $\alpha$ is the sum of the
coefficients from the unitary decomposition previously discussed. In the case
of the RI map, we can use this upper bound almost directly, with a few minor
adjustments. First, we are simulating each individual application of the map
with this algorithm a total of $\kappa$ times, meaning we make the replacement
$t\rightarrow\tau$, and we pick up a multiplicative factor of $\kappa$. Recall
from Corollary 2.1 that $\kappa=m\nu$ is the total number of repeated
interactions in the composite map defined in Equation 14. Also, given our
definition of $A_{N}^{\tau}$, the $\alpha$-norm argument goes through slightly
differently. For RI maps, the component of the $\alpha$-norm that comes from
the interaction Hamiltonian decomposition is necessarily a function of time,
given that the Hamiltonian has been re-scaled to $H_{I_{n}}^{\tau}$. It is
cleaner to instead leave it as a multiplicative factor, and we do the same for
the re-scaled free Hamiltonian $H_{0}$:

$\displaystyle\alpha_{n}$
$\displaystyle=\frac{1}{m}\alpha_{0}+\sqrt{\frac{1}{\tau}}\alpha_{I_{n}}+\alpha_{E_{n}}$
(49)
$\displaystyle=\sqrt{\frac{\nu}{t}}\alpha_{I_{n}}+(\frac{1}{m}\alpha_{0}+\omega_{n}).$
(50)

Now, applying Equation 48 and multiplying by a factor of $m$, and inserting
our defined $\alpha$-norms, we obtain the following bound:

$q_{PS}\in
O\Biggr{(}\max_{n}\biggr{(}t\nu(\alpha_{0}+m\omega_{n})+\alpha_{I_{n}}m\nu^{\frac{3}{2}}\sqrt{t}+\frac{m\nu\log\frac{1}{\epsilon}}{\log\log\frac{1}{\epsilon}}\biggr{)}\Biggr{)}.$
(51)

This completes the proof. ∎

Recall that Corollary 2.1 places an upper bound on $\nu$ to give the overall
time complexity.

###  VI.3 Trotter-Suzuki Formulas

In this Section we implement the RI map using the well known Trotter-Suzuki
product formulas [27, 28], and discuss some technicalities that arise as a
result of redefining the coupling $\lambda\rightarrow 1/\sqrt{\tau}$. Our
motivation for using these formulas are again that we can implement the
Lindbladian dynamics with $\delta=1$, for means of comparison with Ref. [21],
which first implemented a variant of these maps with Trotterization. For an
introduction to these product formulas, see Ref. [29], or the preliminaries
contained in Ref. [30]. Given that we have some exponentiated sum of operators
$e^{x\sum_{i}A_{i}}$ and we define the 1st and 2nd order product formulas as
follows:

$\displaystyle{S_{1}{\left(Ax\right)}}$
$\displaystyle=\prod_{i}^{L}e^{A_{i}x}$ (52)
$\displaystyle{S_{2}{\left(Ax\right)}}$
$\displaystyle=\prod_{i}^{L}e^{\frac{A_{i}x}{2}}\prod^{i}_{L}e^{\frac{A_{i}x}{2}}$
(53)

we can define a $2k$th-order product formula in the following way:

${S_{2k}{\left(Ax\right)}}:=({S_{2k-2}{\left(Axs_{2k}\right)}})^{2}{S_{2k-2}{\left((1-4s_{2k})Ax\right)}}({S_{2k-2}{\left(Axs_{2k}\right)}})^{2},$
(54)

where $s_{2k}=\frac{1}{4-4^{\frac{2k-1}{2}}}$. Then to generalize this to
density matrices we can define a Trotter channel in the following way:

######  Definition 4 (Trotter-Suzuki Channel).

Given a Hamiltonian $H$, a density matrix $\rho$, times $t$ and $t_{0}$
($t>t_{0}\geq 0$), and order 2k, then a Trotter-Suzuki channel
$\mathcal{T}^{2k}{(\rho(t_{0}),\Delta t)}$ performs the operation
$\mathcal{T}^{2k}{(\rho(t_{0}),\Delta t)}\rightarrow\rho(t)$ and can be
defined as:

$\mathcal{T}^{2k}{(\rho(t_{0}),t)}:={S_{2k}{\left(-iH\Delta
t\right)}}\rho{S_{2k}{\left(-iH\Delta t\right)}}^{\dagger}$ (55)

Where ${S_{2k}{\left(-iHt\right)}}$ represent the product formulae from
Equation (53) and $\Delta t=t-t_{0}$.

To control the error in $t$, it is common to implement the channel with for
time $t/r$ and iterate $r$ times. With these channels defined, we can now
bound the error of a RI Trotter implementation and make conclusions about the
cost of such a simulation. However, the input model in this case is different,
considering that we do not need to construct any block encodings, but rather
we implement the unitaries directly. To make this precise, given that we are
provided with a Hamiltonian $H_{n}$ at each interaction that is written as a
sum of operators $H_{n}=H_{0}+{H_{I_{n}}+H_{E_{n}}}=\sum_{j=1}^{l}H_{j}$ ,
with $l\in\text{poly}(\eta)$ where $\eta$ is the number of qubits. Now we
consider the cost of the algorithm as simply the number of exponential gates
needed to implement a product formula. This can likewise be thought of as
queries $q_{TS}(2k)$ to simple quantum circuits that provide the operator of
the Hamiltonian summands $e^{-iH_{j}t/r}$ with the desired time slice. In most
common applications of said decompositions, where we write our Hamiltonian as
a sum of Pauli operators or fermionic creation and annihilation operators,
simple circuits for such a task are well understood [31]. Note $2k$th order
product formulas contain $\Upsilon\equiv 2\cdot 5^{k-1}$ stages, where a stage
can be though of as a string that contains each of the $l$ Hamiltonian
summands $H_{j}$. Therefore, if we have time-slice $t/r$ and we implement the
entire product formula $r$ times, it is equivalent to carrying out
$q_{TS}(2k)=l\Upsilon r$ gates, where for a first order formula we simply have
1 stage. With this in mind, we derive the following bounds on $q_{TS}(2k)$ for
a RI map that implements Lindbladian dynamics:

######  Theorem 5 (Higher Order Trotter-Suzuki Cost).

Given a RI map from Definition 1, an initial state $\rho_{S}$ and total time
$t$, the time evolution $\exp(\mathcal{L}t)$ of the corresponding Lindbladian
from Equation 15 can be implemented with a $2k$th order Trotter-Suzuki channel
$\mathcal{T}^{2k}{(t)}$ from Definition 4, using the following number of
exponential gates $q_{TS}(2k)$:

$q_{TS}(2k)\in O\Biggr{(}l\Upsilon
m^{\frac{1}{2k}}\max_{n}\biggr{(}\frac{(\alpha_{0}/m+\omega_{n})^{1+\frac{1}{2k}}t^{1+\frac{1}{2k}}\nu^{\frac{1}{2k}}+\alpha_{I_{n}}t^{\frac{1}{2}+\frac{1}{4k}}\nu^{\frac{1}{2}+\frac{3}{4k}}}{\epsilon^{\frac{1}{2k}}}\biggr{)}\Biggr{)}.$
(56)

###### Proof.

The procedure consists of bounding the error we define by the induced one
norm, setting this error equal to $\epsilon$ and rearranging for $r$. Equipped
with an upper bound on $r$, we can then upper bound the total number of gates
using $q=l\Upsilon r$ as discussed above.

$\displaystyle\norm{\bigcirc_{n=1}^{m\nu}\mathcal{T}_{n}^{2k}(t/\nu)-e^{t\mathcal{L}}}_{1\rightarrow
1}$
$\displaystyle\leq\norm{\bigcirc_{n=1}^{m\nu}\mathcal{T}_{n}^{2k}(t/\nu)-\bigcirc_{n=1}^{m\nu}\widetilde{\mathcal{R}_{n}}(t)}_{1\rightarrow
1}+\norm{\bigcirc_{n=1}^{m\nu}\widetilde{\mathcal{R}_{n}}(t/\nu)-e^{t\mathcal{L}}}_{1\rightarrow
1}$ (57) $\displaystyle\leq
m\nu\max_{n}\norm{\mathcal{T}_{n}^{2k}(t)-\widetilde{\mathcal{R}_{n}}(t)}_{1\rightarrow
1}+\epsilon_{RI}$ (58)

In the second line, we utilize Lemma 8 from Ref. [30], a minor channel
composition inequality. We also observe that the induced 1-norm on the first
line is the RI error that has already been bounded previously. From here, we
can apply Theorem 2 from [32] to get the form of the Trotter error. The
aforementioned Theorem generalizes bounds from Ref. [28] to channels. In the
case of arbitrarily high order Trotter-Suzuki formulas, an analysis of the
commutator structure is highly non-trivial. In general, given the error
scaling of $m$, and the accuracy of the Lindblad equation itself, higher order
Trotter formulas may not be as useful as they are other applications. For
these reasons, we provide a trivial upper bound for the $2k$th order formulas
by implementing Lemma 1 from Ref. [28].

$\displaystyle\norm{\bigcirc_{n=1}^{m\nu}\mathcal{T}_{n}^{2k}(t)-e^{t\mathcal{L}}}_{1\rightarrow
1}$ $\displaystyle\in
O\Biggr{(}m\nu\max_{n}\biggr{(}\frac{\alpha_{n}^{2k+1}t^{2k+1}}{r^{2k}}\biggr{)}+\epsilon_{RI}\Biggr{)}$
(59) $\displaystyle\in
O\Biggr{(}m\max_{n}\biggr{(}\frac{(\alpha_{0}/m+\omega_{n})^{2k+1}t^{2k+1}\nu+\alpha_{I_{n}}t^{k+1/2}\nu^{k+3/2}}{r^{2k}}\biggr{)}+\epsilon_{RI}\Biggr{)}$
(60)

In both of the cases above, for the purposes of asymptotic analysis, we can
make the promise that the total simulation error
$\epsilon=\frac{1}{2}\epsilon_{RI}$, to avoid solving a non linear equation.
With this, we can simply rearrange the above bounds for $r$:

$\displaystyle r$ $\displaystyle\in
O\Biggr{(}m^{\frac{1}{2k}}\max_{n}\biggr{(}\frac{(\alpha_{0}/m+\omega_{n})^{1+\frac{1}{2k}}t^{1+\frac{1}{2k}}\nu^{\frac{1}{2k}}+\alpha_{I_{n}}t^{\frac{1}{2}+\frac{1}{4k}}\nu^{\frac{1}{2}+\frac{3}{4k}}}{\epsilon^{\frac{1}{2k}}}\biggr{)}\Biggr{)},$
(62)

and and multiplying by $l\Upsilon$ completes the proof. ∎

For simplicity, lets also consider the case of the first-order Trotter-Suzuki
channel. Give the more simple commutator sums that appear in the Trotter error
for low order formulas, we utilize some of the known commutator structure to
give more detailed bounds in the following Theorem.

######  Theorem 6 (Trotterization Cost).

Given a RI map from Definition 1, an initial state $\rho_{S}$ and total time
$t$, the time evolution $\exp(\mathcal{L}t)$ of the corresponding Lindbladian
from Equation 15 can be implemented with a 1st order Trotter-Suzuki channel
$\mathcal{T}^{2k}{(t)}$ from Definition 4, using the following number of
exponential gates $q_{TS}(1)$:

$q_{TS}(1)\leq\
\Biggr{(}\max_{n}\frac{\sqrt{\nu}t^{3/2}}{2\epsilon}\Bigr{(}\sum_{j=1}^{l_{0}}\norm{[H_{0_{j}},H_{I_{n}}]}+m\norm{[H_{E_{n}},H_{I_{n}}]}\Bigr{)}+\frac{\nu
t^{2}}{2m\epsilon}\sum_{p=1,q=1}^{l_{0},l_{0}}\norm{[H_{0_{p}},H_{0_{q}}]}\Biggr{)}$
(63)

###### Proof.

Starting from Equation 58, we can once again apply Theorem 2 from [32] to
insert the form of the Trotter error. Note that given the low order of the
formula, we can more simply avail some of the commutator structure in the RI
map. We have in general $[H_{0},H_{I_{n}}^{\tau}]\neq 0$,
$[H_{E_{n}},H_{I_{n}}^{\tau}]\neq 0$, and $[H_{E_{n}},H_{0}]=0$.

$\displaystyle\norm{\bigcirc_{n=1}^{m\nu}\mathcal{T}_{n}^{1}(t)-e^{t\mathcal{L}}}_{1\rightarrow
1}\leq
m\max_{n}\norm{\mathcal{T}_{n}^{2k}(t)-\widetilde{\mathcal{R}_{n}}(t)}_{1\rightarrow
1}+\epsilon_{RI}$ (64) $\displaystyle\leq
m\nu\max_{n}\frac{t^{2}}{r}\Bigr{(}\sum_{j=1}^{l_{0}}\frac{1}{m}\norm{[H_{0_{j}},H_{I_{n}}^{\tau}]}+\norm{[H_{E_{n}},H_{I_{n}}^{\tau}]}+\sum_{p=1,q=1}^{l_{0},l_{0}}\frac{1}{m^{2}}\norm{[H_{0_{p}},H_{0_{q}}]}\Bigr{)}+\epsilon_{RI}$
(65)
$\displaystyle=m\sqrt{\nu}\max_{n}\frac{t^{3/2}}{2r}\Bigr{(}\sum_{j=1}^{l_{0}}\frac{1}{m}\norm{[H_{0_{j}},H_{I_{n}}]}+\norm{[H_{E_{n}},H_{I_{n}}]}\Bigr{)}+\frac{\nu
mt^{2}}{2r}\sum_{p=1,q=1}^{l_{0},l_{0}}\frac{1}{m^{2}}\norm{[H_{0_{p}},H_{0_{q}}]}+\epsilon_{RI}$
(66)

where in the previous line we pulled the factor of $1/\sqrt{\tau}$ out of the
interaction Hamiltonian. Once again making the promise that
$\epsilon=\frac{1}{2}\epsilon_{RI}$, and rearranging for $r$ we obtain the
following:

$r\leq\sqrt{\nu}\max_{n}\frac{t^{3/2}}{2\epsilon}\Bigr{(}\sum_{j=1}^{l_{0}}m\norm{[H_{0_{j}},H_{I_{n}}]}+\norm{[H_{E_{n}},H_{I_{n}}]}\Bigr{)}+\frac{m\nu
t^{2}}{2m\epsilon}\sum_{p=1,q=1}^{l_{0},l_{0}}\norm{[H_{0_{p}},H_{0_{q}}]},$
(67)

now we can trivially bound the number of operator exponentials $q_{TS}(1)$ by
multiplying by the total number of terms in the Hamiltonian $l$ which
completes the proof. ∎

As in the case of the iterative Qubitization simulation, Corollary 2.1 can be
applied to upper bound the number of applications of the composed RI map
defined in Equation 14.

##  VII Discussion

In this work, we successfully bounded the error between a Lindbladian
evolution and a Repeated Interaction map in Theorem 2. In doing so, we showed
the the error in this approximate correspondence can be made arbitrarily small
by increasing the number of interactions between the quantum system of
interest and environmental subsystems. This also implies that dissipative
Lindbladian dynamics can be efficiently embedded as a sequence of unitary
Hamiltonian simulations with trace out. As a consequence, this analysis also
places an upper bound on the number of ancilla qubits needed, given that we
perform a deferred trace out to maintain unitarity for the entire simulation
time $t$. We then implemented iterated Qubitization and QSVT strategies to
yield an asymptotically optimal implementations of RI maps, as well as 2
Trotter-Suzuki simulations for comparison, and complexity upper bounds were
given in Theorems 3, 5, and 6.

What is somewhat unique to this problem, is that even in generality, some
commutator structure is known a priori. This potentially makes for a good
candidate for composite simulation strategies such as those discussed in Refs.
[32, 30]. Specifically, in Ref. [32] it was proven that if one can find a
partition $H=A+B$ such that the number of non-zero commutators in $A$ scales
like the square root of the number of terms in the set $|A|$, and the terms in
$B$ are sufficiently small, then it is always possible to gain an asymptotic
advantage by simulating $A$ with Trotter and $B$ with QDrift. Whether the RI
map guarantees the construction of the aforementioned sets is an open question
(see Section 4.3 of [32]).

This problem is also interesting due to the two errors that are coming in, one
from the approximate correspondence and the other from the simulation.
Interestingly, the asymptotics suggest that implementing high orders of the RI
map are not guaranteed to reduce $\epsilon_{RI}$. This means that iteratively
implementing just the first two orders of the map, in a non-unitary fashion,
may lead to a more efficient simulation. For instance, constructing a Kraus
map that approximates the first two orders of $\widetilde{\mathcal{R}_{n}}$
and then applying an LCU strategy, similar to the approach of [12], may
outperform the prior approach when it comes to constant factors, given that
the constructions of $P$ and $S$ will be far more simple. However, our goal in
this work was to provide unitary algorithms based on Hamiltonian simulation
with $\delta=1$, so we leave this approach as a future avenue to be explored.

From a broader perspective, our work reveals an issue with current approaches
to simulating open system dynamics. Often the master equations that we aim to
simulate have fundamental errors that make the dynamics only approximately
equivalent to the exact unitary dynamics that describes a system bath
coupling. Understanding these errors in the context of open systems simulation
is important since bounding such discrepancies may be vital to fairly compare
different open systems simulation schemes. In particular, a complete closed
system simulation for a system may only be an approximation to Lindbladian
dynamics but it may be a more accurate description of the quantum dynamics
than the Lindblad master equation because it does not make the Markov
approximation. While our work makes an important step in this direction, we
feel more work is needed in order to understand how best to use approaches
from open quantum systems to simulate quantum dynamics with provable error
guarantees, while also retaining the conceptual simplicity and clarity that
the open systems approach provides.

##  VIII Acknowledgements

We thank Matthew Hagan, Danial Motlagh, and Marlon Brenes for valuable
discussions about this work as well as Lin Lin and Xiantao Li for useful
feedback on our work. MP and NW acknowledge funding from the NSERC Discovery
Program. NW’s work was further supported by the US Department of Energy,
Office of Science, National Quantum Information Science Research Centers, Co-
Design Center for Quantum Advantage under contract number DE-SC0012704. D.S.
acknowledges support from an NSERC Discovery Grant and the Canada Research
Chair program.

## References

  * Feynman  _et al._ [2018] R. P. Feynman _et al._ , Simulating physics with computers, Int. j. Theor. phys 21 (2018).
  * Lloyd [1996] S. Lloyd, Universal quantum simulators, Science 273, 1073 (1996).
  * Low and Chuang [2019] G. H. Low and I. L. Chuang, Hamiltonian simulation by qubitization, Quantum 3, 163 (2019).
  * Childs and Wiebe [2012] A. M. Childs and N. Wiebe, Hamiltonian simulation using linear combinations of unitary operations, arXiv preprint arXiv:1202.5822 (2012).
  * Berry  _et al._ [2015] D. W. Berry, A. M. Childs, R. Cleve, R. Kothari, and R. D. Somma, Simulating hamiltonian dynamics with a truncated taylor series, Physical Review Letters 114, 10.1103/physrevlett.114.090502 (2015).
  * Babbush  _et al._ [2023] R. Babbush, D. W. Berry, R. Kothari, R. D. Somma, and N. Wiebe, Exponential quantum speedup in simulating coupled classical oscillators, arXiv preprint arXiv:2303.13012 (2023).
  * Berkelbach and Thoss [2020] T. C. Berkelbach and M. Thoss, Special topic on dynamics of open quantum systems, The Journal of Chemical Physics 152, 020401 (2020), https://pubs.aip.org/aip/jcp/article-pdf/doi/10.1063/1.5142731/15568495/020401_1_online.pdf .
  * Breuer and Petruccione [2002] H.-P. Breuer and F. Petruccione,  _The Theory of Open Quantum Systems_ (Oxford University Press, 2002).
  * Rivas and Huelga [2012] A. Rivas and S. F. Huelga,  _Open quantum systems_ , Vol. 10 (Springer, 2012).
  * Ciccarello  _et al._ [2022] F. Ciccarello, S. Lorenzo, V. Giovannetti, and G. M. Palma, Quantum collision models: Open system dynamics from repeated interactions, Physics Reports 954, 1 (2022).
  * Bruneau  _et al._ [2014] L. Bruneau, A. Joye, and M. Merkli, Repeated interactions in open quantum systems, Journal of Mathematical Physics 55 (2014).
  * Cleve and Wang [2016] R. Cleve and C. Wang, Efficient quantum algorithms for simulating lindblad evolution, arXiv preprint arXiv:1612.09512 (2016).
  * Li and Wang [2022] X. Li and C. Wang, Simulating markovian open quantum systems using higher-order series expansion, arXiv preprint arXiv:2212.02051 (2022).
  * 201 [2017] Quantum Information and Computation 17, 10.26421/qic17.11-12 (2017).
  * Barthel and Kliesch [2012] T. Barthel and M. Kliesch, Quasilocality and efficient simulation of markovian quantum dynamics, Physical review letters 108, 230504 (2012).
  * Kliesch  _et al._ [2011] M. Kliesch, T. Barthel, C. Gogolin, M. Kastoryano, and J. Eisert, Dissipative quantum church-turing theorem, Physical review letters 107, 120501 (2011).
  * Schlimgen  _et al._ [2021] A. W. Schlimgen, K. Head-Marsden, L. M. Sager, P. Narang, and D. A. Mazziotti, Quantum simulation of open quantum systems using a unitary decomposition of operators, Physical Review Letters 127, 270503 (2021).
  * Schlimgen  _et al._ [2022] A. W. Schlimgen, K. Head-Marsden, L. M. Sager, P. Narang, and D. A. Mazziotti, Quantum simulation of the lindblad equation using a unitary decomposition of operators, Physical Review Research 4, 023216 (2022).
  * Ding  _et al._ [2023] Z. Ding, X. Li, and L. Lin, Simulating open quantum systems using hamiltonian simulations, arXiv preprint arXiv:2311.15533 (2023).
  * Attal and Pautrat [2006] S. Attal and Y. Pautrat, From repeated to continuous quantum interactions, in  _Annales Henri Poincaré_ , Vol. 7 (Springer, 2006) pp. 59–104.
  * Cattaneo  _et al._ [2021] M. Cattaneo, G. De Chiara, S. Maniscalco, R. Zambrini, and G. L. Giorgi, Collision models can efficiently simulate any multipartite markovian quantum dynamics, Physical Review Letters 126, 130403 (2021).
  * Pellegrini and Petruccione [2009] C. Pellegrini and F. Petruccione, Non-markovian quantum repeated interactions and measurements, Journal of Physics A: Mathematical and Theoretical 42, 425304 (2009).
  * Nielsen and Chuang [2010] M. A. Nielsen and I. L. Chuang,  _Quantum computation and quantum information_ (Cambridge university press, 2010).
  * Kliuchnikov [2013] V. Kliuchnikov, Synthesis of unitaries with clifford+ t circuits, arXiv preprint arXiv:1306.3200 (2013).
  * Gilyén  _et al._ [2019] A. Gilyén, Y. Su, G. H. Low, and N. Wiebe, Quantum singular value transformation and beyond: exponential improvements for quantum matrix arithmetics, in  _Proceedings of the 51st Annual ACM SIGACT Symposium on Theory of Computing_ (2019) pp. 193–204.
  * Martyn  _et al._ [2021] J. M. Martyn, Z. M. Rossi, A. K. Tan, and I. L. Chuang, Grand unification of quantum algorithms, PRX Quantum 2, 040203 (2021).
  * Suzuki [1990] M. Suzuki, Fractal decomposition of exponential operators with applications to many-body theories and monte carlo simulations, Physics Letters A 146, 319 (1990).
  * Childs  _et al._ [2021] A. M. Childs, Y. Su, M. C. Tran, N. Wiebe, and S. Zhu, Theory of trotter error with commutator scaling, Physical Review X 11, 011020 (2021).
  * Hatano and Suzuki [2005] N. Hatano and M. Suzuki, Finding exponential product formulas of higher orders, in  _Quantum Annealing and Other Optimization Methods_ (Springer Berlin Heidelberg, 2005) pp. 37–68.
  * Pocrnic  _et al._ [2023] M. Pocrnic, M. Hagan, J. Carrasquilla, D. Segal, and N. Wiebe, Composite qdrift-product formulas for quantum and classical simulations in real and imaginary time, arXiv preprint arXiv:2306.16572 (2023).
  * Whitfield  _et al._ [2011] J. D. Whitfield, J. Biamonte, and A. Aspuru-Guzik, Simulation of electronic structure hamiltonians using quantum computers, Molecular Physics 109, 735 (2011).
  * Hagan and Wiebe [2022] M. Hagan and N. Wiebe, Composite quantum simulations, arXiv preprint arXiv:2206.06409 (2022).
  * Watrous [2018] J. Watrous,  _The theory of quantum information_ (Cambridge university press, 2018).

##  Appendix A General Error Bounding

This Section is dedicated towards bounding the following induced 1-norm
$\norm{e^{\tau\mathcal{L}}-\mathcal{R}_{n}(\tau)}_{1\rightarrow 1}$. In order
to do this, and for simplicity, we assume that the Lindblad generator is
exactly equivalent to the repeated interaction Lindblad form derived in
Section IV. This will allow for the cancellation of some terms in the series
expansion. Rather than yielding the best overall bound in this work, this
section is written to stress the importance of the limiting case of the map.
Note that we are here working with the general RI map $\mathcal{R}_{n}(\tau)$
rather than the limiting case $\widetilde{\mathcal{R}_{n}}(\tau)$ used in the
Lindbladian derivation in the main text. The overall error expression between
the Lindblad dynamics and the General RI is give by

$\displaystyle\norm{e^{\tau\mathcal{L}}-\mathcal{R}_{n}(\tau)}_{1\rightarrow
1}\leq\underbrace{\norm{e^{\tau\mathcal{L}}-\widetilde{\mathcal{R}_{n}}(\tau)}_{1\rightarrow
1}}_{\epsilon_{\tau}}+\underbrace{{\norm{\widetilde{\mathcal{R}_{n}}(\tau)-\mathcal{R}_{n}(\tau)}_{1\rightarrow
1}}}_{\epsilon_{\lambda}}.$ (68)

Here $\widetilde{\mathcal{R}_{n}}(\tau)$ is the repeated interaction map in
the $\lambda^{2}\tau\rightarrow 1$ limit, or, when expanded, the map that
leads to Equation (IV). Here, we split the error into two terms,
$\epsilon_{\tau}$ for the error that arises from the fact that $\tau$ does not
precisely go to zero in the repeated interactions map, and the other
$\epsilon_{\lambda}$ which is the difference arising from $\lambda$ not going
to $\sqrt{\nicefrac{{1}}{{\tau}}}$ in general. Note that if we are dealing
with an RI map composed of the full $m$ interactions, then
$\epsilon_{\tau}=\epsilon_{RI}$, or the error bounded in Theorem 2. Therefore,
we are calculating the error between a Lindblad evolution and a general
repeated interactions map with no additional assumptions:

$\displaystyle\epsilon_{\tau}$
$\displaystyle=\norm{e^{\tau\mathcal{L}}-\widetilde{\mathcal{R}_{n}}(\tau)}_{1\rightarrow
1}$ (69)
$\displaystyle=\norm{\int_{0}^{\tau}\partial_{x}\left(e^{x\mathcal{L}}-\widetilde{\mathcal{R}_{n}}(x)\right)dx+\left(e^{0\mathcal{L}}-\widetilde{\mathcal{R}_{n}}(0)\right)}_{1\rightarrow
1}$ (70)
$\displaystyle=\norm{\int_{0}^{\tau}\partial_{x}\left(e^{x\mathcal{L}}-\widetilde{\mathcal{R}_{n}}(x)\right)dx}_{1\rightarrow
1}$ (71)
$\displaystyle\leq\int_{0}^{\tau}\norm{\partial_{x}\left(e^{x\mathcal{L}}-\widetilde{\mathcal{R}_{n}}(x)\right)}_{1\rightarrow
1}dx$ (72)
$\displaystyle\leq\tau\max_{x\in[0,\tau]}\Biggr{|}\Biggr{|}\mathcal{L}\sum_{n=1}^{\infty}\frac{\mathcal{L}^{n}x^{n}}{n!}-\Tr_{\mathscr{E}_{n}}\Biggl{(}\partial_{x}\sum_{m=3}^{\infty}\frac{\widetilde{\mathcal{H}}^{m}x^{m}}{m!}\Biggr{)}\Biggr{|}\Biggr{|}_{1\rightarrow
1}$
$\displaystyle\quad+\max_{\rho:\norm{\rho}_{1}=1}\Biggr{|}\Biggr{|}\Tr_{\mathscr{E}_{n}}\Biggr{(}\frac{3}{4}\tau^{1/2}\Bigl{(}[(H_{0}+H_{E_{n}}),[H_{I_{n}},\rho_{S}\otimes\rho_{E_{n}}]][H_{I_{n}},[(H_{0}+H_{E_{n}}),\rho_{S}\otimes\rho_{E_{n}}]]\Bigr{)}+$
$\displaystyle\quad+\tau[(H_{0}+H_{E_{n}}),[(H_{0}+H_{E_{n}}),\rho_{S}\otimes\rho_{E_{n}}]]\Biggr{)}\Biggr{|}\Biggr{|}$
(73)

Looking at the above result, we simplify it by using the fact that
$\Tr_{E}[H_{E},\rho_{S}\otimes\rho_{E_{n}}]=0$ and (based on the assumption in
the derivation) $\Tr_{E}[H_{I},\rho_{S}\otimes\rho_{E_{n}}]=0$, as well as the
fact that $H_{0}$ commutes with the partial trace, we have that

$\Tr_{\mathscr{E}_{n}}\frac{3}{4}\tau^{1/2}\Bigl{(}[(H_{0}+H_{E_{n}}),[H_{I_{n}},\rho_{S}\otimes\rho_{E_{n}}]]+[H_{I_{n}},[(H_{0}+H_{E_{n}}),\rho_{S}\otimes\rho_{E_{n}}]]\Bigr{)}=0.$
(74)

After applying the triangle inequality and making some smallness assumptions,
we can see that, to leading order we will have:

$\epsilon_{\tau}\in
O(\tau^{2}\norm{[(H_{0}+H_{E_{n}}),[(H_{0}+H_{E_{n}}),\rho_{S}\otimes\rho_{E_{n}}]]}_{1\rightarrow
1}).$ (75)

This term contributes an error that scales like $\tau^{2}$, which without
additional assumptions seems to be the best one can approximate Lindbladian
evolution with an RI map. We will make these arguments more rigorous in the
next section. However, this is not the case with the other remaining error
term $\epsilon_{\lambda}$. Turning our focus to $\epsilon_{\lambda}$, and
simply expanding each map and cancelling like terms we have:

$\displaystyle\epsilon_{\lambda}$
$\displaystyle=\norm{\widetilde{\mathcal{R}_{n}}(\tau)-\mathcal{R}_{n}(\tau)}_{1\rightarrow
1}$ (76)
$\displaystyle=\max_{\rho:\norm{\rho}_{1}=1}\Biggr{|}\Biggr{|}\Tr_{\mathscr{E}_{n}}\Biggl{(}\sum_{n=3}^{\infty}\frac{\mathcal{H}^{n}x^{n}}{n!}-\sum_{m=3}^{\infty}\frac{\widetilde{\mathcal{H}}^{m}x^{m}}{m!}+\frac{1}{2}\Bigl{(}(\tau^{2}\lambda-\tau^{3/2})([H_{I_{n}},[(H_{0}+H_{E_{n}}),\rho_{S}\otimes\rho_{E_{n}}]]$
$\displaystyle\quad+[(H_{0}+H_{E_{n}}),[H_{I_{n}},\rho_{S}\otimes\rho_{E_{n}}]])+(\tau^{2}\lambda^{2}-\tau)[H_{I_{n}},[H_{I_{n}},\rho_{S}\otimes\rho_{E_{n}}]]\Bigr{)}\Biggr{)}\Biggr{|}\Biggr{|}_{1}$
(77)

We can proceed with repeated applications of the triangle inequality and sub-
multiplicative property. We will also utilize the von Neumann trace inequality
to bound the partial traces. Furthermore, we utilize the property

$\norm{\mathcal{H}}_{1\rightarrow 1}\leq 2\norm{H}\leq
2(\norm{H_{0}}+\norm{H_{E_{n}}}+\lambda\norm{H_{I_{n}}})$ (78)

This bound is saturated when the density matrix support is constrained to the
1-D subspace that coincides with the largest eigenvalue of the Hamiltonian.
Proceeding with applications of the triangle inequality we have:

$\displaystyle\epsilon_{\lambda}$
$\displaystyle\leq\norm{\Tr_{\mathscr{E}_{n}}\sum_{n=3}^{\infty}\frac{\mathcal{H}^{n}x^{n}}{n!}}_{1\rightarrow
1}+\norm{\Tr_{\mathscr{E}_{n}}\sum_{m=3}^{\infty}\frac{\widetilde{\mathcal{H}}^{m}x^{m}}{m!}}_{1\rightarrow
1}$
$\displaystyle\quad+\frac{1}{2}\Bigl{(}(\tau^{2}\lambda+\tau^{3/2})(\norm{\Tr_{\mathscr{E}_{n}}[H_{I_{n}},[(H_{0}+H_{E_{n}}),\rho_{S}\otimes\rho_{E_{n}}]]}_{1\rightarrow
1}$
$\displaystyle\quad+\norm{\Tr_{\mathscr{E}_{n}}[(H_{0}+H_{E_{n}}),[H_{I_{n}},\rho_{S}\otimes\rho_{E_{n}}]]}_{1\rightarrow
1})+(\tau^{2}\lambda^{2}+\tau)\Tr_{\mathscr{E}_{n}}\norm{[H_{I_{n}},[H_{I_{n}},\rho_{S}\otimes\rho_{E_{n}}]]}_{1\rightarrow
1}\Bigr{)},$ (79)

Immediately we notice something problematic with this error. The following
term,

$(\tau^{2}\lambda^{2}+\tau)\Tr_{\mathscr{E}_{n}}\norm{[H_{I_{n}},[H_{I_{n}},\rho_{S}\otimes\rho_{E_{n}}]]}_{1\rightarrow
1}\neq 0,\>\forall\tau>0.$ (80)

as we know the double commutator with the interaction Hamiltonian to
necessarily be non-zero [11]. Further bounding this term, we make the
assumption on the form of the interaction Hamiltonian $H_{I_{n}}=S\otimes
B_{n}$,

$\displaystyle\Tr_{\mathscr{E}_{n}}\norm{[H_{I_{n}},[H_{I_{n}},\rho_{S}\otimes\rho_{E_{n}}]]}_{1\rightarrow
1}$
$\displaystyle\leq\absolutevalue{\Tr_{\mathscr{E}_{n}}(B_{n}^{2}\rho_{E_{n}})}\norm{S^{2}\rho_{S}+\rho_{S}S^{2}-2S\rho_{S}S}_{1\rightarrow
1}$ (81)
$\displaystyle\leq\sum_{j}e_{j}(B^{2})e_{j}(\rho_{E_{n}})\norm{[H_{I_{n}},[H_{I_{n}},\rho_{S}\otimes\rho_{E_{n}}]]}_{1\rightarrow
1}$ (82)
$\displaystyle\leq\norm{B^{2}}_{\infty}\norm{S^{2}\rho_{S}+\rho_{S}S^{2}-2S\rho_{S}S}_{1\rightarrow
1}$ (83)
$\displaystyle\leq\norm{B^{2}}_{\infty}(2\norm{\rho_{S}}_{1}\norm{S^{2}}_{\infty}+2\norm{\rho_{S}}_{1}\norm{S}_{\infty}^{2})$
(84) $\displaystyle\leq 4\norm{B^{2}}_{\infty}\norm{S}_{\infty}^{2}$ (85)

where $e_{j}(B)$ represent the ordered eigenvalues of the matrix $B$. We the
von Neumann trace inequality in the second line, and use identity from Ref.
[33] P.33 (Eq. 1.175) in the third line. We could proceed to bound the
remainder of the terms outside the sums in the same fashion, however, the
takeaway here is that without using the limiting RI map
$\widetilde{\mathcal{R}_{n}}(\tau)$, we obtain an error in the first non-
trivial order in the interaction time $\tau$. To leading order, again with
smallness assumptions, we have:

$\epsilon_{\lambda}\in O(\tau\norm{B^{2}}_{\infty}\norm{S}_{\infty}^{2}).$
(86)

Therefore, implementing an RI map without the limiting case of
$\lambda^{2}\tau\rightarrow 1$, would necessarily impose severe restrictions
on the interaction time $\tau$, and consequentially on number of required
interactions $m$ which would result in poor complexity bounds. For this
reason, we focus only RI maps where our desired limit or replacement of
variables holds.