$\displaystyle\leq\norm{\frac{e^{-\nicefrac{{\beta}}{{N}}\mathcal{L}}}{1+\Tr(\sum_{n=1}^{\infty}\frac{(-\beta\mathcal{L})^{n}}{N^{n}n!})}-\frac{\sum_{i}\frac{h_{i}}{\lambda}e^{-\nicefrac{{\beta\lambda}}{{N}}\mathcal{L}_{i}}}{1+\Tr(\sum_{n=1}^{\infty}\frac{(-\beta\mathcal{L})^{n}}{N^{n}n!})}}_{1\rightarrow
1}$
$\displaystyle+\norm{\frac{\sum_{i}\frac{h_{i}}{\lambda}e^{-\nicefrac{{\beta\lambda}}{{N}}\mathcal{L}_{i}}}{1+\Tr(\sum_{n=1}^{\infty}\frac{(-\beta\mathcal{L})^{n}}{N^{n}n!})}-\frac{\sum_{i}\frac{h_{i}}{\lambda}e^{-\nicefrac{{\beta\lambda}}{{N}}\mathcal{L}_{i}}}{1+\Tr(\sum_{j}\frac{h_{j}}{\lambda}{\sum_{m=1}^{\infty}\frac{(-\beta\mathcal{L}_{j})^{m}}{N^{m}m!}})}}_{1\rightarrow
1}$ (60) $\displaystyle\coloneqq a+b$ (61)

Now we shall proceed by dealing with each of these terms individually. The
denominator of the $a$ can be expanded via its geometric series:

$\displaystyle a$
$\displaystyle=\norm{\left(\sum_{k=0}^{\infty}\bigl{(}-\Tr(\sum_{n=1}^{\infty}\frac{(-\beta\mathcal{L})^{n}}{N^{n}n!})\bigr{)}^{k}\right)(e^{\nicefrac{{-\beta\mathcal{L}}}{{N}}}-\sum_{i}\frac{h_{i}}{\lambda}e^{\nicefrac{{-\beta\mathcal{L}_{i}}}{{N}}})}_{1\rightarrow
1}$ (62)
$\displaystyle=\absolutevalue{\sum_{k=0}^{\infty}(-\Tr\bigl(\sum_{n=1}^{\infty}\frac{(-\beta\mathcal{L})^{n}}{N^{n}n!})\bigr{)}^{k}}\norm{e^{\nicefrac{{-\beta\mathcal{L}}}{{N}}}-\sum_{i}\frac{h_{i}}{\lambda}e^{\nicefrac{{-\beta\lambda\mathcal{L}_{i}}}{{N}}}}_{1\rightarrow
1}$ (63)

We will first focus on the leftmost factor. To bound this term we will apply
the triangle inequality to apply the norm to each term in the sum. We will
also use the fact that the absolute value of the trace is trivially upper-
bounded by the induced 1-norm:
$\absolutevalue{\Tr\mathcal{L}(\rho)}\leq\norm{\mathcal{L}(\rho)}_{1\rightarrow
1}$. Here it is understood that this absolute value is maximized over all
$\rho:\norm{\rho}=1$. We then use the sub-multiplicative property
$\norm{\mathcal{L}^{n}}_{1\rightarrow 1}\leq\norm{\mathcal{L}}_{1\rightarrow
1}^{n}$:

$\displaystyle a$
$\displaystyle\leq\sum_{k=0}^{\infty}\bigl{(}\sum_{n=1}^{\infty}\frac{\beta^{n}\absolutevalue{\Tr(\mathcal{L})^{n}}}{N^{n}n!}\bigr{)}^{k}\norm{e^{\nicefrac{{-\beta\mathcal{L}}}{{N}}}-\sum_{i}\frac{h_{i}}{\lambda}e^{\nicefrac{{-\beta\lambda\mathcal{L}_{i}}}{{N}}}}_{1\rightarrow
1}$ (65)
$\displaystyle\leq\sum_{k=0}^{\infty}\bigl{(}\sum_{n=1}^{\infty}\frac{\beta^{n}\norm{\mathcal{L}}_{1\rightarrow
1}^{n}}{N^{n}n!}\bigr{)}^{k}\norm{e^{\nicefrac{{-\beta\mathcal{L}}}{{N}}}-\sum_{i}\frac{h_{i}}{\lambda}e^{\nicefrac{{-\beta\lambda\mathcal{L}_{i}}}{{N}}}}_{1\rightarrow
1}$ (66)
$\displaystyle\leq\sum_{k=0}^{\infty}\bigl{(}\sum_{n=1}^{\infty}\frac{\beta^{n}2^{n}\lambda^{n}}{N^{n}n!}\bigr{)}^{k}\norm{e^{\nicefrac{{-\beta\mathcal{L}}}{{N}}}-\sum_{i}\frac{h_{i}}{\lambda}e^{\nicefrac{{-\beta\lambda\mathcal{L}_{i}}}{{N}}}}_{1\rightarrow
1}$ (67)
$\displaystyle=\sum_{k=0}^{\infty}(e^{\nicefrac{{2\lambda\beta}}{{N}}}-1)^{k}\norm{e^{\nicefrac{{-\beta\mathcal{L}}}{{N}}}-\sum_{i}\frac{h_{i}}{\lambda}e^{\nicefrac{{-\beta\lambda\mathcal{L}_{i}}}{{N}}}}_{1\rightarrow
1}$ (68)
$\displaystyle=\frac{\norm{e^{\nicefrac{{-\beta\mathcal{L}}}{{N}}}-\sum_{i}\frac{h_{i}}{\lambda}e^{\nicefrac{{-\beta\lambda\mathcal{L}_{i}}}{{N}}}}_{1\rightarrow
1}}{2-e^{\nicefrac{{2\beta\lambda}}{{N}}}}$ (69)

Where we require $|e^{\nicefrac{{2\beta\lambda}}{{N}}}-1|<1$. Next, inspecting
the numerator, this is of the form bounded in the Appendix of [11], therefore,
the following holds trivially:

$a\leq\frac{4\beta^{2}\lambda^{2}}{N^{2}}\frac{e^{\nicefrac{{2\lambda\beta}}{{N}}}}{2-e^{\nicefrac{{2\beta\lambda}}{{N}}}}$
(70)

Next, let us proceed to bound the $b$ term:

$\displaystyle b$
$\displaystyle=\norm{\frac{\sum_{i}\frac{h_{i}}{\lambda}e^{-\beta\mathcal{L}_{i}}}{1+\Tr(\sum_{n=1}^{\infty}\frac{(-\beta\mathcal{L})^{n}}{n!})}-\frac{\sum_{i}\frac{h_{i}}{\lambda}e^{-\beta\mathcal{L}_{i}}}{1+\Tr(\sum_{j}\frac{h_{j}}{\lambda}{\sum_{m=1}^{\infty}\frac{(-\beta\mathcal{L}_{j})^{m}}{m!}})}}_{1\rightarrow
1}$ (71)
$\displaystyle=\norm{\sum_{i}\frac{h_{i}}{\lambda}e^{-\beta\mathcal{L}_{i}}}_{1\rightarrow
1}\absolutevalue{\sum_{k=0}^{\infty}\bigl{(}-\Tr(\sum_{n=1}^{\infty}\frac{(-\beta\mathcal{L})^{n}}{n!})\bigr{)}^{k}-\sum_{l=0}^{\infty}\bigl{(}-\Tr(\sum_{j}\frac{h_{j}}{\lambda}\sum_{n=1}^{\infty}\frac{(-\beta\mathcal{L}_{j})^{n}}{n!})\bigr{)}^{l}}$
(72) $\displaystyle\leq
e^{\nicefrac{{2\beta\lambda}}{{N}}}\absolutevalue{\sum_{k=0}^{\infty}\bigl{(}-\Tr(\sum_{n=1}^{\infty}\frac{(-\beta\mathcal{L})^{n}}{n!})\bigr{)}^{k}-\sum_{l=0}^{\infty}\bigl{(}-\Tr(\sum_{j}\frac{h_{j}}{\lambda}\sum_{n=1}^{\infty}\frac{(-\beta\mathcal{L}_{j})^{n}}{n!})\bigr{)}^{l}}$
(73)

Now, carefully expanding the double sums up to second order:

$=e^{\nicefrac{{2\beta\lambda}}{{N}}}\absolutevalue{\begin{multlined}\sum_{k=2}^{\infty}\bigl{(}-\Tr(\sum_{n=2}^{\infty}\frac{(-\beta\mathcal{L})^{n}}{n!})\bigr{)}^{k}-\sum_{l=2}^{\infty}\bigl{(}-\Tr(\sum_{j}\frac{h_{j}}{\lambda}\sum_{n=2}^{\infty}\frac{(-\beta\mathcal{L}_{j})^{n}}{n!})\bigr{)}^{l}+\\\
\Tr\sum_{j}\frac{h_{j}}{\lambda}\sum_{m=2}^{\infty}\frac{(-\beta\lambda\mathcal{L}_{j})^{m}}{N^{m}m!}-\Tr\sum_{n=2}^{\infty}\frac{(-\beta\lambda\mathcal{L})^{n}}{N^{n}m!}+\Tr\frac{\beta\mathcal{L}}{N}-\Tr\sum_{j}\frac{h_{j}}{\lambda}\frac{\beta\lambda\mathcal{L}_{j}}{N}\end{multlined}\sum_{k=2}^{\infty}\bigl{(}-\Tr(\sum_{n=2}^{\infty}\frac{(-\beta\mathcal{L})^{n}}{n!})\bigr{)}^{k}-\sum_{l=2}^{\infty}\bigl{(}-\Tr(\sum_{j}\frac{h_{j}}{\lambda}\sum_{n=2}^{\infty}\frac{(-\beta\mathcal{L}_{j})^{n}}{n!})\bigr{)}^{l}+\\\
\Tr\sum_{j}\frac{h_{j}}{\lambda}\sum_{m=2}^{\infty}\frac{(-\beta\lambda\mathcal{L}_{j})^{m}}{N^{m}m!}-\Tr\sum_{n=2}^{\infty}\frac{(-\beta\lambda\mathcal{L})^{n}}{N^{n}m!}+\Tr\frac{\beta\mathcal{L}}{N}-\Tr\sum_{j}\frac{h_{j}}{\lambda}\frac{\beta\lambda\mathcal{L}_{j}}{N}}$
(74)

Using the fact that $\sum_{j}h_{j}\mathcal{L}_{j}=\mathcal{L}$, the final two
terms cancel and we have:

$\displaystyle=e^{\nicefrac{{2\beta\lambda}}{{N}}}\Biggr{|}\sum_{k=2}^{\infty}\bigl{(}-\Tr(\sum_{n=2}^{\infty}\frac{(-\beta\mathcal{L})^{n}}{N^{n}n!})\bigr{)}^{k}-\sum_{l=2}^{\infty}\bigl{(}-\Tr(\sum_{j}\frac{h_{j}}{\lambda}\sum_{n=2}^{\infty}\frac{(-\beta\mathcal{L}_{j})^{n}}{N^{n}n!})\bigr{)}^{l}$
$\displaystyle\qquad+\Tr\sum_{j}\frac{h_{j}}{\lambda}\sum_{m=2}^{\infty}\frac{(-\beta\lambda\mathcal{L}_{j})^{m}}{N^{m}m!}-\Tr\sum_{n=2}^{\infty}\frac{(-\beta\lambda\mathcal{L})^{n}}{N^{n}m!}\Biggr{|}$
(75) $\displaystyle\leq e^{\nicefrac{{2\beta\lambda}}{{N}}}\medmath$
$\displaystyle\leq
e^{\nicefrac{{2\beta\lambda}}{{N}}}\medmath{\left(\sum_{k=2}^{\infty}\bigl{(}\sum_{n=2}^{\infty}\frac{\beta^{n}\absolutevalue{\Tr{\mathcal{L}}^{n}}}{N^{n}n!}\bigr{)}^{k}+\sum_{l=2}^{\infty}\bigl{(}\sum_{j}\frac{h_{j}}{\lambda}\sum_{n=2}^{\infty}\frac{\beta^{n}\absolutevalue{\Tr{\mathcal{L}}^{n}}}{N^{n}n!}\bigr{)}^{l}+\sum_{j}\frac{h_{j}}{\lambda}\sum_{m=2}^{\infty}\frac{\beta^{m}\lambda^{m}\absolutevalue{\Tr{\mathcal{L}_{j}}^{n}}}{N^{m}m!}+\sum_{n=2}^{\infty}\frac{\beta^{n}\lambda^{n}\absolutevalue{\Tr{\mathcal{L}}^{n}}}{N^{n}m!}\right)}$
(76)

Now we once again use the property
$\absolutevalue{\Tr\mathcal{L}(\rho)}\leq\norm{\mathcal{L}(\rho)}_{1\rightarrow
1}$:

$\displaystyle\leq
e^{\nicefrac{{2\beta\lambda}}{{N}}}\medmath{\left(\sum_{k=2}^{\infty}\bigl{(}\sum_{n=2}^{\infty}\frac{\beta^{n}\norm{\mathcal{L}}_{1\rightarrow
1}^{n}}{N^{n}n!}\bigr{)}^{k}+\sum_{l=2}^{\infty}\bigl{(}\sum_{j}\frac{h_{j}}{\lambda}\sum_{n=2}^{\infty}\frac{\beta^{n}\norm{\mathcal{L}_{j}}_{1\rightarrow
1}^{n}}{N^{n}n!}\bigr{)}^{l}+\sum_{j}\frac{h_{j}}{\lambda}\sum_{m=2}^{\infty}\frac{\beta^{m}\lambda^{m}\norm{\mathcal{L}_{j}}_{1\rightarrow
1}^{m}}{N^{m}m!}+\sum_{n=2}^{\infty}\frac{\beta^{n}\lambda^{n}\norm{\mathcal{L}}_{1\rightarrow
1}^{n}}{N^{n}m!}\right)}$ (78) $\displaystyle\leq
e^{\nicefrac{{2\beta\lambda}}{{N}}}\left(\sum_{k=2}^{\infty}\bigl{(}\sum_{n=2}^{\infty}\frac{2^{n}\beta^{n}\lambda^{n}}{N^{n}n!}\bigr{)}^{k}+\sum_{l=2}^{\infty}\bigl{(}\sum_{j}\frac{h_{j}}{\lambda}\sum_{n=2}^{\infty}\frac{2^{n}\beta^{n}\lambda^{n}}{N^{n}n!}\bigr{)}^{l}+\sum_{j}\frac{h_{j}}{\lambda}\sum_{m=2}^{\infty}\frac{2^{m}\beta^{m}\lambda^{m}}{N^{m}m!}+\sum_{n=2}^{\infty}\frac{2^{n}\beta^{n}\lambda^{n}}{N^{n}m!}\right)$
(79) $\displaystyle\leq
e^{\nicefrac{{2\beta\lambda}}{{N}}}\left(2\sum_{k=2}^{\infty}\bigl{(}\frac{2\beta^{2}\lambda^{2}}{N^{2}}e^{\nicefrac{{2\lambda\beta}}{{N}}}\bigr{)}^{k}+2\frac{2\beta^{2}\lambda^{2}}{N^{2}}e^{\nicefrac{{2\lambda\beta}}{{N}}}\right)$
(80)
$\displaystyle=2\frac{\frac{4\beta^{4}\lambda^{4}}{N^{4}}e^{\nicefrac{{6\lambda\beta}}{{N}}}}{1-\frac{2\beta^{2}\lambda^{2}}{N^{2}}e^{\nicefrac{{2\lambda\beta}}{{N}}}}+2\frac{2\beta^{2}\lambda^{2}}{N^{2}}e^{\nicefrac{{4\lambda\beta}}{{N}}}$
(81)

Where in the last 3 lines we used both the tail-bound on an exponential sum;
$\sum_{n=2}^{\infty}\frac{x^{n}}{n!}\leq\frac{x^{2}}{2}e^{x}$ and the tail-sum
of the geometric series; $\sum_{n=2}^{\infty}x^{n}=\frac{x^{2}}{1-x}$. The
latter of the two requires that $|x|<1$ for convergence, so we therefore,
require
$|\frac{2\beta^{2}\lambda^{2}}{N^{2}}e^{\nicefrac{{2\lambda\beta}}{{N}}}|<1$,
which is not a strong requirement given we are interested in large $N$. Now
combining the results for $a$ and $b$ we obtain the final result for a channel
that simulates a fraction $\nicefrac{{1}}{{N}}$ of the dynamics:

$\displaystyle\norm{\mathcal{U}_{N}-\Lambda_{QD}}_{1\rightarrow 1}$
$\displaystyle\leq a+b$ (83)
$\displaystyle\leq\frac{4\beta^{2}\lambda^{2}}{N^{2}}\frac{e^{\nicefrac{{2\lambda\beta}}{{N}}}}{2-e^{\nicefrac{{2\beta\lambda}}{{N}}}}+2\frac{\frac{4\beta^{4}\lambda^{4}}{N^{4}}e^{\nicefrac{{6\lambda\beta}}{{N}}}}{1-\frac{2\beta^{2}\lambda^{2}}{N^{2}}e^{\nicefrac{{2\lambda\beta}}{{N}}}}+2\frac{2\beta^{2}\lambda^{2}}{N^{2}}e^{\nicefrac{{4\lambda\beta}}{{N}}}$
(84)

Using the fact that the norm is sub-additive under composition

$\displaystyle\norm{\mathcal{U}_{N}^{\circ N}-\Lambda_{QD}^{\circ
N}}_{1\rightarrow 1}$ $\displaystyle\leq
N\norm{\mathcal{U}_{N}-\Lambda_{QD}}_{1\rightarrow 1},$ (85)

we are left with the following result for the full-time channel

$\displaystyle\norm{\mathcal{U}-\Lambda_{QD}^{\circ N}}_{1\rightarrow 1}$
$\displaystyle\leq\frac{4\beta^{2}\lambda^{2}}{N}\frac{e^{\nicefrac{{2\lambda\beta}}{{N}}}}{2-e^{\nicefrac{{2\beta\lambda}}{{N}}}}+2\frac{\frac{4\beta^{4}\lambda^{4}}{N^{3}}e^{\nicefrac{{6\lambda\beta}}{{N}}}}{1-\frac{2\beta^{2}\lambda^{2}}{N^{2}}e^{\nicefrac{{2\lambda\beta}}{{N}}}}+2\frac{2\beta^{2}\lambda^{2}}{N}e^{\nicefrac{{4\lambda\beta}}{{N}}}.$
(86)

Next, we wish to linearize this bound for a clearer interpretation of the
error. We can do so by utilizing the conditions required to hold in the proof:

$|e^{\nicefrac{{2\beta\lambda}}{{N}}}-1|<1,$ (87)

$|\frac{2\beta^{2}\lambda^{2}}{N^{2}}e^{\nicefrac{{2\lambda\beta}}{{N}}}|<1.$
(88)

The first condition (87) gives $e^{\nicefrac{{2\beta\lambda}}{{N}}}\in(0,2)$
which implies $\frac{2\beta\lambda}{N}\in(-\infty,\ln{2})$, and given that we
require this quantity necessarily be positive we have
$\frac{2\beta\lambda}{N}\in[0,\ln{2})$. Now, our second condition (88), given
the positivity constraint, further implies $\frac{x^{2}}{2}e^{x}\in[0,1)$ with
$x=\nicefrac{{2\beta\lambda}}{{N}}$. Using max value obtained from Equation
(87), that $x<\ln{2}$, Equation (88) reduces to $\frac{x^{2}}{2}e^{x}\leq
2\ln{2}\in[0,1)$. Given that Equation (88) is a product of two monotonic
functions in the domain $(0,\infty)$, any value of $x$ less than the maximum
obtained in the first condition also satisfies the second. Therefore, going
forward we focus only on satisfying the constraint in Equation (87).

Going forward in linearizing Equation 86 we will take $\lambda$ and $N$ to be
constant and consider variable $\beta$, and continue by upper bounding each
term. Beginning with the first term we have:

$\frac{4\beta^{2}\lambda^{2}}{N}\frac{e^{\nicefrac{{2\lambda\beta}}{{N}}}}{2-e^{\nicefrac{{2\beta\lambda}}{{N}}}}=\frac{2\beta^{2}\lambda^{2}}{N}\frac{1}{y-\frac{1}{2}},\;\;\mathbf{s.t.}\;y(\beta)=e^{-\nicefrac{{2\lambda\beta}}{{N}}}\;\text{and}\;\frac{2}{3}\leq
y(\beta)\leq 1$ (89)

Now we rewrite the aforementioned constraint $1\leq
e^{\nicefrac{{2\lambda\beta}}{{N}}}<2$, to be $1\leq
e^{\nicefrac{{2\lambda\beta}}{{N}}}\leq\frac{3}{2}$ to prevent one from
approaching a divergence in the proof arbitrarily closely. Utilizing the
following property:

$f(b)=f(a)+\int_{a}^{b}\frac{df}{dy}dy\leq
f(a)+(b-a)\max_{y\in[a,b]}{\absolutevalue{\frac{df(y)}{dy}}},$ (90)

with $f(\beta)=\frac{1}{y(\beta)-\frac{1}{2}}$ Equation 89 is upper bounded in
the following way:

$\frac{4\beta^{2}\lambda^{2}}{N}\frac{4\beta^{2}\lambda^{2}}{N}e^{\nicefrac{{4\beta\lambda}}{{N}}}\leq\frac{4\beta^{2}\lambda^{2}}{N}(3+12\frac{\lambda}{N})$
(91)

Moving on to the third term in Equation (86), we have:

$\frac{4\beta^{2}\lambda^{2}}{N}e^{\nicefrac{{4\beta\lambda}}{{N}}}=\frac{4\beta^{2}\lambda^{2}}{N}y^{2},\;\;\mathbf{s.t.}\;\;y=e^{\nicefrac{{2\beta\lambda}}{{N}}}\;\text{and}\;1\leq
y(\beta)\leq\frac{3}{2},$ (92)

and using the same integral bound we obtain:

$\frac{4\beta^{2}\lambda^{2}}{N}e^{\nicefrac{{4\beta\lambda}}{{N}}}\leq\frac{4\beta^{2}\lambda^{2}}{N}(1+3\frac{\lambda}{N}).$
(93)

Finally, we are left to bound the second term in Equation 86, which we can
write in the following way:

$\left(\frac{4\beta^{2}\lambda^{2}}{N}e^{\nicefrac{{4\lambda\beta}}{{N}}}\right)\frac{\frac{4\beta^{2}\lambda^{2}}{N^{2}}e^{\nicefrac{{2\lambda\beta}}{{N}}}}{2-\frac{4\beta^{2}\lambda^{2}}{N^{2}}e^{\nicefrac{{2\lambda\beta}}{{N}}}}=2\left(\frac{4\beta^{2}\lambda^{2}}{N}e^{\nicefrac{{4\lambda\beta}}{{N}}}\right)\frac{y\ln{y}}{1-y\ln{y}},\;\;\mathbf{s.t.}\;y(\beta)=e^{\nicefrac{{2\lambda\beta}}{{N}}}\;\text{and}\;1\leq
y(\beta)\leq\frac{3}{2}.$ (94)

The function $f(y)=\frac{y\ln{y}}{1-y\ln{y}}$ is monotonic over the
constrained domain and can be upper bounded by inserting the largest value
that $\ln{y}$ achieves. Also, notice the term in parentheses is precisely the
3rd term in Equation 86 that we bounded above. Using these two observations:

$\displaystyle
2\left(\frac{4\beta^{2}\lambda^{2}}{N}e^{\nicefrac{{4\lambda\beta}}{{N}}}\right)\frac{y\ln{y}}{1-y\ln{y}}$
$\displaystyle\leq\frac{4\beta^{2}\lambda^{2}}{N}2\left(1+3\frac{\lambda}{N}\right)\frac{y\ln{\frac{3}{2}}}{1-y\ln{\frac{3}{2}}}$
(95)
$\displaystyle\leq\frac{4\beta^{2}\lambda^{2}}{N}2\ln{\frac{3}{2}}\left(1+3\frac{\lambda}{N}\right)\left(\frac{1}{\frac{2}{3}-\ln{\frac{3}{2}}}+\frac{2\lambda}{3N}\frac{1}{(\frac{2}{3}-\ln{\frac{3}{2}})^{2}}\right)$
(96)

where in the last step we once again used the integral identity following the
same procedure as we did for the first term. Now that we have linearized each
term, we combine all the results to obtain:

$\norm{\mathcal{U}-\Lambda_{QD}^{\circ N}}_{1\rightarrow
1}\leq\frac{\beta^{2}\lambda^{2}}{N}\left(a+b\frac{\lambda}{N}+c\frac{\lambda^{2}}{N^{2}}\right)$
(97)

where $a\approx 28.41845,\;b\approx 128.95110\;\text{and}\;c\approx 95.08717$.
If we make the promise that $\frac{\lambda}{N}\leq 0.01$, the bound simplifies
to

$\norm{\mathcal{U}-\Lambda_{QD}^{\circ N}}_{1\rightarrow
1}\leq\frac{C\beta^{2}\lambda^{2}}{N}$ (98)

where $C\approx 29.71747$, completing the proof. ∎

This condition is not very restrictive on the imaginary time $\beta$.
Considering the condition given in Equation (87), we have
$\beta<\frac{\ln{2}N}{2\lambda}\approx 34.657$, so $\beta<10\pi$, meaning that
given this promise on the ratio $\frac{\lambda}{N}$, the bound holds for a
simulation time that allows the phase to oscillate up to 5 times if we were
dealing with real-time.

###  A.2 Imaginary Time Trotter-Suzuki Formulas

Following the procedure in [22], we must first convert the existing bound to a
bound on channels and then bound the error with the renormalizing operations.

See III.2

###### Proof.

First consider the induced 1-norm of the difference between the channels after
a single iteration or time-step $r$. As in the previous proof, absolute values
containing a channel operation are understood to be maximized over all inputs
$\rho:\norm{\rho}_{1}=1$:

$\displaystyle\norm{\frac{\mathcal{U}{(\rho,\nicefrac{{\beta}}{{r}},\nicefrac{{\beta}}{{r}})}}{\Tr(\mathcal{U}{(\rho,\nicefrac{{\beta}}{{r}},\nicefrac{{\beta}}{{r}})})}-\frac{\mathcal{T}^{2k}{(\rho,\nicefrac{{\beta}}{{r}},\nicefrac{{\beta}}{{r}})}}{\Tr(\mathcal{T}^{2k}{(\rho,\nicefrac{{\beta}}{{r}},\nicefrac{{\beta}}{{r}})})}}_{1\rightarrow
1}\leq\norm{\frac{\mathcal{U}{(\rho,\nicefrac{{\beta}}{{r}})}-\mathcal{T}^{2k}{(\rho,\nicefrac{{\beta}}{{r}})}}{\Tr\mathcal{U}{(\rho,\nicefrac{{\beta}}{{r}})}}}_{1\rightarrow
1}$
$\displaystyle\qquad+\norm{\mathcal{T}^{2k}{(\rho,\nicefrac{{\beta}}{{r}})}\left(\frac{1}{\Tr\mathcal{T}^{2k}{(\rho,\nicefrac{{\beta}}{{r}})}}-\frac{1}{\Tr\mathcal{U}{(\rho,\nicefrac{{\beta}}{{r}})}}\right)}_{1\rightarrow
1}$ (99)
$\displaystyle=\absolutevalue{\frac{1}{\Tr\mathcal{U}{(\rho,\nicefrac{{\beta}}{{r}})}}}\norm{\mathcal{U}{(\rho,\nicefrac{{\beta}}{{r}})}-\mathcal{T}^{2k}{(\rho,\nicefrac{{\beta}}{{r}})}}_{1\rightarrow
1}+\norm{\mathcal{T}^{2k}{(\rho,\nicefrac{{\beta}}{{r}})}}_{1\rightarrow
1}\absolutevalue{\frac{\Tr(\mathcal{U}{(\rho,\nicefrac{{\beta}}{{r}})}-\mathcal{T}^{2k}{(\rho,\nicefrac{{\beta}}{{r}})})}{\Tr(\mathcal{U}{(\rho,\nicefrac{{\beta}}{{r}})})\Tr(\mathcal{T}^{2k}{(\rho,\nicefrac{{\beta}}{{r}})})}}$
(100)
$\displaystyle\leq\absolutevalue{\frac{1}{\Tr\mathcal{U}{(\rho,\nicefrac{{\beta}}{{r}})}}}\norm{\mathcal{U}{(\rho,\nicefrac{{\beta}}{{r}})}-\mathcal{T}^{2k}{(\rho,\nicefrac{{\beta}}{{r}})}}_{1\rightarrow
1}$
$\displaystyle\qquad+\norm{\mathcal{T}^{2k}{(\rho,\nicefrac{{\beta}}{{r}})}}_{1\rightarrow
1}\norm{\mathcal{U}{(\rho,\nicefrac{{\beta}}{{r}})}-\mathcal{T}^{2k}{(\rho,\nicefrac{{\beta}}{{r}})}}_{1\rightarrow
1}\absolutevalue{\frac{1}{\Tr(\mathcal{U}{(\rho,\nicefrac{{\beta}}{{r}})})\Tr(\mathcal{T}^{2k}{(\rho,\nicefrac{{\beta}}{{r}})})}}$
(101)
$\displaystyle=\norm{\mathcal{U}{(\rho,\nicefrac{{\beta}}{{r}})}-\mathcal{T}^{2k}{(\rho,\nicefrac{{\beta}}{{r}})}}_{1\rightarrow
1}\left(\frac{\norm{\mathcal{T}^{2k}{(\rho,\nicefrac{{\beta}}{{r}})}}_{1\rightarrow
1}}{\absolutevalue{\Tr(\mathcal{U}{(\rho,\nicefrac{{\beta}}{{r}})})\Tr(\mathcal{T}^{2k}{(\rho,\nicefrac{{\beta}}{{r}})})}}+\frac{1}{\absolutevalue{\Tr\mathcal{U}{(\rho,\nicefrac{{\beta}}{{r}})}}}\right)$
(102)

where in the second line we simply cross-multiply the traces of the outputs of
the channels (since they are just numbers), and in third line we use the fact
that the absolute value of the trace is trivially upper bounded by the diamond
norm. Going forward we need to lower bound the traces in the denominator,
which can be done using the well-known von Neumann trace inequality:
$\sum_{i=1}^{N}\sigma^{A}_{i}\sigma^{B}_{N-i+1}\leq\Tr(AB)$ where
$\sigma^{(A,B)}_{i}$ represent the singular values of the matrix $A$ or $B$
respectfully, and they are ordered such that $\sigma_{i}\leq\sigma_{i+1}$. The
absolute values in the denominator also need not be written as these traces
are positive by definition. For the exact evolution channel we have:

$\displaystyle\Tr\mathcal{U}{(\rho,\nicefrac{{\beta}}{{r}})}$
$\displaystyle=\Tr(e^{-\beta H}\rho e^{-\nicefrac{{\beta}}{{r}}H})$ (103)
$\displaystyle=\Tr(e^{-2\nicefrac{{\beta}}{{r}}H}\rho)$ (104)
$\displaystyle\geq\sum_{i}\sigma_{i}(e^{-2\nicefrac{{\beta}}{{r}}H})\sigma_{N-i+1}(\rho)$
(105) $\displaystyle\geq\sigma_{min}(e^{-2\nicefrac{{\beta}}{{r}}H})\Tr(\rho)$
(106) $\displaystyle\geq e^{-2\nicefrac{{\beta}}{{r}}\norm{H}},$ (107)

where we use the fact that the trace of $\rho$ is 1, along with the fact that
$\rho$ is positive semi-definite and hermitian, such that the sum of its
singular values is equal to its trace. Next, we follow a similar procedure in
dealing with the trace of the output of the Trotter-Suzuki channel:

$\displaystyle\Tr\mathcal{T}^{2k}{(\rho,\nicefrac{{\beta}}{{r}})}$
$\displaystyle=\Tr(S_{2k}(\nicefrac{{\beta}}{{r}})\rho
S_{2k}^{\dagger}(\nicefrac{{\beta}}{{r}}))$ (108)
$\displaystyle\geq\sigma_{min}\left(S_{2k}^{\dagger}(\nicefrac{{\beta}}{{r}})S_{2k}(\nicefrac{{\beta}}{{r}})\right)\Tr\rho$
(109)
$\displaystyle=\sigma_{min}\left(S_{2k}^{\dagger}(\nicefrac{{\beta}}{{r}})S_{2k}(\nicefrac{{\beta}}{{r}})\right).$
(110)

In order to understand how this term behaves with respect to
$\nicefrac{{\beta}}{{r}}$ and $H$, we require more careful analysis. To do
this it will be useful to introduce some more convenient notation. Just like
an arbitrary product formula can be defined in the following way (following
the notation of [10]):

$S_{2k}(\nicefrac{{\beta}}{{r}})=\prod_{\upsilon=1}^{\Upsilon}\prod_{l=1}^{L}e^{-\nicefrac{{\beta}}{{r}}s_{\upsilon,l}H_{\pi_{\upsilon(l)}}}$
(111)

Where $s$ controls the pre-factor that modifies the time-step and $\pi$
controls the ordering of the $H$ summands in each stage of the formula. We
utilize the following property of singular values: given two matrices
$A,\;B\in\mathbb{C}$ the following holds
$\sigma_{n}(A)\sigma_{i}(B)\leq\sigma_{i}(AB)\leq\sigma_{1}(A)\sigma_{i}(B)$,
given that singular values are ordered such that
$\sigma_{1}\geq\sigma_{2}\geq...\geq\sigma_{n}$. Also, note that since the
individual evolution operators are no longer anti-hermitian, coupled with the
fact that Trotter-Suzuki formulas are symmetric, $S_{2k}^{\dagger}=S_{2k}$.

$\displaystyle\Tr\mathcal{T}^{2k}{(\rho,\nicefrac{{\beta}}{{r}})}$
$\displaystyle\geq\sigma_{min}\left(S_{2k}^{2}(\nicefrac{{\beta}}{{r}})\right)$
(112)
$\displaystyle\geq\sigma^{2}_{min}\left(\prod_{\upsilon=1}^{\Upsilon}\prod_{l=1}^{L}e^{-\nicefrac{{\beta}}{{r}}s_{\upsilon,l}H_{\pi_{\upsilon(l)}}}\right)$
(113)
$\displaystyle\geq\left(\prod_{\upsilon=1}^{\Upsilon}\prod_{l=1}^{L}\sigma_{min}\left(e^{-\nicefrac{{\beta}}{{r}}s_{\upsilon,l}H_{\pi_{\upsilon(l)}}}\right)\right)^{2}$
(114)
$\displaystyle\geq\left(\prod_{\upsilon=1}^{\Upsilon}\prod_{l=1}^{L}e^{-\nicefrac{{\beta}}{{r}}s_{\upsilon,l}\norm{H_{\upsilon(l)}}}\right)^{2}$
(115)
$\displaystyle=\left(e^{\sum_{\upsilon=1}^{\Upsilon}\sum_{l=1}^{L}-\nicefrac{{\beta}}{{r}}s_{\upsilon,l}\norm{H_{\upsilon(l)}}}\right)^{2}$
(116) $\displaystyle=e^{-2\nicefrac{{\beta}}{{r}}\lambda}$ (117)

where each inequality step is a repeated application of the aforementioned
property. In the last step, we complete the using the property of Trotter-
Suzuki formulas that $\sum_{\upsilon,l}s_{\upsilon,l}=1$ and the definition
$\sum_{l}\norm{H_{l}}=\lambda$.

We also need to deal with the most important term which is the actual induced
1-norm of the difference between the channels. To do this, we will massage
this term into a form where we can apply a result from [10]. To do this we
follow the same procedure from Section II.7:

$\displaystyle\norm{\mathcal{U}(\beta)-\mathcal{T}^{2k}{(\beta)}}_{1\rightarrow
1}:=$
$\displaystyle\max_{\rho:\norm{\rho}=1}\norm{\mathcal{U}(\beta)-\mathcal{T}^{2k}{(\beta)}}_{1}$
(118) $\displaystyle=$
$\displaystyle\max_{\rho:\norm{\rho}_{1}=1}\norm{e^{-H\beta}\rho
e^{-H\beta}-S_{2k}(\beta)\rho S_{2k}(\beta)}_{1}$ (119) $\displaystyle\leq$
$\displaystyle\max_{\rho:\norm{\rho}_{1}=1}\norm{e^{-H\beta}\rho
e^{-H\beta}-e^{-H\beta}\rho S_{2k}(\beta)}_{1}$ (120)
$\displaystyle+\max_{\rho:\norm{\rho}_{1}\leq 1}\norm{e^{-H\beta}\rho
S_{2k}(\beta)-S_{2k}(\beta)\rho S_{2k}(\beta)}_{1}$ (121) $\displaystyle=$
$\displaystyle\max_{\rho:\norm{\rho}_{1}=1}\norm{e^{-\beta
H}\rho\left(e^{-H\beta}-S_{2k}(\beta)\right)}_{1}+\max_{\rho:\norm{\rho}_{1}=1}\norm{\left(e^{-H\beta}-S_{2k}(\beta)\right)\rho
S_{2k}(\beta)}_{1}$ (122) $\displaystyle\leq$
$\displaystyle\left(\norm{e^{-\beta
H}}_{\infty}+\norm{S_{2k}(\beta)}_{\infty}\right)\norm{e^{-H\beta}-S_{2k}(\beta)}_{\infty}\max_{\rho:\norm{\rho}_{1}\leq
1}\norm{\rho}_{1}$ (123) $\displaystyle=$ $\displaystyle
2e^{\beta\lambda}\norm{e^{-H\beta}-S_{2k}(\beta)}_{\infty}$ (124)
$\displaystyle=$ $\displaystyle
2e^{\beta\lambda}\norm{U(\beta)-S_{2k}(\beta)}_{\infty}.$ (125)

Next we apply the following bound from Appendix E in [10]:

$\norm{U(\nicefrac{{\beta}}{{r}})-S_{2k}(\nicefrac{{\beta}}{{r}})}_{\infty}\leq
2\Upsilon^{2k+1}\frac{\alpha_{comm}(H,2k)}{(2k+1)!}\frac{\beta^{2k+1}}{r^{2k+1}}e^{4\frac{\beta}{r}\Upsilon\lambda}$
(126)

Where $\Upsilon$ and $\lambda$ are again the number of stages of the product
formula and the sum of the spectral norms of the Hamiltonian summands
respectively. Now we are simply left with a single term to bound in the
diamond norm of the Trotter-Suzuki channel. We can similarly define the
Trotter-Suzuki channel as a composition of multiple stages of channels, each
containing compositions of evolution channels corresponding to each
Hamiltonian summand. We write this in the following way:

$\mathcal{T}^{2k}{(\nicefrac{{\beta}}{{r}})}\equiv\bigcirc_{\upsilon=1}^{\Upsilon}\bigcirc_{l=1}^{L}e^{-\frac{\beta}{r}s_{\upsilon,l}\mathcal{L}_{\pi_{\upsilon(l)}}},$
(127)

where we once again have the Liouvillian generator for imaginary time
$\mathcal{L}(\rho)=\\{H,\rho\\}$, and all other machinery is defined in the
same was as above. In this way it is more natural to proceed in using the sub-
multiplicative property of the norm:

$\displaystyle\norm{\mathcal{T}^{2k}{(\nicefrac{{\beta}}{{r}})}}_{1\rightarrow
1}$
$\displaystyle\leq\prod_{\upsilon=1}^{\Upsilon}\prod_{l=1}^{L}\norm{e^{-\frac{\beta}{r}s_{\upsilon,l}\mathcal{L}_{\pi_{\upsilon(l)}}}}_{1\rightarrow
1}$ (128)
$\displaystyle\leq\prod_{\upsilon=1}^{\Upsilon}\prod_{l=1}^{L}e^{2\frac{\beta}{r}s_{\upsilon,l}\norm{H_{\pi_{\upsilon(l)}}}}$
(129)
$\displaystyle=e^{\sum_{\upsilon=1}^{\Upsilon}\sum_{l=1}^{L}2\frac{\beta}{r}s_{\upsilon,l}\norm{H_{\upsilon(l)}}}$
(130) $\displaystyle=e^{2\frac{\beta}{r}\lambda}.$ (131)

In the second line we use the Taylor expansion of the exponential of the
Liouvillian to apply the diamond norm to the exponent as seen in the QDrift
proof. We make use of the aforementioned fact that
$\norm{\mathcal{L}_{l}}_{1\rightarrow 1}\leq 2h_{l}$ where $h_{l}$ is the
spectral norm of the $H_{l}$ Hamiltonian term. Now, combining all of our
results we are left with the following:

$\norm{\frac{\mathcal{U}{(\rho,\beta)}}{\Tr(\mathcal{U}{(\rho,\beta)})}-\left(\frac{\mathcal{T}^{2k}{(\rho,\nicefrac{{\beta}}{{r}})}}{\Tr(\mathcal{T}^{2k}{(\rho,\nicefrac{{\beta}}{{r}})})}\right)^{\circ
r}}_{1\rightarrow 1}\leq
4\Upsilon^{2k+1}\frac{\alpha_{comm}(H,2k)}{(2k+1)!}\frac{\beta^{2k+1}}{r^{2k+1}}e^{4\frac{\beta}{r}\Upsilon\lambda}e^{\beta\lambda/r}\left(e^{4\frac{\beta}{r}\lambda}e^{2\frac{\beta}{r}\norm{H}}+e^{2\frac{\beta}{r}\norm{H}}\right)$
(132)

and we wish to compare the exact evolution to the channel after $r$
iterations. In other words, we use the property $\norm{X^{\circ r}-Y^{\circ
r}}_{1\rightarrow 1}\leq r\norm{X-Y}_{1\rightarrow 1}$, proven earlier (for
the diamond norm, but it also holds here) in Equation 35, to obtain the
following:

$\norm{\frac{\mathcal{U}{(\rho,\beta)}}{\Tr(\mathcal{U}{(\rho,\beta)})}-\left(\frac{\mathcal{T}^{2k}{(\rho,\nicefrac{{\beta}}{{r}})}}{\Tr(\mathcal{T}^{2k}{(\rho,\nicefrac{{\beta}}{{r}})})}\right)^{\circ
r}}_{1\rightarrow 1}\leq
4\Upsilon^{2k+1}\frac{\alpha_{comm}(H,2k)}{(2k+1)!}\frac{\beta^{2k+1}}{r^{2k}}e^{4\frac{\beta}{r}\Upsilon\lambda}e^{\beta\lambda/r}\left(e^{4\frac{\beta}{r}\lambda}e^{2\frac{\beta}{r}\norm{H}}+e^{2\frac{\beta}{r}\norm{H}}\right)$
(133)

Similar to the QDrift case, we wish to eliminate the exponential factors to
give a more straightforward interpretation of scaling of the error bound.
However, in this case the proof did not impose restrictions on the size of
these exponentials. Therefore, we instead upper-bound each exponential term in
a common way and then make on observation about the scaling of $r$, to show
that the exponential factors are $\mathcal{O}\left(1\right)$. We begin by
using the fact that $\norm{H}\leq\lambda$, and
$e^{\nicefrac{{\beta\lambda}}{{r}}}\leq
e^{\Upsilon\nicefrac{{\beta\lambda}}{{r}}}$ to obtain the following:

$\norm{\frac{\mathcal{U}{(\rho,\beta)}}{\Tr(\mathcal{U}{(\rho,\beta)})}-\left(\frac{\mathcal{T}^{2k}{(\rho,\nicefrac{{\beta}}{{r}})}}{\Tr(\mathcal{T}^{2k}{(\rho,\nicefrac{{\beta}}{{r}})})}\right)^{\circ
r}}_{1\rightarrow 1}\leq
4\Upsilon^{2k+1}\frac{\alpha_{comm}(H,2k)}{(2k+1)!}\frac{\beta^{2k+1}}{r^{2k}}\left(e^{\frac{11\Upsilon\beta\lambda}{r}}+e^{\frac{7\Upsilon\beta\lambda}{r}}\right)$
(134)

Now, let us posit that
$r\in\Omega(\Upsilon^{a}\beta^{b}\lambda^{c})\;|\;a,b,c\geq 1$. Focusing on
the terms in parentheses above, the terms in the exponents are positive by
definition, so we have
$\left(e^{\nicefrac{{10\Upsilon\beta\lambda}}{{r}}}+e^{\nicefrac{{6\Upsilon\beta\lambda}}{{r}}}\right)\geq
1$. If the posited scaling holds, then
$\left(e^{\nicefrac{{11\Upsilon\beta\lambda}}{{r}}}+e^{\nicefrac{{7\Upsilon\beta\lambda}}{{r}}}\right)\in\mathcal{O}\left(1\right)$
such that this term is asymptotically upper-bounded by some constant $C$. In
fact, if $a,b$ or $c>1$, then C=1 in the infinite limit. From [22], we use the
fact that $\alpha_{comm}(H,2k)\leq 2^{2k}\lambda^{2k+1}$ to verify this
intuition.

$\norm{\frac{\mathcal{U}{(\rho,\beta)}}{\Tr(\mathcal{U}{(\rho,\beta)})}-\left(\frac{\mathcal{T}^{2k}{(\rho,\nicefrac{{\beta}}{{r}})}}{\Tr(\mathcal{T}^{2k}{(\rho,\nicefrac{{\beta}}{{r}})})}\right)^{\circ
r}}_{1\rightarrow 1}\leq
4^{2k+1}\Upsilon^{2k+1}\frac{\lambda^{2k+1}}{(2k+1)!}\frac{\beta^{2k+1}}{r^{2k}}C<\epsilon,$
(135)

where solving for $r$ in terms of $\epsilon$ gives the following:

$\left(4\Upsilon\lambda\beta\right)^{1+\frac{1}{2k}}\left(\frac{C}{\epsilon(2k+1)!}\right)^{\frac{1}{2k}}<r.$
(136)

Asymptotically, this gives the following asymptotic bound:

$r\in\Omega(\frac{\Upsilon^{1+\frac{1}{2k}}\lambda^{1+\frac{1}{2k}}\beta^{1+\frac{1}{2k}}}{\epsilon^{\frac{1}{2k}}}).$
(137)

Here, we enforce $\epsilon\leq 1$, which then confirms the initial intuition
about the asymptotic scaling of $r$, such that that
$r\in\Omega(\Upsilon^{a}\beta^{b}\lambda^{c})$, where
$a=b=c=1+\nicefrac{{1}}{{2k}}$, and the exponential term vanishes in the
infinite limit, yielding the result below:

$\norm{\frac{\mathcal{U}{(\rho,\beta)}}{\Tr(\mathcal{U}{(\rho,\beta)})}-\left(\frac{\mathcal{T}^{2k}{(\rho,\nicefrac{{\beta}}{{r}})}}{\Tr(\mathcal{T}^{2k}{(\rho,\nicefrac{{\beta}}{{r}})})}\right)^{\circ
r}}_{1\rightarrow
1}\in\mathcal{O}\left(\Upsilon^{2k+1}\frac{\alpha_{comm}(H,2k)}{(2k+1)!}\frac{\beta^{2k+1}}{r^{2k}}\right),$
(138)

thus completing the proof. ∎

###  A.3 Imaginary Time Composite Channels

See III.3

###### Proof.

We begin by expanding the imaginary time composite channel
$\mathcal{X}(\rho,\beta)$ in terms of familiar quantities that have been
previously bounded. That being said, this proof relies almost entirely on
Theorems III.1 and III.2. Throughout, notation is used such that
$\mathcal{U}{(\rho,\beta)}_{A}$ represents an exact evolution channel under
the Hamiltonian $A$ from the the theorem statement. It is implied throughout
that terms in $A$ are simulated by the Trotter-Suzuki channel, while terms in
$B$ are simulated by QDrift. Therefore, we write
$\mathcal{T}_{A}^{2k}{(\rho,\nicefrac{{\beta}}{{r}})}$ and
$\mathcal{Q}_{B}{(\rho,\nicefrac{{\beta}}{{r}})}$ as Trotter and QDrift
channels generated by the Hamiltonians $A$ and $B$ for clearer notation, and
once again we initial consider a time-step $\nicefrac{{\beta}}{{r}}$.
Expanding the composite channel we have:

$\norm{\frac{\mathcal{U}{(\rho,\nicefrac{{\beta}}{{r}})}}{\Tr\mathcal{U}{(\rho,\nicefrac{{\beta}}{{r}})}}-\frac{\mathcal{X}^{2k}(\rho,\nicefrac{{\beta}}{{r}})}{\Tr\mathcal{X}^{2k}(\rho,\nicefrac{{\beta}}{{r}})}}_{1\rightarrow
1}=\norm{\frac{\mathcal{U}_{B}(\rho,\nicefrac{{\beta}}{{r}})\circ\mathcal{U}_{A}(\rho,\nicefrac{{\beta}}{{r}})+E_{A,B}(\nicefrac{{\beta}}{{r}})}{\Tr\mathcal{U}(\rho,\nicefrac{{\beta}}{{r}})}-\frac{\mathcal{Q}_{B}{(\rho,\nicefrac{{\beta}}{{r}})}\circ\mathcal{T}_{A}^{2k}{(\rho,\nicefrac{{\beta}}{{r}})}}{\Tr\mathcal{Q}_{B}{(\rho,\nicefrac{{\beta}}{{r}})}\circ\mathcal{T}_{A}^{2k}{(\rho,\nicefrac{{\beta}}{{r}})}}}_{1\rightarrow
1},$ (139)

where $E_{A,B}(\nicefrac{{\beta}}{{r}})$ is an error term that comes from
partitioning the exact evolution $\mathcal{U}$ into two evolution channels
over $A$ and $B$. We leave the denominator un-expanded for convenience. As
before, we can add and subtract a convenient term to give us terms with common
denominators that are more straightforward to bound, and then apply the
triangle inequality.

$\displaystyle=\norm{\frac{\mathcal{U}_{B}(\rho,\nicefrac{{\beta}}{{r}})\circ\mathcal{U}_{A}(\rho,\nicefrac{{\beta}}{{r}})-\mathcal{Q}_{B}{(\rho,\nicefrac{{\beta}}{{r}})}\circ\mathcal{T}_{A}^{2k}{(\rho,\nicefrac{{\beta}}{{r}})}}{\Tr\mathcal{U}(\rho,\nicefrac{{\beta}}{{r}})}}_{1\rightarrow
1}$
$\displaystyle+\norm{\frac{E_{A,B}(\nicefrac{{\beta}}{{r}})}{\Tr\mathcal{U}(\rho,\nicefrac{{\beta}}{{r}})}}_{1\rightarrow
1}+\norm{\frac{\mathcal{Q}_{B}{(\rho,\nicefrac{{\beta}}{{r}})}\circ\mathcal{T}_{A}^{2k}{(\rho,\nicefrac{{\beta}}{{r}})}}{\Tr\mathcal{U}(\rho,\nicefrac{{\beta}}{{r}})}-\frac{\mathcal{Q}_{B}{(\rho,\nicefrac{{\beta}}{{r}})}\circ\mathcal{T}_{A}^{2k}{(\rho,\nicefrac{{\beta}}{{r}})}}{\Tr\mathcal{Q}_{B}{(\rho,\nicefrac{{\beta}}{{r}})}\circ\mathcal{T}_{A}^{2k}{(\rho,\nicefrac{{\beta}}{{r}})}}}_{1\rightarrow
1}.$ (140)

Now we proceed by bounding each of these terms individually. Starting with the
first term:

$\displaystyle\norm{\frac{\mathcal{U}_{B}(\rho,\nicefrac{{\beta}}{{r}})\circ\mathcal{U}_{A}(\rho,\nicefrac{{\beta}}{{r}})-\mathcal{Q}_{B}{(\rho,\nicefrac{{\beta}}{{r}})}\circ\mathcal{T}_{A}^{2k}{(\rho,\nicefrac{{\beta}}{{r}})}}{\Tr\mathcal{U}(\rho,\nicefrac{{\beta}}{{r}})}}_{1\rightarrow
1}\leq\frac{\norm{\mathcal{U}_{B}(\rho,\nicefrac{{\beta}}{{r}})-\mathcal{Q}_{B}{(\rho,\nicefrac{{\beta}}{{r}})}}_{1\rightarrow
1}+\norm{\mathcal{U}_{A}(\rho,\nicefrac{{\beta}}{{r}})-\mathcal{T}_{A}^{2k}{(\rho,\nicefrac{{\beta}}{{r}})}}_{1\rightarrow
1}}{\Tr\mathcal{U}(\rho,\nicefrac{{\beta}}{{r}})}$ (141)
$\displaystyle\leq\left(\norm{\mathcal{U}_{B}(\rho,\nicefrac{{\beta}}{{r}})-\mathcal{Q}_{B}{(\rho,\nicefrac{{\beta}}{{r}})}}_{1\rightarrow
1}+\norm{\mathcal{U}_{A}(\rho,\nicefrac{{\beta}}{{r}})-\mathcal{T}_{A}^{2k}{(\rho,\nicefrac{{\beta}}{{r}})}}_{1\rightarrow
1}\right)e^{2\nicefrac{{\norm{H}\beta}}{{r}}}$ (142)

where in the second line we use the trace bound found in the proof of Theorem
III.2. Each of the diamond norms remaining above have been bounded in Theorems
III.1 and III.2, and we will plug in said results at the end of the proof.
Proceeding to the second term we have:

$\displaystyle\norm{\frac{\mathcal{Q}_{B}{(\rho,\nicefrac{{\beta}}{{r}})}\circ\mathcal{T}_{A}^{2k}{(\rho,\nicefrac{{\beta}}{{r}})}}{\Tr\mathcal{U}(\rho,\nicefrac{{\beta}}{{r}})}-\frac{\mathcal{Q}_{B}{(\rho,\nicefrac{{\beta}}{{r}})}\circ\mathcal{T}_{A}^{2k}{(\rho,\nicefrac{{\beta}}{{r}})}}{\Tr\mathcal{Q}_{B}{(\rho,\nicefrac{{\beta}}{{r}})}\circ\mathcal{T}_{A}^{2k}{(\rho,\nicefrac{{\beta}}{{r}})}}}_{1\rightarrow
1}$ (143)
$\displaystyle\leq\norm{\mathcal{Q}_{B}{(\rho,\nicefrac{{\beta}}{{r}})}}_{1\rightarrow
1}\norm{\mathcal{T}_{A}^{2k}{(\rho,\nicefrac{{\beta}}{{r}})}}_{1\rightarrow
1}\absolutevalue{\frac{\Tr\left(\mathcal{U}(\rho,\nicefrac{{\beta}}{{r}})-\mathcal{Q}_{B}{(\rho,\nicefrac{{\beta}}{{r}})}\circ\mathcal{T}_{A}^{2k}{(\rho,\nicefrac{{\beta}}{{r}})}\right)}{\Tr\left(\mathcal{U}(\rho,\nicefrac{{\beta}}{{r}})\right)\Tr\left(\mathcal{Q}_{B}{(\rho,\nicefrac{{\beta}}{{r}})}\circ\mathcal{T}_{A}^{2k}{(\rho,\nicefrac{{\beta}}{{r}})}\right)}}$
(144)
$\displaystyle=\norm{\mathcal{Q}_{B}{(\rho,\nicefrac{{\beta}}{{r}})}}_{1\rightarrow
1}\norm{\mathcal{T}_{A}^{2k}{(\rho,\nicefrac{{\beta}}{{r}})}}_{1\rightarrow
1}\frac{\absolutevalue{\Tr\left(\mathcal{U}_{B}(\rho,\nicefrac{{\beta}}{{r}})\circ\mathcal{U}_{A}(\rho,\nicefrac{{\beta}}{{r}})+E_{A,B}(\nicefrac{{\beta}}{{r}})-\mathcal{Q}_{B}{(\rho,\nicefrac{{\beta}}{{r}})}\circ\mathcal{T}_{A}^{2k}{(\rho,\nicefrac{{\beta}}{{r}})}\right)}}{\Tr\left(\mathcal{U}(\rho,\nicefrac{{\beta}}{{r}})\right)\Tr\left(\mathcal{Q}_{B}{(\rho,\nicefrac{{\beta}}{{r}})}\circ\mathcal{T}_{A}^{2k}{(\rho,\nicefrac{{\beta}}{{r}})}\right)}$
(145) $\displaystyle\leq
e^{\frac{2\beta\lambda_{B}}{r}}e^{\frac{2\beta\lambda_{A}}{r}}e^{\frac{2\beta\norm{H}}{r}}\frac{\norm{\mathcal{U}_{B}(\rho,\nicefrac{{\beta}}{{r}})\circ\mathcal{U}_{A}(\rho,\nicefrac{{\beta}}{{r}})+E_{A,B}(\nicefrac{{\beta}}{{r}})-\mathcal{Q}_{B}{(\rho,\nicefrac{{\beta}}{{r}})}\circ\mathcal{T}_{A}^{2k}{(\rho,\nicefrac{{\beta}}{{r}})}}_{1\rightarrow
1}}{\Tr\left(\mathcal{Q}_{B}{(\rho,\nicefrac{{\beta}}{{r}})}\circ\mathcal{T}_{A}^{2k}{(\rho,\nicefrac{{\beta}}{{r}})}\right)}$
(146) $\displaystyle\leq
e^{\frac{2\beta\lambda}{r}}e^{\frac{2\beta\norm{H}}{r}}\left(\frac{\norm{\mathcal{U}_{B}(\rho,\nicefrac{{\beta}}{{r}})-\mathcal{Q}_{B}{(\rho,\nicefrac{{\beta}}{{r}})}}_{1\rightarrow
1}+\norm{\mathcal{U}_{A}(\rho,\nicefrac{{\beta}}{{r}})-\mathcal{T}_{A}^{2k}{(\rho,\nicefrac{{\beta}}{{r}})}}_{1\rightarrow
1}}{\Tr\left(\mathcal{Q}_{B}{(\rho,\nicefrac{{\beta}}{{r}})}\circ\mathcal{T}_{A}^{2k}{(\rho,\nicefrac{{\beta}}{{r}})}\right)}+\frac{\norm{E_{A,B}(\nicefrac{{\beta}}{{r}})}_{1\rightarrow
1}}{\Tr\left(\mathcal{Q}_{B}{(\rho,\nicefrac{{\beta}}{{r}})}\circ\mathcal{T}_{A}^{2k}{(\rho,\nicefrac{{\beta}}{{r}})}\right)}\right),$
(147)

where in the last line we once again apply the triangle inequality, and also
observe that $\lambda_{A}+\lambda_{B}=\lambda$. Now we only have 2 unfamiliar
terms left to bound, the partitioning error $E_{A,B}(\nicefrac{{\beta}}{{r}})$
and the trace of the composition of the QDrift and Trotter channels. To
analyze the trace, we expand both of the channels into their operator form:

$\displaystyle\Tr\left(\mathcal{Q}_{B}{(\rho,\nicefrac{{\beta}}{{r}})}\circ\mathcal{T}_{A}^{2k}{(\rho,\nicefrac{{\beta}}{{r}})}\right)$
$\displaystyle=\Tr\left(\sum_{\mathbf{j}}P_{\mathbf{j}}\mathbf{V_{j}}S_{2k}\rho
S_{2k}\mathbf{V_{j}}^{\dagger}\right)$ (148)
$\displaystyle=\sum_{\mathbf{j}}P_{\mathbf{j}}\Tr\left(\mathbf{V_{j}}S_{2k}\rho
S_{2k}\mathbf{V_{j}}^{\dagger}\right)$ (149)
$\displaystyle=\sum_{\mathbf{j}}P_{\mathbf{j}}\Tr\left(S_{2k}\mathbf{V_{j}}^{\dagger}\mathbf{V_{j}}S_{2k}\rho\right)$
(150)
$\displaystyle\geq\sum_{\mathbf{j}}P_{\mathbf{j}}\sigma_{min}\left(S_{2k}\mathbf{V_{j}}^{\dagger}\mathbf{V_{j}}S_{2k}\right)\Tr\rho$
(151)
$\displaystyle\geq\sigma_{min}^{2}\left(S_{2k}\right)\sum_{\mathbf{j}}P_{\mathbf{j}}\sigma_{min}^{2}\left(\mathbf{V_{j}}\right)$
(152) $\displaystyle\geq
e^{-\frac{2\beta\lambda_{A}}{r}}\sum_{\mathbf{j}}P_{\mathbf{j}}\sigma_{min}^{2}\left(\mathbf{V_{j}}\right)$
(153)

Now recall the definition of the product formula from the multi-sample QDrift
channel defined in 3 $\mathbf{V_{j}}=\prod^{N}_{k=1}e^{-iH_{j_{k}}\tau}$ with
$\tau=\nicefrac{{\beta\lambda}}{{N}}$. To reiterate, the index $\mathbf{j}$ is
a vector of length $N$, and the index set contains every possible length-$N$
permutation of the $L_{B}$ Hamiltonian terms in the QDrift channel. Since we
already have exponential terms in our bound, we can use a naive lower bound
here, which coincides with the most likely case of drawing the largest term
$N$ times. This gives:

$\displaystyle\geq
e^{-\frac{2\beta\lambda_{A}}{r}}\sum_{\mathbf{j}}P_{\mathbf{j}}\prod^{N}_{k=1}\sigma_{min}^{2}\left(e^{-B_{j_{k}}\frac{\beta\lambda_{B}}{Nr}}\right)$
(154) $\displaystyle\geq
e^{-\frac{2\beta\lambda_{A}}{r}}e^{-\frac{2\norm{B}\beta\lambda_{B}}{r}}$
(155) $\displaystyle=e^{-\frac{2\beta\lambda}{r}},$ (156)

where in the last line $\norm{B}=1$ based by definition of QDrift, and
$\lambda_{A}+\lambda_{B}=\lambda$.

Now combining all of our results we have the following:

$\left(\norm{\mathcal{U}_{B}(\rho,\nicefrac{{\beta}}{{r}})-\mathcal{Q}_{B}{(\rho,\nicefrac{{\beta}}{{r}})}}_{1\rightarrow
1}+\norm{\mathcal{U}_{A}(\rho,\nicefrac{{\beta}}{{r}})-\mathcal{T}_{A}^{2k}{(\rho,\nicefrac{{\beta}}{{r}})}}_{1\rightarrow
1}+\norm{E_{A,B}(\nicefrac{{\beta}}{{r}})}_{1\rightarrow
1}\right)e^{\frac{2\beta\norm{H}}{r}}\left(1+e^{\frac{4\beta\lambda}{r}}\right)$
(157)

This term is a sum of the errors of each respective algorithm, with the
addition of the partitioning error, multiplied by some exponential terms that
come from lower bounding traces, which have been common throughout. Therefore,
this is exactly what we might expect a priori. Going forward we wish to
represent this error for the entire simulation rather than just a single time
step, and then show that the exponential terms become irrelevant in the limit
of interest. For the full time simulation we have

$\norm{\left(\frac{\mathcal{U}{(\rho,\nicefrac{{\beta}}{{r}})}}{\Tr\mathcal{U}{(\rho,\nicefrac{{\beta}}{{r}})}}\right)^{\circ
r}-\left(\frac{\mathcal{X}^{2k}(\rho,\nicefrac{{\beta}}{{r}})}{\Tr\mathcal{X}^{2k}(\rho,\nicefrac{{\beta}}{{r}})}\right)^{\circ
r}}_{1\rightarrow 1}\leq
r\norm{\frac{\mathcal{U}{(\rho,\nicefrac{{\beta}}{{r}})}}{\Tr\mathcal{U}{(\rho,\nicefrac{{\beta}}{{r}})}}-\frac{\mathcal{X}^{2k}(\rho,\nicefrac{{\beta}}{{r}})}{\Tr\mathcal{X}^{2k}(\rho,\nicefrac{{\beta}}{{r}})}}_{1\rightarrow
1},$ (158)

where the quantity on the right is that which has been bounded. Inserting the
imaginary time QDrift bound from Theorem III.1, and the single step Trotter
bound from Equation 132 we have:

$\left(\frac{C\beta^{2}\lambda_{B}^{2}}{N_{B}r}+2\Upsilon^{2k+1}\frac{\alpha_{comm}(A,2k)}{(2k+1)!}\frac{\beta^{2k+1}}{r^{2k}}e^{4\nicefrac{{\beta}}{{r}}\Upsilon\lambda_{A}}\left(e^{4\nicefrac{{\beta}}{{r}}\lambda_{A}}e^{2\nicefrac{{\beta}}{{r}}\norm{A}}+e^{2\nicefrac{{\beta}}{{r}}\norm{A}}\right)+r\norm{E_{A,B}(\nicefrac{{\beta}}{{r}})}_{1\rightarrow
1}\right)e^{\frac{2\beta\norm{H}}{r}}\left(1+e^{\frac{4\beta\lambda}{r}}\right).$
(159)

Next we deal with the partitioning error
$\norm{E_{A,B}(\nicefrac{{\beta}}{{r}})}_{1\rightarrow
1}=\norm{\mathcal{U}{(\rho,\nicefrac{{\beta}}{{r}})}-\mathcal{U}_{B}(\rho,\nicefrac{{\beta}}{{r}})\circ\mathcal{U}_{A}(\rho,\nicefrac{{\beta}}{{r}})}_{1\rightarrow
1}$. We know from Theorem III.2 that
$\norm{\mathcal{U}{(\beta,\rho)}-\mathcal{U}_{B}(\rho,\nicefrac{{\beta}}{{r}})\circ\mathcal{U}_{A}(\rho,\nicefrac{{\beta}}{{r}})}_{1\rightarrow
1}\leq
2e^{\frac{\beta\lambda}{r}}\norm{\mathcal{U}{(\beta,\rho)}-\mathcal{U}_{B}(\rho,\nicefrac{{\beta}}{{r}})\circ\mathcal{U}_{A}(\rho,\nicefrac{{\beta}}{{r}})}_{\infty}$,
then once again using the bound from Appendix E of [10], we have:

$\left(\frac{C\beta^{2}\lambda_{B}^{2}}{N_{B}r}+2\Upsilon^{2k+1}\frac{\alpha_{comm}(A,2k)}{(2k+1)!}\frac{\beta^{2k+1}}{r^{2k}}e^{4\nicefrac{{\beta}}{{r}}\Upsilon\lambda_{A}}\left(e^{4\nicefrac{{\beta}}{{r}}\lambda_{A}}e^{2\nicefrac{{\beta}}{{r}}\norm{A}}+e^{2\nicefrac{{\beta}}{{r}}\norm{A}}\right)+\norm{[A,B]}e^{\frac{\beta\lambda}{r}}\frac{2\beta^{2}}{r}e^{\frac{4\beta\lambda}{r}}\right)e^{\frac{2\beta\norm{H}}{r}}\left(1+e^{\frac{4\beta\lambda}{r}}\right)$
(160)

We are now left with an equation non-linear in $r$, which makes it more
difficult to make asymptotic arguments similar to those in made Theorem III.2.
Since, we wish to eliminate the exponentials from this bound, we first upper-
bound all exponentials by $e^{\nicefrac{{n\Upsilon\beta\lambda}}{{r}}}$ where
we keep the coefficients $n$ for each of them. This allows us to collect the
exponentials and simplify the expression to the following:

$\left(\frac{C\beta^{2}\lambda_{B}^{2}}{N_{B}r}+2\Upsilon^{2k+1}\frac{\alpha_{comm}(A,2k)}{(2k+1)!}\frac{\beta^{2k+1}}{r^{2k}}\left(e^{10\nicefrac{{\Upsilon\beta\lambda}}{{r}}}+e^{6\nicefrac{{\Upsilon\beta\lambda}}{{r}}}\right)+\norm{[A,B]}\frac{2\beta^{2}}{r}e^{\nicefrac{{5\Upsilon\beta\lambda}}{{r}}}\right)\left(e^{\nicefrac{{2\Upsilon\beta\lambda}}{{r}}}+e^{\nicefrac{{6\Upsilon\beta\lambda}}{{r}}}\right)$
(161)

Now, we introduce the constraint that $r\geq\Upsilon\beta\lambda$. Imposing
this constraint over $r$ allows all of the exponentials to be upper-bounded by
a constant such that in the limit that all other parameters approach infinity,
this bound becomes the following:

$\norm{\frac{\mathcal{U}{(\rho,\beta)}}{\Tr\mathcal{U}{(\rho,\beta)}}-\frac{\mathcal{X}^{2k}(\rho,\beta)}{\Tr\mathcal{X}^{2k}(\rho,\beta)}}_{1\rightarrow
1}\in\mathcal{O}\left(\frac{\beta^{2}\lambda_{B}^{2}}{N_{B}r}+\Upsilon^{2k+1}\frac{\alpha_{comm}(A,2k)}{(2k+1)!}\frac{\beta^{2k+1}}{r^{2k}}+\norm{[A,B]}\frac{\beta^{2}}{r}\right),$
(162)

thus completing the proof. ∎

##  Appendix B Numerical Error Analysis

###  B.1 Types of Error Measures

Here, we propose two different error measures and discuss the advantages
associated with each.

####  B.1.1 Infidelity of Quantum States

The fidelity is a quantity that serves as a measure of the closeness between
two quantum states, of which can be either pure or mixed states and is
summarized in the following definition.

######  Definition 10 (Fidelity).

Let $\rho,\sigma\in\mathbb{C}^{2^{n}\times 2^{n}}$ be density matrices that
represent quantum states with equal Hilbert space dimension. The fidelity
between the states is then defined as:

Where $\mathscr{F}(\rho,\sigma)\in[0,1]$ and
$\mathscr{F}(\rho,\sigma)=\mathscr{F}(\sigma,\rho)$.  The infidelity is then
simply defined as $1-\mathscr{F}$ and represents a measure of distance between
two states with the latter properties in the definition still holding. The
computational advantage of using infidelity over some other measures is that
it is relatively cheap to compute, especially for large systems. This does not
appear obvious on first glance given the need to compute 2 matrix square roots
which often involves diagonalization. However, we can take advantage of a nice
property: $\mathscr{F}=\Tr(\rho,\sigma)$ if one of $\rho$ or $\sigma$ is pure.
A simple proof is provided in the following lemma.

######  Lemma 11 (Pure-Non-Pure Fidelity Property).

Given a two density matrices $\rho$ and $\sigma$ where at least one of which
is pure, then the expression for the fidelity simplifies to
$\mathscr{F}(\rho,\sigma)=\Tr(\rho\sigma)$.

###### Proof.

$\displaystyle\mathscr{F}(\rho,\sigma):=$
$\displaystyle\Tr\left[\sqrt{\sqrt{\rho}\sigma\sqrt{\rho}}\right]^{2}$ (164)
$\displaystyle=\Tr\left[\sqrt{\sqrt{\outerproduct{\psi}{\psi}}\sigma\sqrt{\outerproduct{\psi}{\psi}}}\right]^{2}$
(165)
$\displaystyle=\Tr\left[\sqrt{\outerproduct{\psi}{\psi}\sigma\outerproduct{\psi}{\psi}}\right]^{2}$
(166)
$\displaystyle=\langle\psi|\sigma|\psi\rangle\Tr(\outerproduct{\psi}{\psi})$
(167) $\displaystyle=\Tr(\rho\sigma)$ (168)

Here we choose $\rho$ to be a pure state $\outerproduct{\psi}{\psi}$ which is
a projector and thus $\rho=\rho^{2}=\sqrt{\rho}$. We use these properties
along with the fact that the trace of a pure state is equal to one to complete
the proof. Note, $\mathscr{F}$ is also known to be symmetric
$\mathscr{F}(\rho,\sigma)=\mathscr{F}(\sigma,\rho)$ but we do not show this
here. This implies the statement is true for either $\rho$ or $\sigma$ pure. ∎

The important part of this lemma is that we only require one of the density
matrices to be a pure state. To evaluate the error of a composite channel, we
will be taking the infidelity between channel outputs, and due to the nature
of QDrift the channel will always output a mixed state. However, the exact
evolution channel providing the means for comparison only produces pure states
given pure state input. With these definitions and properties we can then
present the following strategy: Reformulate the fidelity expression such that
it can be written as an expectation value that only requires matrix vector
multiplication to compute, and then use Monte-Carlo sampling to approximate
this quantity. This is shown via the following derivation:

$\displaystyle\mathscr{F}(\mathcal{U}(\rho,t),\mathcal{X}^{2k}(\rho,t))$
$\displaystyle=\Tr(e^{-iHt}\outerproduct{\psi}{\psi}e^{iHt}e^{-iAt}e^{-iBt}\outerproduct{\psi}{\psi}e^{iBt}e^{iAt})$
(169)
$\displaystyle=\Tr(e^{-iHt}\outerproduct{\psi}{\psi}e^{iHt}{S_{2k}{\left(-iAt\right)}}\sum_{j}P_{j}e^{-iB_{j}\tau}\outerproduct{\psi}{\psi}e^{iB_{j}t}{S_{2k}{\left(-iAt\right)}}^{\dagger})$
(170)
$\displaystyle=\sum_{j}P_{j}\Tr(\bra{\psi}e^{iB_{j}t}{S_{2k}{\left(-iAt\right)}}^{\dagger}e^{-iHt}\outerproduct{\psi}{\psi}e^{-iHt}{S_{2k}{\left(-iAt\right)}}e^{-iB_{j}\tau}\ket{\psi})$
(171)
$\displaystyle=\sum_{j}P_{j}|\bra{\psi}e^{iB_{j}\tau}{S_{2k}{\left(-iAt\right)}}^{\dagger}e^{-iHt}\ket{\psi}|^{2}$
(172)
$\displaystyle=\mathbb{E}\left[|\bra{\psi}e^{iB_{j}\tau}{S_{2k}{\left(-iAt\right)}}^{\dagger}e^{-iHt}\ket{\psi}|^{2}\right]$
(173)
$\displaystyle\approx\frac{1}{N}\sum_{j}^{N}|\bra{\psi}e^{iB_{j}\tau}{S_{2k}{\left(-iAt\right)}}^{\dagger}e^{-iHt}\ket{\psi}|^{2}$
(174)

Now when numerically evaluating the error, only matrix-vector multiplication
is needed as opposed to matrix-matrix before. Here we no longer need to build
exact density matrices outputted by the composite channel, and the step
involving the matrix square root has been eliminated. However, the trade-off
is that we have to use Monte-Carlo sampling to estimate this expectation
value. Given that the naive methods for matrix-vector and matrix-matrix
multiplication scale like $\mathcal{O}\left(n^{2}\right)$ and
$\mathcal{O}\left(n^{3}\right)$ respectively, for small systems the overhead
of Monte-Carlo sampling may lead to a loss of any computational advantage.
However, for large systems where $n^{3}\gg n^{2}$, this sampling method is
quite advantageous. It is also important to note that sampling makes QDrift
channel much cheaper to study as well. For a QDrift channel with $L$ terms and
$N$ samples, building the exact channel yields additional matrix
multiplication overhead $\mathcal{O}\left(NL\right)$ (not considering
arithmetic overhead). While this is does not seem prohibitive, it does become
expensive when we have a Hamiltonian with many terms that we wish to sample
many times for a high degree of accuracy. It should be noted, however, that
this is not as bad as it initially seems given that there are $L^{N}$ terms
QDrift channel sum. The numerical routine that computes this sum was optimized
to achieve this, by using the fact that channel samples are i.i.d. With
intermediate summation to produce a new input $\rho$ after each sample vector,
we do not need to compute each term in this sum individually:

$\sum_{j}p_{j}e^{-iH_{j}t}\rho_{i}e^{iH_{j}t}\rightarrow\rho_{i+1}$ (175)

repeating this procedure $N$ times until $\rho_{i}\rightarrow\rho_{N}$ which
constructs the exact QDrift channel output with $\mathcal{O}\left(NL\right)$
matrix multiplications. We mention this to highlight the advantages of Monte-
Carlo sampling where it can be done reliably, as well as to highlight how
QDrift is numerically computed exactly in our work (which will be made
necessary by results in the following section).

####  B.1.2 Trace Distance over Density Matrices

The next error measure we wish to investigate is the trace distance, which is
equivalent to the Schatten 1-norm. The trace distance is a metric over the
space of density matrices that can be defined in the following way:

######  Definition 12 (Trace Distance).

Let $\rho,\sigma\in\mathbb{C}^{2^{n}\times 2^{n}}$ be density matrices that
represent quantum states with equal Hilbert space dimension. The trace
distance between the states is then defined as:

Where $\mathscr{T}(\rho,\sigma)\in[0,1]$.

Given that $\rho$ and $\sigma$ are necessarily hermitian, the trace distance
simplifies to
$\mathscr{T}(\rho,\sigma)=\sum_{i}\absolutevalue{\lambda_{i}}^{2}$ where
$\lambda_{i}$ are eigenvalues of the matrix $(\rho-\sigma)$. There are no
simplifications in our approach to computing this numerically. Here, we must
compute the exact density matrix output of the composite channel (a
necessarily mixed state outputted by the QDrift channel component) and the
output of the unitary channel. Exact diagonalization is performed to obtain
the eigenvalues, and the sum above is then computed. Providing a rigorous
argument for crossover points in the compute time between the fidelity and
trace distance (as a function of $N,L$) is beyond the scope of this paper,
however, given the arguments from the previous section our intuition is that
for small systems the trace distance will be computable in a reasonable amount
of time whereas Monte-Carlo sampling the infidelity will become more
advantageous as the system size grows.

###  B.2 Comparison of Error Measures

In the following section we provide short proofs regarding how the error of
each of the QDrift and Trotter-Suzuki algorithms scale with time, both in the
infidelity and trace distance framework. It is important to note that we are
specifically using the Schatten 1-norm over density matrices and not an
induced 1-norm over channels. Therefore, we are computing this quantity using
the outputs of the two channels, and the same goes for the infidelity which is
only defined on the density matrices that are the outputs of the channels in
question. We begin with the infidelity:

See V.1

###### Proof.

$\displaystyle
1-\mathscr{F}(\mathcal{Q}{(\rho,t)},\mathcal{U}{(\rho,\nicefrac{{t}}{{N}})}):=$
$\displaystyle
1-\Tr(\mathcal{Q}{(\rho,t)}\mathcal{U}(\rho,\nicefrac{{t}}{{N}}))$ (177)
$\displaystyle=$ $\displaystyle
1-\Tr(\sum_{j}P_{j}e^{-iH_{j}\tau}\outerproduct{\psi}{\psi}e^{iH_{j}\tau}e^{-iHt/N}\outerproduct{\psi}{\psi}e^{iHt/N})$
(178) $\displaystyle=$ $\displaystyle
1-\sum_{j}P_{j}\bra{\psi}e^{iHt/N}e^{-iH_{j}\tau}\outerproduct{\psi}{\psi}e^{iH_{j}\tau}e^{-iHt/N}\ket{\psi}$
(179) $\displaystyle=$ $\displaystyle 1-\sum_{j}P_{j}|\bra{\psi}(\openone-
iHt/N-\frac{H^{2}t^{2}/N^{2}}{2!}+...)(\openone+iH_{j}\tau-\frac{H_{j}^{2}\tau^{2}}{2!}+...)\ket{\psi}|^{2}$
(180) $\displaystyle=$ $\displaystyle 1-\sum_{j}P_{j}|\bra{\psi}(\openone-
iHt/N+iH_{j}\tau+HH_{j}t\tau/N-\frac{H^{2}t^{2}/N^{2}}{2!}-\frac{H_{j}^{2}\tau^{2}}{2!}+...)\ket{\psi}|^{2}$
(181) $\displaystyle=$ $\displaystyle 1-\sum_{j}P_{j}|\bra{\psi}(\openone-
iHt/N+iH_{j}t\lambda/N+HH_{j}t^{2}\lambda/N^{2}$ (182)
$\displaystyle-\frac{H^{2}t^{2}/N^{2}}{2!}-\frac{H_{j}^{2}t^{2}\lambda^{2}/N^{2}}{2!}+...)\ket{\psi}|^{2}$
(183) $\displaystyle=$ $\displaystyle
1-\sum_{j}P_{j}(1-2\langle\psi|H|\psi\rangle\langle\psi|H_{j}|\psi\rangle
t^{2}\lambda/N^{2}+\langle\psi|H|\psi\rangle^{2}t^{2}/N^{2}+\langle\psi|H_{j}|\psi\rangle^{2}t^{2}\lambda^{2}/N^{2}+ct^{2}...)$
(184) $\displaystyle=$ $\displaystyle
t^{2}\sum_{j}P_{j}(2\langle\psi|H|\psi\rangle\langle\psi|H_{j}|\psi\rangle\lambda/N^{2}-\langle\psi|H|\psi\rangle^{2}/N^{2}-\langle\psi|H_{j}|\psi\rangle^{2}\lambda^{2}/N^{2}-c...)$
(185) $\displaystyle\leq$ $\displaystyle
t^{2}\sum_{j}P_{j}(|2\langle\psi|H|\psi\rangle\langle\psi|H_{j}|\psi\rangle|\lambda/N^{2}+|\langle\psi|H|\psi\rangle^{2}|/N^{2}+|\langle\psi|H_{j}|\psi\rangle^{2}|\lambda^{2}/N^{2}+|c|...)\in\mathcal{O}\left(t^{2}\right)$
(186)

We use the cyclic property of the trace and the fact that $\rho$ is a pure
state $\outerproduct{\psi}{\psi}$ and then replace the operator exponentials
with their respective Taylor series expansions, keeping terms only up to
$t^{2}$ throughout. In the final step, the modulus is taken and the triangle
inequality applied to upper bound the sum. This is done for the purpose of
ensuring the signs of the terms are “well behaved” such that
$\mathcal{O}\left(\right)$ notation can be aptly applied. As well, $c$ is used
for brevity where
$c=\langle\psi|(HH_{j}\lambda/N^{2}-\frac{H^{2}}{2!}/N^{2}-\frac{H_{j}^{2}\lambda^{2}/N^{2}}{2!})|\psi\rangle+h.c$.
Here the notation $\mathcal{O}\left(t^{2}\right)$ is understood in the limit
as $t\rightarrow 0$. ∎

Given the diamond distance bound in [11], the $\mathcal{O}\left(t^{2}\right)$
infidelity scaling is not surprising. Not only this, but recall from Equation
180, where the largest order in commutator error from the product formula was
quadratic, and given that QDrift applies no symmetrization strategy, we did
not expect this to disappear. We would like to remind the reader that the
shining feature of QDrift is that the cost is independent of the number of
terms in the Hamiltonian. Next, we examine how a Trotter-Suzuki channel
behaves in the infidelity measure:

See V.2

###### Proof.

$\displaystyle
1-\mathscr{F}(\mathcal{T}^{2k}{(\rho,t)},\mathcal{U}(\rho,t)):=$
$\displaystyle 1-\Tr(\mathcal{T}^{2k}{(\rho,t)}\mathcal{U}{(\rho,t)})$ (187)
$\displaystyle=$ $\displaystyle
1-\Tr({S_{2k}{\left(-iHt\right)}}\outerproduct{\psi}{\psi}{S_{2k}{\left(-iHt\right)}}^{\dagger}e^{-iHt}\outerproduct{\psi}{\psi}e^{iHt})$
(188) $\displaystyle=$ $\displaystyle
1-\bra{\psi}e^{iHt}{S_{2k}{\left(-iHt\right)}}\outerproduct{\psi}{\psi}{S_{2k}{\left(-iHt\right)}}^{\dagger}e^{-iHt}\ket{\psi}$
(189) $\displaystyle=$ $\displaystyle
1-|\bra{\psi}e^{iHt}{S_{2k}{\left(-iHt\right)}}\ket{\psi}|^{2}$ (190)
$\displaystyle=$ $\displaystyle
1-|\bra{\psi}({S_{2k}{\left(-iHt\right)}}^{\dagger}\pm it^{2k+1}R_{2k+1}\mp
it^{2k+3}R_{2k+3}+...){S_{2k}{\left(-iHt\right)}}\ket{\psi}|^{2}$ (191)
$\displaystyle=$ $\displaystyle 1-|\bra{\psi}\openone+\pm
it^{2k+1}R_{2k+1}{S_{2k}{\left(-iHt\right)}}\ket{\psi}\mp
it^{2k+3}R_{2k+3}{S_{2k}{\left(-iHt\right)}}\ket{\psi}+...|^{2}$ (192)
$\displaystyle=$ $\displaystyle 1-|1\pm
it^{2k+1}\langle\psi|R_{2k+1}{S_{2k}{\left(-iHt\right)}}|\psi\rangle\mp
it^{2k+3}\langle\psi|R_{2k+3}{S_{2k}{\left(-iHt\right)}}|\psi\rangle\pm...|^{2}$
(193) $\displaystyle\leq$ $\displaystyle
1-|1|^{2}+(t^{2k+1})^{2}|\langle\psi|R_{2k+1}{S_{2k}{\left(-iHt\right)}}|\psi\rangle|^{2}+...$
(194) $\displaystyle=$
$\displaystyle(t^{2k+1})^{2}|\langle\psi|R_{2k+1}{S_{2k}{\left(-iHt\right)}}|\psi\rangle|^{2}+...\in\mathcal{O}\left((t^{2k+1})^{2}\right)$
(195)

We again use the cyclic property of the trace and the fact that $\rho$ is a
pure state $\outerproduct{\psi}{\psi}$. The evolution operator $e^{-iHt}$ is
rewritten in terms of the product formula and corrections at significant
orders. $R_{n}$ are remainder operators of order $n$ of which we keep to order
$2k+3$. In the second last step, we use the triangle inequality and ignore
terms of higher order. Keeping higher order terms is unnecessary in that we
cannot use them to generate terms of order $\leq(t^{2k+1})^{2}$. Once again,
$\mathcal{O}\left((t^{2k+1})^{2}\right)$ is understood in the limit as
$t\rightarrow 0$. ∎

This theorem provides a peculiar result. The expected error in the output of a
Trotter-Suzuki channel from section III.2 is squared by the nature of the
infidelity calculation. Given that $\epsilon^{2}\leq\epsilon\leq 1$, error
squaring can lead to a much lower gate cost, and given that Trotter
experiences this phenomenon while QDrift does not, a numerical optimizer may
heavily favour Trotter and the potential advantages from the composite
approach will likely fade. This is the main issue with choosing infidelity as
an error measure whilst hybridizing algorithms with different error scaling.
The following two theorems show that the trace distance suffers no such side
effect:

See V.3

###### Proof.

$\displaystyle\mathscr{T}(\mathcal{Q}{(\rho,t)},\mathcal{U}(\rho,\nicefrac{{t}}{{N}})):=$
$\displaystyle\Tr\left[\sqrt{(\mathcal{Q}{(\rho,t)}-\mathcal{U}(\rho,\nicefrac{{t}}{{N}}))(\mathcal{Q}{(\rho,t)}-\mathcal{U}(\rho,t/N))^{\dagger}}\right]$
(196) $\displaystyle=$
$\displaystyle\Tr\left[\sqrt{(\mathcal{Q}{(\rho,t)}-\mathcal{U}(\rho,\nicefrac{{t}}{{N}}))^{2}}\right]$
(197) $\displaystyle=$
$\displaystyle\Tr\left[\sqrt{((\sum_{j}P_{j}e^{-iH_{j}\tau}\rho
e^{iH_{j}\tau})-(e^{-iHt/N}\rho e^{iHt/N}))^{2}}\right]$ (198)
$\displaystyle=$ $\displaystyle\Tr\left[\sqrt{((\sum_{j}P_{j}(\openone-
iH_{j}\tau-H_{j}^{2}\tau^{2}+...)\rho(\openone+iH_{j}\tau-
H_{j}^{2}\tau^{2}+...)-(e^{iHt/N}\rho e^{-iHt/N}))^{2}}\right]$ (199)
$\displaystyle=$ $\displaystyle\Tr\left[\sqrt{(\sum_{j}P_{j}(\rho+i\lambda
t/N[H_{j},\rho]+\lambda^{2}t^{2}/N^{2}(H_{j}\rho
H_{j}-i\\{H_{j},\rho\\})+...)-(e^{iHt/N}\rho e^{-iHt/N}))^{2}}\right]$ (200)
$\displaystyle=$
$\displaystyle\Tr\left[\sqrt{((\rho+\sum_{j}\frac{h_{j}}{\lambda}(i\lambda
t/N[H_{j},\rho]+\lambda^{2}t^{2}/N^{2}(H_{j}\rho
H_{j}-i\\{H_{j},\rho\\})+...)-(e^{iHt/N}\rho e^{-iHt/N}))^{2}}\right]$ (201)
$\displaystyle=$
$\displaystyle\Tr\left[\sqrt{(\rho+it/N[H,\rho]+\sum_{j}\frac{h_{j}}{\lambda}(\lambda^{2}t^{2}/N^{2}(H_{j}\rho
H_{j}-i\\{H_{j}^{2},\rho\\})+...)-(e^{iHt/N}\rho e^{-iHt/N}))^{2}}\right]$
(202) $\displaystyle=$
$\displaystyle\Tr\left[\sqrt{\sum_{j}\frac{h_{j}}{\lambda}(\lambda^{2}t^{2}/N^{2}(H_{j}\rho
H_{j}-i\\{H_{j}^{2},\rho\\}...)+(t^{2}(H\rho
H+i\\{H^{2},\rho\\})+...))^{2}}\right]$ (203) $\displaystyle=$ $\displaystyle
t^{2}\Tr\left[\sqrt{(R(t))^{2}}\right]\in\mathcal{O}\left(t^{2}\right)$ (204)

Here we expand the Qdrift channel with the Taylor series of the exponential
operators and show that the terms of this series cancel with that of the
Taylor series of $\mathcal{U}_{N}$ up to order $t^{2}$. To achieve this we
simply use the definitions of $P_{j}=h_{j}/\lambda$ and $\sum_{j}H_{j}=H$. In
the last line $R(t)$ is considered a remainder operator which is only a
function of polynomials $t^{n},n\geq 0$. The notation
$\mathcal{O}\left(t^{2}\right)$ is understood in the limit as $t\rightarrow
0$. ∎

See V.4

###### Proof.

$\displaystyle\mathscr{T}(\mathcal{T}^{2k}{(\rho,t)},\mathcal{U}(\rho,t)):=\Tr\left[\sqrt{(\mathcal{T}^{2k}{(\rho,t)}-\mathcal{U}(\rho,t))(\mathcal{T}^{2k}{(\rho,t)}-\mathcal{U}(\rho,t))^{\dagger}}\right]$
(205)
$\displaystyle=\Tr\left[\sqrt{(\mathcal{T}^{2k}{(\rho,t)}-\mathcal{U}(\rho,t))^{2}}\right]$
(206)
$\displaystyle=\Tr\left[\sqrt{(({S_{2k}{\left(-iHt\right)}}\rho{S_{2k}{\left(-iHt\right)}}^{\dagger})-(e^{-iHt}\rho
e^{iHt}))^{2}}\right]$ (207)
$\displaystyle=\Tr\left[\sqrt{(({S_{2k}{\left(-iHt\right)}}\rho{S_{2k}{\left(-iHt\right)}}^{\dagger})-(({S_{2k}{\left(-iHt\right)}}^{\dagger}\pm
it^{2k+1}R_{2k+1}+...)\rho({S_{2k}{\left(-iHt\right)}}^{\dagger}\mp
it^{2k+1}R_{2k+1}^{\dagger}+...))^{2}}\right]$ (208)
$\displaystyle=\Tr\left[\sqrt{((\mp it^{2k+1}{S_{2k}{\left(-iHt\right)}}\rho
R^{\dagger}_{2k+1}\pm
it^{2k+1}R_{2k+1}{S_{2k}{\left(-iHt\right)}}^{\dagger}+...))^{2}}\right]$
(209)
$\displaystyle=t^{2k+1}\Tr\left[\sqrt{(R(t))^{2}}\right]\in\mathcal{O}\left(t^{2k+1}\right)$
(210)

Here expand the exponential operators in the evolutionary channel
$\mathcal{U}$ in terms of their error terms for some general $2k$th order
Trotter-Suzuki channel, where $R_{n}$ are remainder operators of order $n$. As
the algebra proceeds we only keep terms of order $\leq t^{2k+1}$ throughout.
Keeping higher order terms is unnecessary in that we cannot use them to
generate terms of order $\leq t^{2k+1}$. In the last line $R(t)$ is considered
a general remainder operator which is only a function of polynomials
$t^{n},n\geq 0$. Once again, $\mathcal{O}\left(t^{2k+1}\right)$ is understood
in the limit as $t\rightarrow 0$. ∎

(a) QDrift Infidelity

(b) Trotter Infidelity

(c) TS Order $2k=2$ Infidelity

(d) QDrift Trace Distance

(e) Trotter Trace Distance

(f) TS order $2k=2$ Trace Distance

Figure 12: Numerical plots of log infidelity and trace distance scaling with
log time for both Trotter-Suzuki and QDrift channels. Plots are all
constructed using 1 iteration or sample respectively. Note that the slope of
the plots follow the analytical predictions exactly, with the infidelity
squaring the expected Trotter error in (b) and (c). Note that in (c) we run
into machine precision error at small times and for long times in (a) and (d)
QDrift no longer converges. The plots were compute using 100 Monte-Carlo
samples for QDrift Infidelity, and exact density matrices elsewhere for a
Graph Hamiltonian with 5 spins.

There are no surprises with the trace distance error as it yields the expected
time scaling from the results of section B.1.2. Numerical investigations also
support our proofs (see Figure 12). Therefore, while it may be more expensive
to compute for larger systems, this makes it favourable to work with as it
provides a more “fair” measure, treating the algorithms on more equal footing.
Thus, it is expected that the composite channel will appear more performant
using this metric. In addition, using this framework requires no sampling,
which means that monotonicity of the cost will be guaranteed with respect to
channel iterations. This will prove useful in the cases in which we choose to
optimize over possible partitions with respect to the cost. Now, out of
interest, we also wish to numerically investigate the trace distance of the
composite channel. Analytically, this is likely a messy problem given the
partition, but intuitively we expect the scaling to be some linear combination
$\alpha t^{2}+\beta t^{3}+\gamma$ with the slope of the log plot being between
2 and 3 (see Figure 13).

Figure 13: Trace distance time scaling for a composite channel in inner order
$2k=2$. This plot emits an interesting structure in that the time dependence
of the trace distance goes through a “phase transition” between QDrift (red)
and 2nd order Trotter-Suzuki (blue) dominated regions. Given that QDrift
excels in short time simulations and whereas the higher order Trotter-Suzuki
dominates at longer simulation times this is not entirely unexpected.


# CrossStack: A 3-D Reconfigurable RRAM Crossbar Inference Engine

Jason K. Eshraghian1, Kyoungrok Cho2 and Sung Mo Kang3
1School of Electrical, Electronic and Computer Engineering, University of
Michigan, Ann Arbor, MI 48109 USA 2College of Electrical and Computer
Engineering, Chungbuk National University, Cheongju 362763, South Korea 3Jack
Baskin School of Engineering, University of California, Santa Cruz, Santa
Cruz, CA 95064 USA

###### Abstract

Deep neural network inference accelerators are rapidly growing in importance
as we turn to massively parallelized processing beyond GPUs and ASICs. The
dominant operation in feedforward inference is the multiply-and-accumlate
process, where each column in a crossbar generates the current response of a
single neuron. As a result, memristor crossbar arrays parallelize inference
and image processing tasks very efficiently. In this brief, we present a 3-D
active memristor crossbar array ‘CrossStack’, which adopts stacked pairs of
Al/TiO2/TiO_2-x/Al devices with common middle electrodes. By designing CMOS-
memristor hybrid cells used in the layout of the array, CrossStack can operate
in one of two user-configurable modes as a reconfigurable inference engine: 1)
expansion mode and 2) deep-net mode. In expansion mode, the resolution of the
network is doubled by increasing the number of inputs for a given chip area,
reducing IR drop by 22%. In deep-net mode, inference speed per-10-bit
convolution is improved by 29% by simultaneously using one TiO2/TiO_2-x layer
for read processes, and the other for write processes. We experimentally
verify both modes on our $10\times 10\times 2$ array.

###### Index Terms:

deep learning, in-memory computing, memristors, neural network, RRAM.

##  I Introduction

Increasing the sizes of artificial neural networks (ANNs) has been the most
common response to the copious amount of data being continuously generated.
Where training sets are in excess of billions of inputs processed through
hundreds of millions of parameters in a neural network [1], new ways to speed
up the processing of all this information must be developed. Since 2012, the
training runtime of neural networks has doubled every 3–4 months. It is
equally important to develop hardware that is not only optimized for running
very large scale networks, but is adaptable to the unceasing wave of emerging
ANN topologies.

Memristors are now ubiquitous in neuromorphic computing literature due to
their long retention [2, 3, 4], excellent scalability [5, 6], fast read and
write speeds [7, 8], compatibility with CMOS technology [9, 10, 11, 12], and
precise weight updates [13, 14]. The development of dense integrated
structures with 3-D stacked crossbar arrays enables an increase in the
throughput for a given chip area, but thus far, target applications of 3-D
RRAM have mostly been limited to digital memory [16, 17, 18, 19, 20, 21, 22,
23].

It can be difficult for ASIC designs to keep pace with the latest developments
in machine learning due to the lag time between algorithm development and the
full IC design cycle. With popular machine learning methods in a rapidly
evolving state, reconfigurability is of paramount importance to ensure
hardware does not become obsolete the moment new network architectures and
topologies are introduced. In response to this, we present a reconfigurable
stacked pair of memristor crossbars dubbed ‘CrossStack’ that can be operated
in one of two modes: 1) expansion mode, and 2) deep-net mode. In expansion
mode, the resolution of the network is doubled by increasing the number of
inputs for a given chip area, thus reducing IR drop by 22% of an equivalent
array. In deep-net mode, inference speed per-10-bit convolution is improved by
29% by simultaneously using one array for read processes, and the other for
write processes. This brief will demonstrate how to selectively isolate and
couple the two layers using CMOS cell design. We experimentally verify this on
our in-house fabricated crossbar stack, using separately controlled CMOS
circuitry in the SK Hynix 180nm process.

##  II Matrix-Vector Multiplication

To perform analog-domain multiply-and-accumulate (MAC) using RRAM arrays, a
voltage vector $\boldsymbol{V}^{i}$ is applied at the input, multiplied by a
conductance matrix $\boldsymbol{G}$, to generate a current vector of
$\boldsymbol{i}=\boldsymbol{V}^{i}\boldsymbol{G}$ in accordance with Ohm’s Law
and Kirchhoff’s Current Law. On a pre-trained network, the conductance of each
memristor is programmed prior to read-out. For further detail on MAC on a
crossbar, we recommend referring to [26, 27]. In a conventional crossbar, the
memristors must be programmed prior to read-out. Where the number of
parameters exceed the memory resources available, RRAM cells must be
reprogrammed while the data flow of the activations is stalled. CrossStack
avoids possible stalling that may occur by adding the option to pipeline the
read and write processes simultaneously across the two layers. How this is
achieved will be described in the following section.

##  III Circuit Operation Modes

(a)

(b)

(c)

(d)

(e)

(f)

Figure 1: (a) Simplified 3-D memristor crossbar array in expansion mode with
cumulative current through the shared electrode (b) Memristor-CMOS cell (c)
Current flow in read mode when read-enable is set high (d) Current flow in
write mode when read-enable is set low (e) A stacked pair of cells during the
read cycle in expansion mode (f) A stacked pair of cells in deep-net mode.

A simplified structure of CrossStack is depicted in Fig. 1(a), and the
memristor-CMOS cell schematic is given in Fig. 1(b). This work presents two
modes in which CrossStack may operate in. These modes are controlled by an
active-high read-enable signal RE.

###  III-A Expansion Mode

Expansion mode enables access to a shared column line from memristors both
above and below the wire. Vertical stacking of memristors doubles the number
of possible inputs and weights to each neuron for a given length of column
wire which is illustrated in Fig. 1(a) when compared to conventional
crossbars. This can be formalized by the following equation:

$\begin{bmatrix}i_{1}&i_{2}&...&i_{m}\end{bmatrix}=\begin{bmatrix}V_{1}^{i}\\\
V_{2}^{i}\\\ \vdots\\\
V_{n}^{i}\end{bmatrix}^{T}\begin{bmatrix}G_{1,1}&G_{1,2}&...&G_{1,m}\\\
G_{2,1}&G_{2,2}&...&G_{2,m}\\\ \vdots&\vdots&\ddots&\vdots\\\
G_{n,1}&G_{n,2}&...&G_{n,m}\end{bmatrix}$ (1)

where $m$ is the number of columns in the crossbar and $n$ is the number of
rows. In expansion mode, $n$ is double that of a 2-D array, as there are rows
both above and below the shared column contributing to output current.

To activate a cell in expansion mode, the read-enable signal of all cells must
be identical. The current pathway from input to output of a single cell is
depicted in Fig. 1(c), and a pair of stacked cells is shown in Fig. 1(e). To
generate a read-out current at each shared column, RE must be set high (in our
case, $V\geq V_{Th}=0.4V$; described in further detail in our experimental
results). If transistors N1 and N2 are treated as switches, then N2 would be
off and N1 would be on. Therefore, a current pathway is formed from both upper
and lower crossbar arrays to the column line. To program a memristor, RE is
set low as in Fig. 1(d). This causes N2 to switch on and N1 off, which forms a
pathway from the memristor to ground and prevents current from flowing to the
shared column. Therefore, the two crossbars can be programmed independently of
one another by isolating them. The transistors are sized to ensure a
negligible leakage current, and to sustain a sufficiently low ON resistance in
comparison to the memristance while it is operating in the linear region such
that it behaves as an ideal switch.

###  III-B Deep-net Mode

Deep-net mode ensures both pairs of arrays are isolated from one another at
all times. Isolation biasing enables each layer to operate independently. The
two arrays must have complementary RE signals as distinct from expansion mode,
depicted in Fig. 1(f). This means that while one array generates a read-out
current, the other array of memristors are being programmed (in write mode).
As described above, when RE is low the cell does not contribute to read-out
current and input voltage $V^{i}$ for the write layer is applied such that the
conductances written to the memristors correspond to the weights of the next
hidden layer in the neural network. Once the analog output current is
digitized as a voltage, the write-layer has been pre-programmed and there is
no need to buffer the current or store it in memory, as is required by most
other pipelines [27, 29, 30]. This process is repeated, but now the roles of
the crossbars are reversed. The original read-layer is now programmed to the
next hidden layer of weights, and the original write-layer generates the read-
out current. Thus, read-write processes run in parallel. One layer is
programmed in anticipation of the output from the other layer, which enables a
novel in-situ pipeline.

In describing how to program a memristor, the key difference between expansion
and deep-net modes is that in expansion, all cells are identically biased for
either read or write at any given time. In deep-net mode, 50% of the cells are
biased for read processes and the remaining 50% are biased for write
processes, only switching once each operation is complete. The shorter read-
out time is subsumed within the programming cycle, but at the expense of half
the number of inputs $n$ in (1). We quantitatively demonstrate this in our
experimental results.

##  IV Experimental Results

###  IV-A Crossbar Fabrication

Figure 2: A $10\times 10\times 2$ prototype of CrossStack, with a cross-
sectional view of the active layer taken using a focus ion beam analyzer.

CrossStack was constructed based on two monolithically integrated crossbar
arrays with a shared central electrode which make up the column line, based on
a sandwich structure of Al/TiO2/TiO_2-x/Al layers. A layer of Al (200-nm-thick
and 20-$\mu$m-wide) was deposited using photolithography on a glass wafer as
the bottom electrode (irradiated using mask alignment for 100 s, subsequently
developed at 23∘C for 120 s). Any excess Al outside of the channel region was
removed via wet etching (H3PO4 : HNO3 : CH3COOH : H2O = 80 ml : 5 ml : 5 ml :
10 ml) at a rate of $\Delta$d/t = 300 nm/min. TiO2 (5-nm-thick) and TiO_2-x
(15-nm-thick) thin films were formed by atomic layer deposition and magnetron
sputtering to fabricate the memristor. Another 200-nm-thick layer of Al was
sputtered as the top electrode using photolithography to create a 20 $\mu$m
$\times$ 20 $\mu$m mask. After a planarization step, the top stack of active
layer and metal were also deposited. Note that the polarity of the pair of
active layers are mirrored, as distinguishable from [31]. This allows for
identical input voltages to be applied when programming the memristors. Fig. 2
shows a working $10\times 10\times 2$ prototype of CrossStack, and a cross-
sectional view of the memristor taken using a focus ion beam analyzer.

TABLE I: CrossStack Characteristics |  |
---|---|---
[HTML]FFFFFFSymbol |  [HTML]FFFFFFParameter |  [HTML]FFFFFFValue
[HTML]FFFFFF$R_{s}$ |  [HTML]FFFFFFstatic resistance of set; $R_{s}=\frac{1}{g_{s}}$ |  [HTML]FFFFFF10K$\Omega~{}\pm$ 7%
[HTML]FFFFFF$R_{r}$ |  [HTML]FFFFFFstatic resistance of reset;$R_{r}=\frac{1}{g_{r}}$ |  [HTML]FFFFFF100K$\Omega~{}\pm$ 10%
[HTML]FFFFFF$V_{DD}$ |  [HTML]FFFFFFsupply voltage |  [HTML]FFFFFF1.8 V
[HTML]FFFFFF$V_{read}$ |  [HTML]FFFFFFread voltage |  [HTML]FFFFFF0.5 V
[HTML]FFFFFF$V_{write}$ |  [HTML]FFFFFFwrite voltage |  [HTML]FFFFFF1.2 V
[HTML]FFFFFF$t_{read}$ |  [HTML]FFFFFFcurrent-read out time |  [HTML]FFFFFF10 ns
[HTML]FFFFFF$t_{write}$ |  [HTML]FFFFFFprogramming time |  [HTML]FFFFFF250 ns
[HTML]FFFFFF$n$ |  [HTML]FFFFFFnumber of memristors |  [HTML]FFFFFF200
[HTML]FFFFFF$v_{Th}$ |  [HTML]FFFFFFthreshold voltage |  [HTML]FFFFFF$|$0.4$|$ V
[HTML]FFFFFF$P_{critical}$ |  [HTML]FFFFFFworst case power consumption |  [HTML]FFFFFF2.9 mW
[HTML]FFFFFF$R_{wire}$ |  [HTML]FFFFFFwire resistance |  [HTML]FFFFFF3.2$\Omega$ p/cell
[HTML]FFFFFF$A_{cell}$ |  [HTML]FFFFFFcell area |  [HTML]FFFFFF20 $\mu$m $\times$ 20$\mu$m
[HTML]FFFFFF$W/L$ |  [HTML]FFFFFFtransistor sizing |  [HTML]FFFFFF450nm/180nm = 2.5
|  |

###  IV-B Cell Test

(a)

(b)

(c)

(d)

Figure 3: Experimental results (a) Pinched hysteresis loop of a single
memristor (b) IR loss comparison in expansion mode (c) Worst-case leakage
current through transistor N1 during a write cycle in deep-net mode with Monte
Carlo parameter sweep of memristance (d) Extreme test case under large input
write signal resulting in a nonlinear voltage drop across the memristor in
deep-net mode.

The CMOS cell was designed in the SK Hynix 180-nm process where $v_{DD}$ = 1.8
V, $v_{Th}$ = 0.4 V, and our parameters are summarized in Table I. We used a
read voltage in the range of {0V, 0.5V}, a write voltage of 1.2 V, both
applied at $V_{IN}$, with measurements taken with a Micromanipulator tungsten
probe tip. The pinched hysteresis loop measured with a 50Hz source is shown in
Fig. 3(a).

First, we tested the circuit in expansion mode. In the critical case of a
write voltage applied to all devices $V_{write}$=1.2V and RE is set HIGH
(1.8V), we show that IR drop is decreased by approximately 22% compared to a
similar planar inference engine in [26]. This is shown by the slower decline
of current output across columns in Fig. 3(b) for CrossStack, where the gold
standard would be a perfectly straight line. Therefore, we verified that
expansion mode reduces line losses for a given number of inputs due to the
shorter length of column wire required. The trade-off is that column wires
must handle twice the current capacity of an equivalent 2-D array. This
demands wide column lines to handle such current capacity without risk of
electromigration. But given that RRAM is integrated in the back end of the
line, minimum thickness of higher layer routing wires may mitigate this.

Designing for deep-net mode opens up susceptibility to leakage currents
through N1 during write mode, concurrently with N2 in read mode (see Fig. 1).
The worst case leakage occurs along the shared column line, when there is a
minimal read current and a maximal write current. The read array memristors
will all be OFF, $R_{r}=100K\Omega$, and all write array memristors will be
ON, $R_{s}=10K\Omega$. To calculate the minimum value for $V_{read}$, we use
0.5 V as the maximum read voltage and assume an input of 7-bit resolution
which requires increments of approximately $500mV/128\approx 4mV$. The output
current of a single cell under these conditions was measured to be 39.6 nA,
which is 1% off the ideal 40 nA. The accumulated leakage current through a
column of 10 memristors being programmed ($V_{IN}$=1.2V; RE=LOW) was
negligible in our experiments, and simulated to be approximately 2.5 pA per
cell (i.e., 2.5pA$\times$10 cells = 25pA column current). This is $6.3\times
10^{-2}$% of the worst-case read-current, and so in the 180-nm process used,
we are able to employ minimum transistor dimensions (W = 450nm, L = 180nm, W/L
= 2.5). Leakage from a single cell is shown as a function of a DC sweep at the
input is measured in Fig. 3(c) and (d), with a Monte Carlo parametric sweep of
resistance overlaid (10k$\Omega\pm$7%, Gaussian distribution across 200
trials).

Figure 4: Transient analysis of a single cell current output during deep-net
mode in a read cycle.

Current read-out measurements taken from a transient analysis under a
switching input during a read-cycle are given in Fig. 4, where current
deviates by 8% from the measured value in the worst-case. This suggests that a
single device should be limited to 3.5-bits. We note this is a limitation not
of the CrossStack architecture, but rather of memristor variability. In the
conservative case of a 1-bit cell across 10 columns, the 10ns read-out is
subsumed within the 25ns programming time as opposed to being added to a
separate cycle after programming. This confirms improvement in speed due to
our in-situ pipelining mechanism by 29%.

##  V Discussion and Conclusion

With two different modes available to the user, how would you decide to use
one over another? In short, deep-net is suited for tasks where speed is
paramount and expansion mode is for tasks involving a large number of inputs,
be it for increased resolution, or to process a very large number of
parameters in a shallow network.

There are three primary motivations in the use of expansion mode: 1) fully-
connected networks typically have a far larger number of connections than in
convolutional layers. Expansion mode doubles the number of possible inputs for
a fixed area which enables increases the number of possible inputs; 2) a
larger number of inputs requires a larger length of column wire. By using
expansion mode to double up on inputs, we reduce wire resistance to half the
original amount for a given number of inputs. More inputs per unit length of
wire enables our crossbar to be resilient to write failures arising from line
resistance IR drops by reducing line losses by 22%. 3) The lack of reliable
analog memory technology makes it hard to perform hardware multiplexing in
analog, and transmitting analog values over long distances or at high speed is
not efficient. Restricting each memristor to one of two conductance values
(i.e., single-bit memristors) means that one would require $log_{2}(n)$
memristors for $n$ bits of precision. For crossbars that use conductance
states of memristors conservatively, digital computing is more desirable but
requires more memristors in a crossbar than analog for the same precision.
Expansion mode facilitates this increase in devices whilst halving crossbar
area. The drawback is that a larger current must be carried through the wire
with more vias.

Deep-net mode is a novel processing scheme where the two crossbar layers are
isolated from one another by appropriately biasing RE, in order to parallelize
read and write operations. It is engaged where speed is of greater importance
than precision. In the most simplistic way (ignoring max-pooling and dropout),
a crossbar performs inference in the following way:

  1. 1.

Write weights to memristor conductances,

  2. 2.

Apply a sub-threshold read voltage,

  3. 3.

Buffer or store read-out signal in memory,

  4. 4.

Write the next hidden layer weights to the crossbar,

  5. 5.

Repeat steps 2-5 until output is generated.

In deep-net mode, CrossStack performs steps 2 and 4 simultaneously which
enables read and write operations to occur together. By the time an output is
generated from step 3, it is ready to be processed by the next hidden layer in
step 4 for a speed increase of 29% over an equivalent 2-D array.

The most prevailing issues in realizing large scale crossbar arrays beyond our
working prototype are mismatch and endurance. At this stage, our current error
rate was 8% which limits the number of bits that can be represented by a
single cell. Subthreshold current was hardly an issue in our design, but Fig.
3(d) shows the nonlinearity of voltage across the memristor under high write
voltages ($V_{IN}>3V$), and leakages will become more prevalent at shorter
channel lengths. Capacitive coupling as a result of high programming voltages
into bit-lines that are reading out poses a degree of risk, but should be
tolerable for higher metal layers given the wider spacing. Allowing for heat
dissipation in a 3D structure is an ongoing challenge and the subject of
significant process-related research, often calling for external heat sinks.
In general, the advantages seem to outweigh the drawbacks and CrossStack
presents a promising methodology for reconfigurable inference acceleration to
adapt to the various types of ANNs being deployed.

## Acknowledgment

This work was supported by the National Research Foundation of Korea grant
funded by the Korean government (MSIT) (No.2020R1F1A1069381).

## References

  * [1] D. Mahajan, R. Girshick, V. Ramanathan, K. He, M. Paluri, Y. Li, A. Bharambe and L. van der Maaten, “Exploring the limits of weakly supervised pretraining,” _Proc. of the European Conf. on Computer Vision (EECV)_ , pp. 181–196, 2018. A. G. Howard, et al., “MobileNets: Efficient convolutional neural networks for mobile vision applications”, arXiv preprint, arXiv:1704.04861.
  * [2] K. H. Kim, S. Jo, S. Gaba and W. Lu, “Nanoscale resistive memory with intrinsic diode characteristics and long endurance,” _Applied Physics Letters_ , vol. 96, no. 5, p. 053106, Feb. 2010.
  * [3] A. Kumar, M. Das, V. Garg, B. S. Sengar, M. T. Htay, S. Kumar, A. Kranti and S. Mukherjee, “Forming-free high-endurance Al/ZnO/Al memristor fabricated by dual ion beam sputtering,” _Applied Physics Letters_ , vol 110, no. 25, p. 253509, Jun. 2017.
  * [4] C.-Y. Lin et al., “Adaptive synaptic memory via lithium ion modulation in RRAM devices”, _Small_ , vol. 16, no. 42, p. 2003964, Oct. 2020.
  * [5] S. Pi, C. Li, H. Jiang, W. Xia, H. Xin, J. J. Yang and Q. Xia, “Memristor crossbar arrays with 6-nm half-pitch and 2-nm critical dimension,” _Nature nanotechnology_ , vol. 14, no. 1, p. 35, Jan. 2019.
  * [6] E. J. Fuller, S. T. Keene, A. Melianas, Z. Wang, S. Agarwal, Y. Li, Y. Tuchman, C. D. James, M. J. Marinella, J. J. Yang, A. salleo and A. A. Talin, “Parallel programming of an ionic floating-gate memory array for scalable neuromorphic computing”, _Science_ , vol. 364, no. 6440, pp. 570–574, May 2019.
  * [7] E. J. Merced-Grafals, N. Dávila, N. Ge, R. S. Williams and J. P. Strachan, “Repeatable, accurate, and high speed multi-level programming of memristor 1T1R arrays for power efficient analog computing applications,” _Nanotechnology_ , vol. 27, no. 36, p. 365202, Aug. 2016.
  * [8] C.-Y. Lin et al., “A high-speed MIM resistive memory cell with an inherent vanadium selector”, _Applied Materials Today_ , vol. 21, p. 100848, Dec. 2020.
  * [9] J. K. Eshraghian, K. Cho, C. Zheng, M. Nam, H. H. C. Iu, W. Lei, and K. Eshraghian, “Neuromorphic Vision Hybrid RRAM-CMOS Architecture,” _IEEE Trans. Very Large Scale Integration (VLSI) Systems_ , vol. 26, no. 12, pp. 2816–-2829, Dec. 2018
  * [10] B. Chakrabarti, M. A. Lastras-Montaño, G. Adam, M. Prezioso, B. Hoskins, M. Payvand, A. Madhavan, A. Ghofrani, L. Theogarajan, K. T. Cheng, D. B. Strukov, “A multiply-add engine with monolithically integrated 3D memristor crossbar/CMOS hybrid circuit,” _Scientific reports_ , vol. 14, no. 7, p. 42429, Feb. 2017.
  * [11] F. Cai, J. M. Correll, S. H. Lee, Y. Lim, V. Bothra, Z. Zhang, M. P. Flynn and W. D. Lu, “A fully integrated reprogrammable memristor-CMOS system for efficient multiply-accumulate operations”, _Nature Electronics_ , vol. 15, no. 1, Jul. 2019.
  * [12] M. R. Azghadi et al., “Complementary Metal-Oxide Semiconductor and Memristive Hardware for Neuromorphic Computing”, _Advanced Intelligent Systems_ , vol. 2, no. 5, Mar. 2020.
  * [13] C. Lammie and M. R. Azghadi, “MemTorch: A simulation framework for deep memristive Cross-Bar architectures”, _2020 IEEE Int. Symp. Circuits and Systems (ISCAS), Sevilla, Spain, Oct. 2020._
  * [14] O. Krestinskaya, K. N. Salama and A. P. James, “Analog Backpropagation Learning Circuits for Memristive Crossbar Neural Networks”, _2018 IEEE Int. Symp. Circuits and Systems (ISCAS), Florence, Italy, May 2018_
  * [15] A. Serb, J. Bill, A. Khiat, R. Berdan, R. Legenstein and T. Prodromakis, “Unsupervised learning in probabilistic neural networks with multi-state metal-oxide memristive synapses,” _Nature communications_ , vol. 7, p. 12611, Sep. 2016.
  * [16] C. J. Chevallier, C. H. Siau, S. F. Lim, S. R. Namala, M. Matsuoka, B. L. Bateman and D. A. Rinerson, “A 0.13$\mu$m 64Mb multi-layered conductive metal-oxide memory”, _2010 IEEE Int. Solid-State Circuits Conf. (ISSCC)_ , pp. 260–261, Feb. 2010.
  * [17] C. H. Wang, Y. H. Tsai, K. C. Lin, M. F. Chang, Y. C. King, C. J. Lin, S. S. Sheu, Y. S. Chen, Y. H. Lee, F. T. Chen and M. J. Tsai, “Three-dimensional $4F^{2}$ ReRAM cell with vertical BJT driver by CMOS logic compatible process,” _IEEE Trans. on Electron Devices_ , vol. 58, no. 8, pp. 2466–2472, Jul 2011.
  * [18] S. H. Jo, T. Kumar, S. Narayanan, W. D. Lu and H. Nazarian, “3D-stackable crossbar resistive memory based on field assisted superlinear threshold (FAST) selector,” _2014 IEEE Int. Electron Devices Meeting_ , Washington, DC, USA, pp. 6–7, Dec. 2014.
  * [19] J. Hong, M. Stone, B. Navarrete, K. Luongo, Q. Zheng, Z. Yuan, K. Xia, N. Xu, J. Bokor, L. You and S. Khizroev, “3D multilevel spin transfer torque devices,” _Applied Physics Letters_ , vol. 112, no. 11, p. 112402, Mar. 2018.
  * [20] M. R. Azghadi et al., “Hardware implementation of deep network accelerators towards healthcare and biomedical applications”, _IEEE Trans. Biomedical Circuits and Systems_ , vol. 14, no. 6, pp. 1138-1159, Dec. 2020.
  * [21] S. Baek, J. K. Eshraghian, S.-H. Ahn, A. James and K. Cho, “A memristor-CMOS Braun multiplier array for arithmetic pipelining”, _2019 26th IEEE Int. Conf. Electronics, Circuits and Systems (ICECS), Genoa, Italy, pp. 735–738, Nov. 2019._
  * [22] R. Fastow, K. Hasnat, P. Majhi and O. Jungroth, “Three-dimensional (3D) memory with shared control circuitry using wafer-to-wafer bonding,” _United States Patent Application_ , US 16/011, 139, Feb. 2019.
  * [23] J. K. Eshraghian, K. R. Cho, H. H. C. Iu, T. Fernando, N. Iannella, S. M. Kang and K. Eshraghian, “Maximization of crossbar array memory using fundamental memristor theory,” _IEEE Trans. on Circuits and Syst. II: Express Briefs_ , vol. 64, no. 12, pp. 1402–1406, Oct. 2017.
  * [24] C. Li, M. Hu, Y. Li, H. Jiang, N. Ge, E. Montgomery, J. Zhang, W. Song, N. Davila, C. E. Graves, Z. Li, J. P. Strachan, P. Lin, Z. Wang, M. Barnell, Q. Wu, R. S. Williams, J. J. Yang, and Q. Xia, “Analogue signal and image processing with large memristor crossbars”, _Nature Electronics_ , vol. 1, no. 1, pp. 52-–59, Jan. 2018.
  * [25] J. K. Eshraghian, S. M. Kang, S. Baek, G. Orchard, H. H. C. Iu, and W. Lei, “Analog weights in ReRAM DNN Accelerators”, _2019 IEEE Artificial Circuits and Systems Conference_ , Mar. 2019.
  * [26] M. Hu, J. P. Strachan, Z. Li, E. M. Grafals, N. Davila, C. Graves, S. Lam, N. Ge, J. J. Yang, R. S. Williams, “Dot-product engine for neuromorphic computing: Programming 1T1M crossbar to accelerate matrix-vector multiplication,” _Proc. of the 53rd Annual Design Automation Conference_ , p. 19, Jun. 2016.
  * [27] A. Shafiee, A. Nag, N. Muralimanohar, R. Balasubramonian, J. P. Strachan, M. Hu, R. S. Williams and V. Srikumar, “ISAAC: A convolutional neural network accelerator with in-situ analog arithmetic in crossbars,” _ACM SIGARCH Computer Architecture News_ , vol. 44, no. 3, pp. 14–26, Oct. 2016.
  * [28] S. Stathopoulos, A. Khiat, M. Trapatseli, S. Cortese, A. Serb, I. Valov and T. Prodromakis, “Multibit memory operation of metal-oxide bi-layer memristors,” _Scientific Reports_ , vol. 7, no. 17532, Dec. 2017.
  * [29] L. Song, X. Qian, H. Li and Y. Chen, “Pipelayer: A pipelined ReRAM-based accelerator for deep learning,” _2017 IEEE Int. Symp. on High Performance Computer Architecture (HPCA)_ , pp. 541–552, Feb. 2017.
  * [30] H. Valavi, P. J. Ramadge, E. Nestler and N. Verma, “A 64-Tile 2.4-Mb In-Memory-Computing CNN accelerator employing charge-domain computer,” _IEEE Journal of Solid-State Circuits_ , vol. 54, no. 6, pp. 1789-1799, Mar. 2019.
  * [31] G. C. Adam, B. D. Hoskins, M. Prezioso, F. Merrikh-Bayat, B. Chakrabarti and D. B. Strukov, “3-D memristor crossbars for analog and neuromorphic computing applications,” _IEEE Transactions on Electron Devices_ , vol. 64, no. 1, pp. 312–318, Jan. 2017.