\section{Introduction}
For any two matrices $A$ and $B$, Lie proved the celebrated product formula
\begin{equation}\label{Lie}
\ep{A+B} = \lim_{n\to\infty}\left(\ep{A/n}\ep{B/n}\right)^n.
\end{equation}
There is a long line of similar formulae of increasing generality, pioneered by Trotter~\cite{trotter1959product}, simplified by Chernoff~\cite{chernoff1968note} for semigroups on Banach spaces. In the particular setting of quantum mechanics where $A,B$ are densely defined unbounded self-adjoint operators and $\exp(i t A),\exp(i t B)$ and $\exp(i t (A+B))$ are the corresponding unitary groups, the product formula was proved under general assumptions by Kato~\cite{kato1978trotter} and Ichinose~\cite{ichinose1980product}. It plays a crucial role in functional integration, see in particular~\cite{nelson1964feynman}. For related results in the context of the quantum Zeno effect, we refer to~\cite{exner2005product}, and point further to~\cite{neidhardt1999trotter} for the case of Gibbs semigroups. 

The recent interest in proving general product formulae with explicit control of the rate of convergence has been motivated by two related developments in many-body quantum systems. On the one hand in quantum information theory, operator products arise as quantum circuits and a product formula is interpreted as a \emph{simulation algorithm} for the time evolution of a quantum system~\cite{Llo96,BerAhoCleSan06,wiebe2011simulating}. On the other hand in condensed matter physics, operator products are referred to as \emph{finite depth quantum circuits} and play a central role in the classification of gapped phases~\cite{chen2010local}, as they can be used to define the very notion of equivalence of states. In both cases, the concepts have recently been tested experimentally, see e.g.~\cite{NatureSimulationExp,SmiKimPolKno19}.

In both applications, the rate of convergence of the product formula to the full dynamics is of crucial importance: for quantum simulation because it determines the number of quantum gates required to simulate to a given error, for gapped phases because it relates to the degree of entanglement of ground states. Beyond the mere scaling, sharp constants are essential and may prove fatal in a many-body setting. Indeed, for a lattice system having $N$ degrees of freedom, the error diverges as $N\to\infty$, even in the strong operator topology, which is the natural topology as soon as $A,B$ are unbounded. This is related to the infrared catastrophe: If two states are locally close to each other but the error extends to spatial infinity, then they are in fact orthogonal.

In this work, we consider $d$-dimensional quantum lattice systems in the \emph{infinite volume limit}. The dynamics is an automorphism group $t\mapsto\tau^\Phi_t$ generated by a local Hamiltonian formally given by
\begin{equation}\label{Hamiltonians}
H = \sum_{X}\Phi(X) = \sum_{j=1}^k K_j
\end{equation}
where the $K_j$'s correspond to an arbitrary grouping of the interaction terms. We provide product formulae and prove explicit bounds in the strong topology of the operator algebra and for sufficiently localized observables: For any $m\in{\mathbb N}$ there is a product automorphism denoted $\pi_{t,n}^{(m)}$ such that
\begin{equation*}%
\left\Vert  \tau^\Phi_t(O) - \pi_{t,n}^{(m)}(O) \right\Vert
\leq C_{m,t}(O) n^{-m}
\end{equation*}
for any almost local observable. The constant $C_{m,t}(O)$ depends on both the observable~$O$ and the Hamiltonian. As in~(\ref{Lie}), $\pi_{t,n}^{(m)}$ are compositions of the individual dynamics generated by each $K_j$ individually. The general form of the product $\pi_{t,n}^{(m)}$ was proposed by Suzuki~\cite{suzuki1991general} although in the Hilbert space setting, see also~\cite{hatano2005finding}, and used recently by Childs et al.~\cite{PhysRevX.11.011020}, but only for finite systems with bounds that diverge in the volume.

From a technical point of view, we find it convenient to consider almost exclusively the \emph{interaction norms}, see Definiton~(\ref{def:int norm}) below, that measure the local size of extensive observables of the type~(\ref{Hamiltonians}). We extend their definition to be able to consider interactions localized in (possibly infinite) subsets of the lattice (see also~\cite{monaco2019adiabatic}) and remark that this construction is well-suited to discuss \emph{almost local observables}. Crucially, the $*$-derivation formally given by
\begin{equation*}
\sum_{X}i[\Phi(X),\cdot]
\end{equation*}
extends to a bounded operator on the set of almost local observables equipped with the interaction norm. 

We further wish to point out two related results. Firstly, a slightly different approach to a product formula was taken in~\cite{sahinoglu2021hamiltonian}, focussing on the `quasi-adiabatic' properties of product formulae, namely the error when projected onto the ground state space. Secondly, a similar spatial product factorization with sharp error bounds was derived in~\cite{haah2021quantum}: it is not based on the Trotter strategy but it uses rather directly the Lieb-Robinson bound, see also~\cite{bachmann2017local}.

While the results hold for a general decomposition~(\ref{Hamiltonians}), in applications the factors $K_j$ will be chosen so as to be commuting Hamiltonians, namely each $K_j$ is a sum of \emph{mutually commuting} interaction terms $\Phi_j(X)$, i.e. $[\Phi_j(X),\Phi_j(X')] = 0$ for all $X,X'$. Such Hamiltonians have the property that the corresponding automorphism $\tau^{\Phi_j}_t(O)$ is strictly local in that the support of the observable $O$ grows at most by the range of the interaction, uniformly in the time $t$. Propagation, which is obviously present in the full dynamics $\tau_t^{\Phi}(O)$ arises then through the alternating action of the automorphisms $\{\tau^{\Phi_j}:j=1,\ldots,k\}$. While the Lieb-Robinson bound is at the heart of the proofs, the product formulae provide a very clear picture of the mechanism of propagation. 

In the context of quantum simulation, much attention is given to the error made in the approximation upon truncation of the product formula to a finite number of terms. As we shall see, the error has a complicated dependence on a number of parameters and we shall discuss this in detail later. We already point out now (i) that the error is exponential
in the total time $t$, 
(ii) that the number of factors in~$\pi_{t,n}^{(m)}$ is proportional to $n$ and to the number $k$ of factors in the decomposition of the Hamiltonian, and that it is exponential in the order $m$ of approximation, 
and (iii) that unlike in the original Trotter product formula, the times involved in the various factors of~$\pi_{t,n}^{(m)}$ are not all equal, although they are all of order $\frac{t}{n}$; in fact, the time evolution runs backwards for a fraction of the factors, giving rise to a fractal path, see Figure~\ref{fig:Nice plots} at the end of Section~\ref{sec:ArbOrder}. 

\section{Setting} \label{sec:Setting}

Let $(\Gamma,\textbf{d})$ be a metric graph, where $\textbf{d}$ is the graph distance. We assume that $\Gamma$ is $d$-dimensional in the sense that $\sup_{x\in\Gamma}\vert \{y\in\Gamma: \textbf{d}(x,y) = r\}\vert =\kappa (1+r)^{d-1}$. For any subset $X\subset\Gamma$, we define for any $r>0$ the set $X^{(r)} = \{x\in\Gamma:\textbf{d}(x,X)\leq r\}$ which is an $r$-fattening of the set $X$.

To each site $x \in \Gamma$ we associate a complex Hilbert space ${\mathcal H}_x$ and define for any finite $\Lambda \subset \Gamma$, 
\begin{align*}
	{\mathcal H}_\Lambda := \bigotimes_{x \in \Lambda} {\mathcal H}_x \, \quad \textrm{ and } \quad {\mathcal A}_\Lambda := \bigotimes_{x \in \Lambda} {\mathcal B}({\mathcal H}_x) \,,
\end{align*}
where ${\mathcal B}({\mathcal H}_x)$ denotes the bounded linear operators over ${\mathcal H}_x$. Moreover, we identify $A \in {\mathcal A}_{\Lambda_0}$ with $A \otimes \bbI_{\Lambda\setminus\Lambda_0} \in {\mathcal A}_\Lambda$ whenever $\Lambda_0\subset\Lambda$. 
With this we can inductively define the algebra of local observables 
\begin{align*}
	{\mathcal A}_{\rm loc} := \bigcup_{\Lambda \in {\mathcal P}_{\rm fin}(\Gamma)} {\mathcal A}_{\Lambda}\,,
\end{align*}
where the union is taken over ${\mathcal P}_{\rm fin}(\Gamma)$, the set of all finite subsets of $\Gamma$. If $O\in{\mathcal A}_{\rm loc}$, then $\mathrm{supp}(O)$ is the smallest set $X$ such that $O\in{\mathcal A}_X$. The completion of ${\mathcal A}_{\rm loc}$ with respect to the norm topology is a $C^*$-algebra which is called the \textit{quasi-local algebra} and we denote it by ${\mathcal A}$. 

\subsection{Interactions and Hamiltonians}
\begin{defn}
An \textit{interaction} is a map $\Phi : {\mathcal P}_{\rm fin}(\Gamma) \to {\mathcal A}_{\rm loc}$ such that 
$$
\Phi(X)\in{\mathcal A}_X,\qquad	\Phi(X) = \Phi(X)^*,
$$ 
for all $X \in {\mathcal P}_{\rm fin}(\Gamma)$. 
\end{defn}

We turn the set of interactions into a Banach space in the following way. Let $0 < p < 1$ and let
$$
	\xi_{b} : [0,\infty) \to (0,\infty)\,, \qquad \xi_{b}(x) = e^{-bx^p}\,,
$$ 
for any $b > 0$. The function $\xi_{b}$ is a decreasing, logarithmically superadditive function, namely $\xi_b(x+y)\geq \xi_b(x)\xi_b(y)$, that is summable in the following sense	
\begin{equation}\label{xi is integrable}
\|\xi_{b}\|_1 := \sup_{y \in \Gamma} \sum_{x \in \Gamma} \xi_{b}(\textbf{d}(x,y)) < \infty,
\end{equation}
since $\Gamma$ is finite dimensional.
\begin{defn}\label{def:int norm}
 Let $D(X):= \max\{\textbf{d}(x,y)\,:\, x, y \in X\}$ denote the diameter of the set $X\subset \Gamma$.  The \textit{interaction norm} of an interaction $\Phi$ is given by 
\begin{align}\label{eq:StandardNormInteraction}
	\vvvert\Phi\vvvert_b := \sup_{x \in \Gamma} \sum_{\substack{X \in {\mathcal P}_{\rm fin}(\Gamma): \\ x \in X}} \frac{\|\Phi(X)\|}{\xi_{b}(D(X))}.
\end{align}
For fixed $b>0$ we denote the Banach space of interactions with finite $\vvvert\cdot\vvvert_b$-norm
by ${\mathcal B}_b$ and set 
\begin{align*}
	{\mathcal B} := \bigcup_{b >0} {\mathcal B}_b\,.

\end{align*}
\end{defn}
The norm indicates both the rate of decay of the interaction as well as its intensity in the sense that the total interaction at any given $x \in \Gamma$ is bounded by the interaction norm:
\begin{equation*}
	\sup_{x\in\Gamma}\Big\|\sum_{\substack{X \in {\mathcal P}_{\rm fin}(\Gamma): \\ x \in X}} \Phi(X) \Big\|
	\leq \vvvert \Phi\vvvert_b\,.	
\end{equation*}
We will also need the notion of an interaction that is \textit{almost localized} in some set $Z \subset \Gamma$.
\begin{defn}
Let $Z \subset \Gamma$ and let $D_{Z}(X) := D(X) + \textbf{d}(X,Z)$. Let
\begin{align}\label{eq:RefinedNormInteraction}
	\vvvert\Phi\vvvert_{b,Z} := \sup_{x \in \Gamma} \sum_{\substack{X \in {\mathcal P}_{\rm fin}(\Gamma): \\ x \in X}} \frac{\|\Phi(X)\|}{\xi_b(D_Z(X))}.
\end{align}
We denote ${\mathcal B}_b(Z)$ the corresponding Banach space of interactions and define 
\begin{align*}
	{\mathcal B}(Z) := \bigcup_{b > 0} {\mathcal B}_{b}(Z) \,.
\end{align*} 
\end{defn}
Note that the replacement of the diameter by $D_{Z}(X)$ captures the decay of the interaction $\Phi(X)$ in the size of $X$ and in the distance from $X$ to $Z$. If $\Phi \in {\mathcal B}_b(Z)$ then the total contribution of $\Phi$ at a given point $x$ is not only finite, but it decays with the distance of $x$ to $Z$,
\begin{align}\label{eq:TotalIntAtxFiniteAndDecays}
	\Big\|\sum_{\substack{X \in {\mathcal P}_{\rm fin}(\Gamma): \\ x \in X}} \Phi(X) \Big\| \leq \vvvert \Phi\vvvert_{b,Z} \,\xi_b(\textbf{d}(x,Z))\,.
\end{align}
Indeed, we first write
\begin{equation}\label{usual}
\Big\Vert \sum_{\substack{X\in{\mathcal P}_{\mathrm{fin}}(\Gamma):\\x\in X}}\Phi(X)\Big\Vert
\leq \sum_{\substack{X\in{\mathcal P}_{\mathrm{fin}}(\Gamma):\\x\in X}} \frac{\Vert \Phi(X)\Vert }{\xi_b(D_Z(X))} \xi_b(D_Z(X)).
\end{equation}
Given $x\in X$, let $z,z_0\in Z$ and $x_0\in X$ be such that $\textbf{d}(x,Z) = \textbf{d}(x,z)$, $\textbf{d}(X,Z) = \textbf{d}(x_0,z_0)$. Then,
\begin{equation}\label{dxZ vx dXZ}
\textbf{d}(x,Z) = \textbf{d}(x,z) \leq \textbf{d}(x,z_0) \leq \textbf{d}(x,x_0) + \textbf{d}(x_0,z_0) \leq D_Z(X).
\end{equation}
Since $\xi_b$ is a decreasing function, we have that $\xi_b(D_Z(X)) \leq \xi_b(\textbf{d}(x,Z))$ which yields~(\ref{eq:TotalIntAtxFiniteAndDecays}) when plugged in~(\ref{usual}).

In general, the sum of an interaction is not convergent in ${\mathcal A}$, but such a \emph{local Hamiltonian} defines a dynamics and a *-derivation on ${\mathcal A}$. Moreover, if an interaction is almost localized in a finite set, then the sum is convergent and defines an \emph{almost local observable} in ${\mathcal A}$. 

We start with the latter.
\begin{lemma}\label{lemma: finite Z interaction}
Let $Z\in{\mathcal P}_{\mathrm{fin}}(\Gamma)$ and $\Phi\in{\mathcal B}(Z)$. 
The sum
\begin{equation}\label{G of Phi}
G^\Phi := \sum_{X\in{\mathcal P}_{\mathrm{fin}}(\Gamma)}\Phi(X)
\end{equation}
is convergent in ${\mathcal A}$ and 
\begin{equation*}
\Vert G^\Phi\Vert \leq C \vvvert \Phi \vvvert_{b,Z}\vert Z\vert 
\end{equation*}
for all $b>0$.
\end{lemma}
\begin{proof}
Let $b>0$ be such that $\Phi\in{\mathcal B}_b(Z)$. Then for any $\Lambda\in {\mathcal P}_{\mathrm{fin}}(\Gamma)$,
\begin{equation*}
\sum_{\substack{X\in{\mathcal P}_{\mathrm{fin}}(\Gamma):\\X\cap\Lambda^c \neq \varnothing}}\Vert \Phi(X) \Vert \leq \sum_{x\in\Lambda^c}\sum_{X\ni x }\Vert \Phi(X) \Vert 
\end{equation*}
and we conclude by~(\ref{eq:TotalIntAtxFiniteAndDecays}) that
\begin{equation*}
\Big \Vert \sum_{\substack{X\in{\mathcal P}_{\mathrm{fin}}(\Gamma):\\X\cap\Lambda^c \neq \varnothing}} \Phi(X) \Big \Vert\leq \vvvert \Phi \vvvert_{b,Z} \sum_{x\in\Lambda^c} \xi_b(\textbf{d}(x,Z)).
\end{equation*}
The integrability~(\ref{xi is integrable}) of $\xi_b$ and the finiteness of $Z$ imply that $\lim_{\Lambda\nearrow\Gamma}\sum_{x\in\Lambda^c} \xi_b(\textbf{d}(x,Z)) = 0$. The second claim follows from
\begin{equation*}
\Vert G^\Phi\Vert
\leq \sum_{x\in\Gamma}\sum_{\substack{X\in{\mathcal P}_{\mathrm{fin}}(\Gamma):\\x\in X}}\frac{\Vert \Phi(X)\Vert}{\xi_b(D_Z(X))}\xi_b(D_Z(X))
\leq \vvvert \Phi \vvvert_{b,Z} \sum_{x\in\Gamma} \xi_b(\textbf{d}(x,Z))
\end{equation*}
where we used~(\ref{dxZ vx dXZ}). We decompose $\sum_{x\in\Gamma} = \sum_{n=0}^\infty\sum_{x:\textbf{d}(x,Z) = n}$ to finally get the bound
\begin{equation*}
\Vert G^\Phi\Vert
\leq\vvvert \Phi \vvvert_{b,Z}\sum_{n=0}^\infty \sum_{z\in Z}\sum_{x:\textbf{d}(x,z) = n}  \xi_b(n)
\leq C \vvvert \Phi \vvvert_{b,Z} \vert Z\vert
\end{equation*}
where $C = \kappa\sum_{n=0}^\infty(1+n)^{d-1}\xi_b(n)$ is convergent.
\end{proof}

\begin{defn}
An \emph{almost local observable} is an elements $O\in{\mathcal A}$ for which there exists $Z\in{\mathcal P}_{\mathrm{fin}}(\Gamma)$ and a $\Phi\in{\mathcal B}(Z)$ such that $O = G^\Phi$. 
We denote the set of almost local observables by ${\mathcal L}(Z)$, respectively ${\mathcal L}_b(Z)$ whenever the rate $b$ is fixed.
\end{defn}

Slightly abusing language, we shall call $Z$ the (almost) support of $O\in {\mathcal L}(Z)$. We further point out that ${\mathcal A}_Z\subset\bigcap_{b>0}{\mathcal L}_b(Z)$. Lemma~\ref{lemma: finite Z interaction} shows that the ratio of the operator norm of an almost local observable to the interaction norm of the corresponding interactions is proportional to the volume of the support.  

Let us now turn to interactions that are supported in the whole set $\Gamma$. They are locally finite, but the fact that they are extensive implies that a sum as in~(\ref{G of Phi}) is divergent. This suggests the following definition.
\begin{defn}\label{def:Hamiltonian}
A family of self-adjoint operators $H = \{H_\Lambda \, : \, \Lambda \in {\mathcal P}_{\rm fin}(\Gamma)\}$ with $\supp(H_\Lambda) = \Lambda$ is a \textit{local Hamiltonian} if there exists an interaction $\Phi \in {\mathcal B}$ such that 
\begin{align*
	H_\Lambda = \sum_{X \subset \Lambda} \Phi(X)\,.
\end{align*}
We denote the set of local Hamiltonians by ${\mathcal L}$.
\end{defn}
Let $H$ be a local Hamiltonian associated with an interaction $\Phi\in{\mathcal B}_b$. Let $(\Lambda_n)_{n\in{\mathbb N}}$ be an increasing and absorbing sequence of finite sets. For an observable $O\in{\mathcal A}_Z$, we have that if $n>m$
\begin{align*}
 \Vert [H_{\Lambda_n} - H_{\Lambda_m}, O] \Vert
&\leq  \sum_{\substack{X\subset\Lambda_n: \\ X\cap(\Lambda_m^c\cap Z)\neq\varnothing}}2\Vert O\Vert \Vert \Phi(X)\Vert 
\leq 2\Vert O\Vert \sum_{x\in Z}\sum_{\substack{X\ni x: \\ X\cap\Lambda_m^c\neq\varnothing}}
\frac{\Vert \Phi(X)\Vert}{\xi_b(D_Z(X))}\xi_b(D_Z(X)) \\
&\leq 2\Vert O\Vert \vert Z\vert\vvvert \Phi\vvvert_b \xi_b(\textbf{d}(Z,\Lambda_m^c))
\end{align*}
since $D_Z(X) = D(X) \geq \textbf{d}(Z,\Lambda_m^c)$. Hence $(i [H_{\Lambda_n} , O ])_{n\in{\mathbb N}}$ is a Cauchy sequence in ${\mathcal A}$. 

\subsection{Strongly continuous dynamics and derivations}

Let $H\in{\mathcal L}$ with interaction $\Phi\in{\mathcal B}$. Since the finite volume dynamics $\ep{itH_\Lambda}O \ep{-itH_\Lambda}$ satisfyies a Lieb-Robinson bound, it follows \cite{lieb1972finite} that the limit  
\begin{align*}
	\tau_t^{\Phi}(O) = \lim_{\Lambda \nearrow \Gamma} \ep{it H_\Lambda} O \ep{-it H_\Lambda}	
\end{align*}
exists on ${\mathcal A}_{\rm loc}$ and that it extents to a strongly continuous group of $^*$-automorphisms of ${\mathcal A}$. 
The corresponding generator $\delta^\Phi$ of the dynamics $\tau_t^{\Phi}$ is given by 
\begin{align*}
	\frac{d}{dt} \tau_t^\Phi(O) = \tau_t^\Phi(\delta^\Phi(O))\,.
\end{align*}
A core for $\delta^\Phi$ is the local algebra ${\mathcal A}_{\rm loc}$ where $\delta^\Phi$ is explicitly given as the limit 
\begin{equation}\label{delta on Aloc}
\delta^\Phi(O) = \lim_{\Lambda\to\Gamma}i[H_\Lambda,O] = \sum_{X\in{\mathcal P}_{\mathrm{fin}}(\Gamma)}i[\Phi(X),O],
\end{equation}
by the remark immediately after Definition~\ref{def:Hamiltonian}.

For a general Hamiltonian in ${\mathcal L}$, neither $\tau_t^\Phi(O)$ nor $\delta^\Phi(O)$ is strictly local even if $O\in{\mathcal A}_{\mathrm{loc}}$. However, we shall prove later that for any finite $Z$, ${\mathcal L}(Z)$ belongs to the domain of $\delta^\Phi$ and it is invariant under the action of the derivation. Similarly, the Lieb-Robinson bound implies that ${\mathcal L}(Z)$ is invariant under the action of~$\tau_t^\Phi$.

Let $r\in{\mathbb N}_0$. For any local observable $O \in {\mathcal A}_Z$ we write
\begin{align}\label{eq:decompDynamics}
	\tau_t^\Phi(O) = {\mathbb E}_{Z^{(r)}}(\tau_t^\Phi(O)) + \sum_{n=r}^\infty \left({\mathbb E}_{Z^{(n+1)}}(\tau_t^\Phi(O)) - {\mathbb E}_{Z^{(n)}}(\tau_t^\Phi(O))\right)
\end{align}
where  ${\mathbb E}_{X}$ is the projection onto the subalgebra ${\mathcal A}_{X}$. Now, if $\Phi \in {\mathcal B}_{b'}$, then the Lieb-Robinson bound for $\tau_t^\Phi$ implies that, for any $0<b''<b'$ and $n \in {\mathbb N}$,
\begin{align}\label{eq:AppOfLiebRobBound}
	\|\left({\mathbb E}_{Z^{(n)}} - \textrm{id}\right)(\tau_t^\Phi(O))\|\leq  C \|O\| |Z| \ep{v|t|}\xi_{b''}(n)\,,
\end{align} 
where $C$ depends only on $b'$ and the choice of $b''$ and the Lieb-Robinson velocity $v$ is proportional to $\vvvert \Phi\vvvert_{b'}$.

\begin{prop}\label{prop:DynamicsIsAlmostLocalObservable}
Let $Z\subset\Gamma$ and $O\in{\mathcal L}_b(Z)$. Let $\Phi \in {\mathcal B}_{b'}$ and let $\tau_t^\Phi$ be the corresponding dynamics. Then $\tau^\Phi_t(O)\in{\mathcal L}_{b''}\left(Z\right)$ for any $b''<\min(b,2^{-p} b')$. If $\Psi$, respectively $\tau_t^\Phi(\Psi)$, is the interaction associated with $O$, respectively $\tau^\Phi_t(O)$, then 
\begin{equation*}
\vvvert \tau^\Phi_t(\Psi)\vvvert_{b'',Z} \leq C_\tau \ep{ v\vert t\vert}  \vvvert \Psi\vvvert_{b,Z}
\end{equation*}
where $v$ is the Lieb-Robinson velocity of $\tau^\Phi_t$ and $C_\tau$ depends only on $b,b',b''$ and the spatial dimension $d$.
\end{prop}
\begin{proof}
Let $\Psi\in{\mathcal B}_{b}(Z)$ be the interaction associated with $O$, namely $O = \sum_{X\in{\mathcal P}_{\mathrm{fin}}(\Gamma)}\Psi(X)$. We construct an interaction, denoted $\tau_t^\Phi(\Psi)$, such that $\tau_t^\Phi(O) = \sum_{X\in{\mathcal P}_{\mathrm{fin}}(\Gamma)}\tau_t^\Phi(\Psi)(X)$ as follows. We decompose each $\tau_t^\Phi(\Psi(Y))$ according to~(\ref{eq:decompDynamics}) with $r=0$ and gather contributions supported in a given set $X$ to get
\begin{equation*}
\tau_t^\Phi(\Psi)(X) := {\mathbb E}_{X}(\tau_t^\Phi(\Psi(X))) +
\sum_{n=1}^\infty\sum_{\substack{Y\in{\mathcal P}_{\mathrm{fin}}(\Gamma): \\ X=Y^{(n)}}}\left({\mathbb E}_{Y^{(n)}}-{\mathbb E}_{Y^{(n-1)}}\right)(\tau_t^\Phi(\Psi(Y))).
\end{equation*}
Using~(\ref{eq:AppOfLiebRobBound}), all terms in the sum are bounded above by
\begin{equation}\label{E-E}
\left\Vert ({\mathbb E}_{Y^{(n)}}-{\mathbb E}_{Y^{(n-1)}})(\tau_t^\Phi(\Psi(Y)))\right\Vert
\leq C\Vert \Psi(Y)\Vert \left\vert Y^{(n)}\right\vert\ep{v|t|}\xi_{\tilde b}(n),
\end{equation}
for any $\tilde b<b'$. The first term is bounded above by $\Vert \Psi(X)\Vert$ since ${\mathbb E}_{X}$ is a projection; for simplicity, we shall rather use the estimate~(\ref{E-E}) with $n=0$ in the following. 

Let $x\in\Gamma$. We claim that $\sum_{\substack{X \in {\mathcal P}_{\rm fin}(\Gamma) : \\ x\in X}}\frac{\Vert \tau_t^\Phi(\Psi)(X)\Vert }{\xi_{b''}(D_{Z}(X))}$ is uniformly bounded. For any $r\in{\mathbb N}$, we can first replace $\sum_{X\ni x}\sum_{Y: Y^{(r)} = X}$ by $\sum_{Y: Y^{(r)}\ni x}$. We then note that such sets $Y$ are so that $B_r(x)\cap Y\neq\varnothing$ so that
\begin{equation*}
\sum_{\substack{X \in {\mathcal P}_{\rm fin}(\Gamma) : \\ x\in X}}\frac{\Vert \tau_t^\Phi(\Psi)(X)\Vert }{\xi_{b''}(D_{Z}(X))}
\leq C\sum_{n=0}^\infty \sum_{y\in B_{n}(x)}\sum_{\substack{Y \in {\mathcal P}_{\rm fin}(\Gamma) : \\ y\in Y}}
\frac{ \Vert \Psi(Y)\Vert}{\xi_{b}(D_Z(Y))}\vert Y\vert \xi_{b-b''}(D_Z(Y))\frac{\xi_{\tilde b}(n)}{\xi_{b''}(2n)}(n+1)^d\ep{v\vert t\vert}.
\end{equation*}
where we used that $D_Z(Y^{(n)})\leq D_Z(Y) + 2n$ and the logarithmic superadditivity of $\xi_{b''}$. Hence,
\begin{equation*}
\sup_{x\in\Gamma}\sum_{\substack{X \in {\mathcal P}_{\rm fin}(\Gamma) : \\ x\in X}}\frac{\Vert \tau_t^\Phi(\Psi)(X)\Vert }{\xi_{b''}(D_{Z}(X))}\leq C \vvvert\Psi\vvvert_{b,Z}\ep{v\vert t\vert}
\end{equation*}
where the constant is proportional to $\sup_{Y\in{\mathcal P}_{\mathrm{fin}}(\Gamma)}\{\vert Y\vert \xi_{b-b''}(D_Z(Y))\}$ and $\sum_{n=0}^\infty\frac{\xi_{\tilde b}(n)}{\xi_{b''}(2n)}(1+n)^{2d}$ which are both finite, because $b''<b$ for the first one and we can pick $b''<2^{-p}\tilde b$ since $b''<2^{-p}b'$ to ensure the convergence of the series.
\end{proof}
\begin{rem}
(i) A less detailed but clearer way to formulate the result would be that the *-subalgebra of almost local observables supported in $Y$ is an invariant subspace for $\tau_t^\Phi$ for $t$ in a compact interval.\\
(ii) One could wish to take the propagation into account in this result and prove rather that $\tau_t^\Phi({\mathcal L}(Z)) \subset {\mathcal L}(Z^{(v\vert t\vert)})$, at least in the case of an interaction $\Phi$ that decays exponentially. This is of course true as ${\mathcal L}(Z^{(v\vert t\vert)})$ and ${\mathcal L}(Z)$ are equal as sets, but equipped with different norms. Since however the bound would still be superpolynomially large in time (because the support of each individual interaction term grows with time and hence the decay rate of the interaction does worsen), there is no real gain in doing so.
\end{rem}

The derivation $\delta^\Phi$ associated with a local interaction is in general unbounded on ${\mathcal A}$ and accordingly not everywhere defined. As pointed out earlier, ${\mathcal A}_{\mathrm{loc}}$ is a core on which it is given explicitly as the limit of a commutator. We prove that $\delta^\Phi$ extends to the set of almost local observables and that, as for the automorphism $\tau^\Phi_t$, the sets ${\mathcal L}(Y)$ are invariant under the action of $\delta^\Phi$. 

Instead of considering $\delta^\Phi$ as an unbounded operator on observables, we find it more convenient to define it on the Banach space of interactions ${\mathcal B}$ on which it is a bounded linear operator.
\begin{defn}\label{def:delta of psi}
Let $\Phi,\Psi\in{\mathcal B}$. The interaction $\delta^\Phi(\Psi)$ is defined by
\begin{equation}\label{eq:DerivationForTwoInteractions}
	\delta^\Phi(\Psi)(X) := \sum_{\substack{Y,Y' \in {\mathcal P}_{\rm fin}(\Gamma): \\ Y \cap Y' \ne \varnothing, \, Y\cup Y' = X}} i [\Phi(Y),\Psi(Y')]
\end{equation}
for any $X\in {\mathcal P}_{\rm fin}(\Gamma)$. 
\end{defn}
\noindent Note that the condition $Y \cap Y' \ne \varnothing$ is only for clarity since the commutator vanishes if it is not satisfied.
\begin{rem}
If $O\in{\mathcal A}_Z$ and $\Psi$ is the interaction trivially associated with it, namely $\Psi(Z) = O$  and $\Psi(X) = 0$ otherwise, then the definition above yields an interaction such that 
\begin{equation*}
\sum_{X\in {\mathcal P}_{\rm fin}(\Gamma)}\delta^\Phi(\Psi)(X)
=\sum_{Y\in {\mathcal P}_{\rm fin}(\Gamma)} i [\Phi(Y),O] = \delta^\Phi(O)
\end{equation*}
as in~(\ref{delta on Aloc}), justifying the notation $\delta^\Phi$. 
\end{rem}
\begin{prop}\label{prop:ActionOfDeltas}
Let $Z\subset \Gamma$ and let $\Psi\in{\mathcal B}_b(Z)$. Let $\Phi \in {\mathcal B}_{b'}$. If $\delta^\Phi(\Psi)$ is defined as in~(\ref{eq:DerivationForTwoInteractions}), then $\delta^\Phi(\Psi)\in{\mathcal B}_{b''}(Z)$ for any $b''<\min\{b,b'\}$ and 
\begin{equation}\label{eq:InteractionNormOfDerviationDelta}
\vvvert \delta^\Phi(\Psi)\vvvert_{b'',Z}\leq C_{\delta}\vvvert \Phi\vvvert_{b'} \: \vvvert  \Psi \vvvert_{b,Z}\,,
\end{equation}
where $C_{\delta}$ depends only on $b,b',b''$ and the dimension $d$.
\end{prop}
\begin{proof}
For $x\in\Gamma$, we wish to estimate
\begin{equation}\label{eq:FirstEstimateDerivationNorm}
	\sum_{\substack{X \in {\mathcal P}_{\rm fin}(\Gamma): \\ x \in X}} 
	 \sum_{\substack{Y,Y' \in {\mathcal P}_{\rm fin}(\Gamma): \\ Y \cap Y' \ne \varnothing, \, Y\cup Y' = X}} \frac{2\|\Phi(Y)\| \|\Psi(Y')\|}{\xi_{b''}(D_Z(X))}.
\end{equation}
There are two possibilities for the second sum, either $x\in Y$  or $x\in Y'\setminus Y$. In the first case we can bound the sum by
\begin{align*}
 \sum_{\substack{Y \in {\mathcal P}_{\rm fin}(\Gamma): \\ x \in Y}} \frac{\|\Phi(Y)\|}{\xi_{b'}(D(Y))}& \sum_{y \in Y} \sum_{\substack{Y' \in {\mathcal P}_{\rm fin}(\Gamma): \\ y \in Y'}} \frac{\|\Psi(Y')\|}{\xi_{b}(D_Z(Y'))} \frac{\xi_{b'}(D(Y)) \xi_{b}(D_Z(Y'))}{\xi_{b''}(D_Z(Y\cup Y'))}
\end{align*}
and in the second case the bound is similar. Now $\textbf{d}(Z,Y\cup Y')\leq \textbf{d}(Z,Y')$. What is more, since $Y,Y'$ are not disjoint, $D(Y\cup Y') \leq D(Y) + D(Y')$ so that monotonicity and superadditivity yield
\begin{equation*}
\frac{\xi_{b'}(D(Y)) \xi_{b}(D_Z(Y'))}{\xi_{b''}(D_Z(Y\cup Y'))}
\leq \frac{\xi_{\min\{b,b'\}}(D(Y) + D(Y') + \textbf{d}(Z,Y'))}{\xi_{b''}(D(Y) + D(Y') + \textbf{d}(Z,Y\cup Y'))}
\leq \xi_{\min\{b,b'\}-b''}(D(Y) + D_Z(Y'))
\end{equation*}
Since $b''<\min\{b,b'\}$, we conclude that
\begin{align*}
		\vvvert\delta^{\Phi}(\Psi)\vvvert_{b'',Z} 
		\leq C \vvvert \Phi\vvvert_{b'} \vvvert \Psi\vvvert_{b,Z}\,,
\end{align*}
where the constant is proportional to $\sup_{Y,Y' \in {\mathcal P}_{\rm fin}(\Gamma) } \left\{|Y|  \xi_{\min\{b,b'\}-b''}(D(Y))\right\}$, as announced.
\end{proof}
\begin{rem}
We note that this is valid for any set $Z$, not necessarily finite. If $Z\in {\mathcal P}_{\rm fin}(\Gamma)$, then by the proposition both $\Psi$ and $\delta^\Phi(\Psi)$ correspond to almost local observables $G^\Psi$ and $G^{\delta^\Phi(\Psi)}$ in ${\mathcal L}(Z)$ and the map $G^\Psi \mapsto \delta^\Phi(G^\Psi):= G^{\delta^\Phi(\Psi)}$ provides the announced extension of $\delta^\Phi$ from ${\mathcal A}_{\mathrm{loc}}$ to the set of almost local observables in ${\mathcal A}$. 
\end{rem}

We conclude this section with a joint corollary of Proposition~\ref{prop:DynamicsIsAlmostLocalObservable} \& \ref{prop:ActionOfDeltas}. They state firstly that the set ${\mathcal L}(Z)$ of almost local observables belongs to the domain of $\tau_t^\Phi$ and $\delta^\Phi$ and secondly that it is invariant under the action of $\tau_t^\Phi$ and $\delta^\Phi$ for any $t$ in a compact interval and arbitrary $Z \in {\mathcal P}_{\rm fin}(\Gamma)$. It follows in particular that Duhamel's formula and its iterates to arbitrary order are well-defined.

\begin{cor}
Let $H\in{\mathcal L}$ with interaction $\Phi\in{\mathcal B}$.
The function $t\mapsto \tau_t^\Phi$ is infinitely often strongly differentiable on the algebra of almost local observables. 
In particular if $O\in{\mathcal L}(Z)$ for some $Z\in{\mathcal P}_{\mathrm{fin}}(\Gamma)$, then Duhamel's formula 
\begin{equation}
\tau^\Phi_t(O) 
= O + \sum_{j=1}^{n-1} \frac{t^j}{j!}\left(\delta^\Phi\right)^j(O)
+ \int_{\Sigma^n_t}\tau^\Phi_{s_n}\left((\delta^\Phi)^n(O)\right)d^n\!s \label{Duahmel}
\end{equation}
is well-defined for any $n\in{\mathbb N}$. We denoted $\Sigma^n_t :=\{0\leq s_1\leq \ldots \leq s_n\leq t\}$ and $d^n\!s = ds_n\ldots ds_1$.
\end{cor}

\section{A product automorphism of lowest order} \label{Sec:Anchor}

With these preliminaries at hand, we now prove the validity of what is sometimes referred to as the symmetric Trotter product formula in the context of an infinite quantum lattice system. Let $H\in{\mathcal L}$ be a local Hamiltonian with interaction $\Phi\in{\mathcal B}$. We assume that
\begin{equation}\label{K's}
H_\Lambda = \sum_{j=1}^k K_{j,\Lambda}
\end{equation}
where $K_{j,\Lambda}\in{\mathcal L}$ are local Hamiltonians with corresponding interactions $\Phi_{j}\in{\mathcal B}$. If $b:=\max\{\tilde b: \Phi\in{\mathcal B}_{\tilde b}\}$, we let 
\begin{equation*}
{\mathcal N} := \vvvert \Phi\vvvert_b
\end{equation*}
and observe that $\Phi_j\in{\mathcal B}_b$ with $\vvvert \Phi_j\vvvert_b\leq{\mathcal N}$ for all $j=1,\ldots,k$. Finally, we let $v$ be the Lieb-Robinson velocity associated with $\Phi$ and note that $v_j\leq v$.

We denote $\tau^j_t = \tau_t^{\Phi_j}$ and $\delta^j = \delta^{\Phi_j}$. Let us first consider the automorphism of ${\mathcal A}$ defined by
\begin{equation}\label{Lowest order symmetric}
\sigma^{(1)}_t(O) := \tau^1_{t/2}\circ\cdots\circ\tau^k_{t/2}\circ\tau^k_{t/2}\circ\cdots\circ\tau^1_{t/2}(O).
\end{equation}
\begin{thm}\label{Thm: anchor}
Let $\mu = \frac{t}{n}$ and let $\pi^{(1)}_{t,n}(O) := \left(\sigma^{(1)}_\mu\right)^n(O)$. Let $Z\in{\mathcal P}_{\mathrm{fin}}(\Gamma)$. Then for any $O \in {\mathcal L}(Z)$,
\begin{equation*}
\left\Vert  \tau^\Phi_t(O) - \pi_{t,n}^{(1)}(O) \right\Vert
\leq C \vvvert O \vvvert \vert Z\vert {\mathcal N}^3 \frac{t^3 \ep{v t }}{n^2}.
\end{equation*}
Here $\vvvert O \vvvert$ denotes the interaction norm of an interaction $\Psi$ such that $O = G^\Psi$.
\end{thm}
\noindent As the proof will show, the constant $C$ is exponentially large in the number of terms $k$, and depends otherwise only on the decay rates of $\Phi,\Psi$ and on the spatial dimension.

Note that in the case $k=2$, namely $H = A+B$, and in finite volume, the product automorphism reduces to the adjoint action of $\left(\ep{i\frac{t}{2n}B}\ep{i\frac{t}{n}A}\ep{i\frac{t}{2n}B}\right)^n$, which is indeed well-known to converge to $\mathrm{Ad}_{\ep{i t (A+B)}}$
as $n\to\infty$. 
\begin{proof}[Proof of Theorem~\ref{Thm: anchor}]
We first decompose the time interval $[0,t]$ in $n$ subintervals of width $\mu = \frac{t}{n}$ to get the following telescopic sum:
\begin{equation}\label{Telescopic}
 \tau^\Phi_t(O) - \pi_{t,n}^{(1)}(O)
= \sum _{j = 0}^{n-1} \left(\sigma^{(1)}_\mu\right)^j\left((\tau^\Phi_\mu -\sigma^{(1)}_\mu)\left( (\tau^\Phi_{\mu})^{n-j-1}(O) \right)\right).
\end{equation}
For any almost local observable $\widetilde O\in {\mathcal L}(Z)$, we see that
\begin{equation*}
\frac{d}{ds}\left.\left(\sigma^{(1)}_{s}\circ\tau^\Phi_{-s}\right)(\widetilde O)\right\vert_{s=0}
= (2\sum_{j=1}^k \frac{\delta^j}{2} - \delta^\Phi)(\widetilde O) = 0
\end{equation*}
by~(\ref{Lowest order symmetric}) and~(\ref{K's}). Similarly, but with a little more algebra, 
\begin{align*}
\frac{d^2}{ds^2}\left.\left(\sigma^{(1)}_{s}\circ\tau^\Phi_{-s}\right)(\widetilde O)\right\vert_{s=0}
&= \bigg(\frac{1}{4}\sum_{j=1}^k \Big\{\sum_{l=1}^j\delta^l \delta^j + \sum_{l=j+1}^k\delta^j \delta^l + \delta^j \sum_{l=1}^k \delta^l-2\delta^j\delta^\Phi \\
&\qquad\qquad\quad+ \sum_{l=1}^k \delta^l \delta^j + \sum_{l=j+1}^k\delta^l \delta^j + \sum_{l=1}^j\delta^j \delta^l  -2\delta^j\delta^\Phi\Big\} \\
&\qquad-2\sum_{j=1}^k \frac{1}{2}\delta^j\delta^\Phi
+\delta^\Phi\delta^\Phi
\bigg)(\widetilde O )
\end{align*}
Writing $\sum_{j=1}^k\sum_{l=1}^j \delta^l \delta^j = (\delta^\Phi)^2 - \sum_{j=1}^k\sum_{l=j+1}^k \delta^l \delta^j$ and proceeding similarly for the second-to-last term of the second line, we conclude that this derivative vanishes again by $\delta^\Phi = \sum_{j=1}^k\delta^j$. Thus,
\begin{align}%
\tau^\Phi_\mu(\widetilde O) -\sigma^{(1)}_\mu(\widetilde O)
&=-\left.\left(\sigma^{(1)}_{s}\circ\tau^\Phi_{-s}\right)(\tau^\Phi_\mu(\widetilde O))\right\vert_{s=0}^{s=\mu} \nonumber \\
&= - \int_{\Sigma_\mu^3} \frac{d^3}{ds_3^3}\left(\sigma^{(1)}_{s_3}\circ\tau^\Phi_{-s_3}\right)(\tau^\Phi_\mu(\widetilde O)) d^3s.\label{Infinitesimal steps}
\end{align}

Distributing the three derivatives across the $2k$ factors of $\sigma^{(1)}_{s_3}\circ\tau^\Phi_{-s_3}$ 
inserts three derivations to the product $\tau^1_{s_3/2}\circ\cdots\circ\tau^k_{s_3}\circ\cdots\circ\tau^1_{s_3/2}\circ\tau^\Phi_{-s_3}(\tau^\Phi_\mu(\widetilde O))$. By Propositions~\ref{prop:DynamicsIsAlmostLocalObservable} and~\ref{prop:ActionOfDeltas} all terms are well-defined and belong to ${\mathcal L}(Z)$ since $\tau^\Phi_\mu(\widetilde O) \in{\mathcal L}(Z)$. By the same propositions, the interaction norm of each of them is bounded by 
$C_{\tau}^{2k+1}(C_{\delta} {\mathcal N})^3 \ep{v((k+1)s_3 + \mu)}\vvvert\widetilde O\vvvert$, where $(k+1)s_3 + \mu$ is the total time (in absolute value) involved in~$\left(\sigma^{(1)}_{s_3}\circ\tau^\Phi_{-s_3}\right)(\tau^\Phi_\mu(\widetilde O))$. Finally, we recall from~\eqref{Telescopic} that $\widetilde O = (\tau^\Phi_{\mu})^{n-j-1}(O) = \tau^\Phi_{(n-j-1)\mu}(O)$ with $O\in{\mathcal L}(Z)$, so that $\vvvert\widetilde O\vvvert\leq C_{\tau}\ep{v(n-j-1)\mu}\vvvert O\vvvert$ and hence 
\begin{align*}
\left\vvvert \left(\tau^\Phi_\mu -\sigma^{(1)}_\mu\right)\left( (\tau^\Phi_{\mu})^{n-j-1}(O) \right) \right\vvvert
&\leq \vvvert O \vvvert (2k)^3C_{\delta}^3 C_{\tau}^{2k+2}{\mathcal N}^3 \ep{v (n-j)\mu}\int_{\Sigma_\mu^3} \ep{v(k+1)s_3} d^3s \\
&\leq C \vvvert O \vvvert \ep{v (n-j)\mu}\frac{({\mathcal N}\mu)^3}{3!}\ep{v(k+1)\mu}.
\end{align*}
Since $\sigma^{(1)}_\mu$ preserves the operator norm, each term of~(\ref{Telescopic}) is bounded by
\begin{equation*}
\left\Vert (\sigma^{(1)}_\mu)^j((\tau^\Phi_\mu -\sigma^{(1)}_\mu) ((\tau^\Phi_{\mu})^{n-j-1}(O) ))\right\Vert
\leq C \vvvert O \vvvert \vert Z\vert \ep{v t }{\mathcal N}^3 \mu^3
\end{equation*}
by Lemma~\ref{lemma: finite Z interaction}, where we used that $(n-j)\mu\leq n\mu =t$ and absorbed the remaining exponential in $k$ in the constant. This estimate being uniform across the $n$ terms of~(\ref{Telescopic}), we immediately conclude with the claim of the theorem.
\end{proof}

\section{Arbitrary order}\label{sec:ArbOrder}

The symmetric Trotter formula discussed in the previous section has an error of order $n^{-2}$. As pioneered by Suzuki in, e.g.~\cite{suzuki1991general}, a recursive construction can be build upon it to generate higher order product formulae. We now show that they too extend to the infinite volume setting. 

\subsection{Time reversal}

Since $\sigma^{(1)}_\mu$ is a composition of automorphisms, it is an automorphism, but the fact that the individual factors do not commute with each other breaks the group property of the fundamental time evolution $\tau^\Phi_t\circ\tau^\Phi_{s} = \tau^\Phi_{t+s}$. However, the specific `symmetric' form of~(\ref{Lowest order symmetric}) implies that
\begin{equation*}
\sigma^{(1)}_{-\mu}\circ\sigma^{(1)}_\mu = \mathrm{id}.
\end{equation*}
A product automorphism having this property shall be called time-reversal symmetric. Theorem~\ref{Thm: anchor} shows that, despite its label $(1)$, the corresponding product automorphism $\pi_{t,n}^{(1)}$ is in fact a second order approximation of $\tau_t^\Phi$. This improvement from any odd order to the next even order is in fact general for time-reversal symmetric product approximations.
\begin{prop}\label{prop: TRI}
Let $m\in{\mathbb N}$ and let $\sigma^{(2m-1)}_\mu$ be an $(2m-1)$-th order product approximation of $\tau^\Phi_\mu$ in the sense that
\begin{equation}\label{lower orders}
\left.\frac{d^j}{d\mu^j} \left(\tau^\Phi_\mu(\widetilde O) - \sigma^{(2m-1)}_\mu(\widetilde O)\right)\right\vert_{\mu = 0} = 0\qquad(j\in\{0,\ldots,2m-1\})
\end{equation}
for any $\widetilde O \in {\mathcal L}(Y)$. If $\sigma^{(2m-1)}_\mu$ is time-reversal symmetric,
\begin{equation*}
\sigma^{(2m-1)}_{-\mu}\circ\sigma^{(2m-1)}_\mu = \mathrm{id},
\end{equation*}
then it is a $(2m)$-th order approximation of $\tau^\Phi_\mu$.
\end{prop}
\begin{proof}
The identity
\begin{equation*}
\widetilde O = \tau^\Phi_{-\mu}\circ(\tau^\Phi_\mu-\sigma^{(2m-1)}_\mu)(\widetilde O) + (\tau^\Phi_{-\mu}-\sigma^{(2m-1)}_{-\mu})\circ\sigma^{(2m-1)}_\mu(\widetilde O) + \sigma^{(2m-1)}_{-\mu}\circ\sigma^{(2m-1)}_\mu(\widetilde O)
\end{equation*}
and time-reversal symmetry imply that
\begin{equation*}
\tau^\Phi_{-\mu}\circ(\tau^\Phi_\mu-\sigma^{(2m-1)}_\mu)(\widetilde O) + (\tau^\Phi_{-\mu}-\sigma^{(2m-1)}_{-\mu})\circ\sigma^{(2m-1)}_\mu(\widetilde O) = 0.
\end{equation*}
The derivative of order $2m$ of this equation at $\mu=0$ reduces by~(\ref{lower orders}) to
\begin{equation*}
\left.\frac{d^{2m}}{d\mu^{2m}}\left(\tau_\mu^\Phi-\sigma^{(2m-1)}_\mu\right)(\widetilde O)\right\vert_{\mu=0}
+ \left.\frac{d^{2m}}{d\mu^{2m}}\left(\tau_{-\mu}^\Phi-\sigma^{(2m-1)}_{-\mu}\right)(\widetilde O)\right\vert_{\mu=0} =0,
\end{equation*}
which concludes the proof since the two derivatives of even order are equal.
\end{proof}

\subsection{Suzuki's Ansatz} 

We now recall Suzuki's inductive construction~\cite{suzuki1991general} of higher order product formulae, translated in the present language of automorphisms. Since Section~\ref{Sec:Anchor} provided a time reversal symmetric approximation of order $2$, we shall use it to anchor the induction. 

Let $\sigma_\mu^{(2m)}$ be a time reversal symmetric $(2m)$-th order product approximation of $\tau_\mu^\Phi$ in the sense of~(\ref{lower orders}). A higher order approximation can be constructed as follows. Let $r\geq 3$ be an odd integer, $r = 2\ell + 1$ and let
\begin{equation}\label{s coefficients}
s_{m} := \frac{1}{(r-1)-(r-1)^{\frac{1}{2m+1}}}.
\end{equation}
Then 
\begin{equation}\label{Higher order construction}
\sigma_\mu^{(2m+1)} :=\left(\sigma_{s_{m}\mu}^{(2m)}\right)^\ell\circ\sigma_{(1-(r-1)s_{m})\mu}^{(2m)}\circ\left(\sigma_{s_{m}\mu}^{(2m)}\right)^\ell,
\qquad \sigma_\mu^{(2m+2)}:=\sigma_\mu^{(2m+1)}.
\end{equation}
We immediately point out firstly that $2\ell s_m + (1-(r-1)s_m) = 1$ and secondly that $ -1+(r-1)s_m =s_m(r-1)^{\frac{1}{2m+1}}$, and so
\begin{equation}\label{scalings}
(r-1)s_m^{2m+1} + \left(1-(r-1)s_{m}\right)^{2m+1} = 0.
\end{equation}
\begin{thm}\label{Thm: TRI-Trotter}
Let $\sigma_\mu^{(1)}$ be defined by~(\ref{Lowest order symmetric}) and let $\sigma_\mu^{(2)}:=\sigma_\mu^{(1)}$. For $m\in{\mathbb N}$, let $\sigma_\mu^{(m)}$ be defined recursively by~(\ref{Higher order construction}). Then for all $m\geq 1$, $\sigma_\mu^{(m)}$ is time-reversal symmetric. In addition, let
\begin{equation*}
\pi^{(m)}_{t,n} := \left(\sigma^{(m)}_\mu\right)^n\qquad\left(\mu = \frac{t}{n}\right).
\end{equation*}
Then, for any $O \in {\mathcal L}(Z)$, 
\begin{equation}\label{main formula}
\left\Vert  \tau^\Phi_t(O) - \pi_{t,n}^{(m)}(O) \right\Vert
\leq C \vvvert O\vvvert \vert Z \vert \frac{t^{\alpha+1}\ep{vt}}{n^{\alpha}},
\end{equation}
with 
$
	\alpha = \begin{cases}\;\;\,m \quad\,, \textrm{ if $m$ is even} \\
								m+1 \,, \textrm{ if $m$ is odd} 
			\end{cases}
$.
\end{thm}
\begin{proof}
The symmetry for all $m$ is immediate by~(\ref{Higher order construction}) since $\sigma_\mu^{(1)}$ is symmetric. The estimate holds by Theorem~\ref{Thm: anchor} for $m=1,2$, so we proceed by induction. We assume that $\sigma^{(2m)}_\mu$ is a $(2m)$-th order approximation of $\tau_\mu^\Phi$ and that the claim of the theorem holds for $2m$. We write as in \eqref{Telescopic}
\begin{equation}\label{eq:Telescopic-General}
 \tau^\Phi_t(O) - \pi_{t,n}^{(2m+1)}(O)
= \sum _{j = 0}^{n-1} \left(\sigma^{(2m+1)}_\mu\right)^j\left((\tau^\Phi_\mu -\sigma^{(2m+1)}_\mu)\left( (\tau^\Phi_{\mu})^{n-j-1}(O) \right)\right),
\end{equation}
and proceed with an estimate on $(\tau^\Phi_\mu -\sigma^{(2m+1)}_\mu)(\widetilde O)$ for an almost local observable $\widetilde O\in{\mathcal L}(Y)$. Here again, we decompose the interval $[0,\mu]$ into $r = 2\ell + 1$ intervals according to~(\ref{Higher order construction}) and obtain
\begin{align*}
\tau^\Phi_\mu(\widetilde O) -\sigma_\mu^{(2m+1)}(\widetilde O)
&= \sum_{j=0}^{\ell-1} \left(\sigma_{s_{m}\mu}^{(2m)}\right)^j \circ \left(\tau_{s_{m}\mu}^\Phi - \sigma_{s_{m}\mu}^{(2m)}\right) \circ \left(\tau_{((2\ell-j-1)s_{m} + \tilde s_m)\mu}^{\Phi}\right) (\widetilde O) \\
&\quad + \left(\sigma_{s_{m}\mu}^{(2m)}\right)^\ell \circ \left(\tau^\Phi_{\tilde s_m\mu} - \sigma_{\tilde s_m \mu}^{(2m)}\right) \circ \left(\tau_{\ell s_{m}\mu}^\Phi\right) (\widetilde O)  \\
&\quad + \sum_{j=1}^{\ell} \left(\sigma_{s_{m}\mu}^{(2m)}\right)^{\ell} \circ \sigma_{\tilde s_m\mu}^{(2m)}\circ  \left(\sigma_{s_{m}\mu}^{(2m)}\right)^{j-1} \circ \left(\tau_{s_{m}\mu}^\Phi - \sigma_{s_{m}\mu}^{(2m)}\right) \circ \left(\tau_{(\ell-j)s_{m}\mu}^\Phi\right) (\widetilde O)
\end{align*}
where we denoted $\tilde s_m = 1-(r-1)s_{m}$. By the induction hypothesis, $\left.\frac{d^{j}}{d\mu^{j}}\left(\tau^\Phi_\mu(\widetilde O) - \sigma_\mu^{(2m)}(\widetilde O)\right)\right\vert_{\mu=0}=0$ for all $j=0,\ldots,2m$. This and the above identity imply first of all that the same holds with $\sigma_\mu^{(2m+1)}$ instead of $\sigma_\mu^{(2m)}$, and secondly that
\begin{equation*}
\left.\frac{d^{2m+1}}{d\mu^{2m+1}}\left(\tau^\Phi_\mu(\widetilde O)-\sigma_\mu^{(2m+1)}(\widetilde O)\right)\right\vert_{\mu=0}
= \left(2\ell s_m^{2m+1} + \tilde s_m^{2m+1}\right)\left.\frac{d^{2m+1}}{d\nu^{2m+1}}\left(\tau^\Phi_{\nu}(\widetilde O)-\sigma_{\nu}^{(2m)}(\widetilde O)\right)\right\vert_{\nu=0}.
\end{equation*}
Since $2\ell = r-1$, this vanishes by~(\ref{scalings}), so that $\sigma_\mu^{(2m+1)}$ is a $(2m+1)$-th order approximation. Since $\sigma^{(2m+1)}_\mu$ is time reversal symmetric we get from Proposition~\ref{prop: TRI} that the $(2m+2)$-th derivative similarly vanishes at $\mu=0$. Hence, 
\begin{equation*}
\tau^\Phi_\mu(\widetilde O)-\sigma_\mu^{(2m+1)}(\widetilde O)
=\int_{\Sigma^{2m+3}_\mu}\frac{d^{2m+3}}{ds^{2m+3}}\left(\sigma_{s}^{(2m+1)}\circ\tau^\Phi_{-s}\right)\left(\tau^\Phi_\mu(\widetilde O) \right) d^{2m+3}\!s,
\end{equation*}
and we can proceed as in the proof of~Theorem~\ref{Thm: anchor}. 

There are $r^m(2k-2)+1$ factors in $\sigma_s^{(2m+1)}$ and hence a total of $2^{2m+3}(r^m(k-1)+1)^{2m+3}$ 
terms in the derivative, each of them involving a combination of $2m+3$ derivations from $\{\delta^\Phi\}\cup\{\delta^{j}:j=1,\ldots k\}$. If $s^{(m)}$ is the total time (in absolute value) involved in $\sigma_s^{(m)}$ (for example, $s^{(1)} = s^{(2)} = ks$), then $s^{(2m+1)} = (r-1)s_{m}s^{(2m)} + \left\vert (1-(r-1)s_m)\right\vert s^{(2m)} = (2(r-1)s_m -1)s^{(2m)}$ since $1-(r-1)s_m = -s_m(r-1)^{\frac{1}{2m+1}}<0$. Hence,
\begin{equation*}
s^{(2m+1)} = \bigg(\prod_{j=1}^m(2(r-1)s_j -1)\bigg) ks.
\end{equation*}
Setting $\widetilde O = (\tau^\Phi_{\mu})^{n-j-1}(O)$, we conclude by Propositions~\ref{prop:DynamicsIsAlmostLocalObservable} and~\ref{prop:ActionOfDeltas} that
\begin{multline*}
\left\vvvert\frac{d^{2m+3}}{ds^{2m+3}}\left(\sigma_{s}^{(2m+1)}\circ\tau^\Phi_{-s}\right)\left(\tau^\Phi_{(n-j)\mu}(O) \right)\right\vvvert \\
\leq \vvvert  O\vvvert\left(2(r^m(k-1)+1)  C_\delta \right)^{2m+3}  C_\tau^{\left(r^m(2k-2)+3\right)} {\mathcal N}^{2m+3} \ep{v \mu(n-j)} \ep{v (s^{(2m +1)} +s)}.
\end{multline*}
Integrating this over the simplex $\Sigma^{2m+3}_\mu$ and gathering all constants yields
\begin{equation*}
\left\Vert 
\left(\sigma^{(2m+1)}_\mu\right)^j\left((\tau^\Phi_\mu -\sigma^{(2m+1)}_\mu)\left( (\tau^\Phi_{\mu})^{n-j-1}(O) \right)\right)
 \right\Vert
\leq C \vvvert  O\vvvert \vert Z\vert \mu^{2m+3} \ep{v t}.
\end{equation*}
Since there are $n$ such terms in~(\ref{eq:Telescopic-General}), we have now proved that~(\ref{main formula}) holds for $2m+1$ and therefore also for $2m+2$ by the definition~(\ref{Higher order construction}) of $\sigma_\mu^{(2m+2)}$, concluding the induction.
\end{proof}
\begin{rem}
The theorem should not be misinterpreted as an invitation to take a limit in $m$. Rather, it provides for each fixed $m$ a formula that scales as $n^{-m}$ as $n\to\infty$. As can be read from the proof, the constant $C$ scales as $\frac{r^{m^2}\ep{r^m}}{m!}$, underlying the importance of picking a small possible $r$, namely $r=3,5$. 
\end{rem}
As is clear in the proof (see also~\cite{suzuki1991general}), the choice of time scalings by $s_m$, respectively $1-(r-1)s_{m}$ is largely arbitrary. The claim of the theorem would continue to hold for any $2\ell + 1$ real coefficients $\{p_{m,j}:j\in\{1,\ldots,r\}\}$ in~(\ref{Higher order construction}), provided $p_{m,j} = p_{m,r+1-j}$ for all $j=1,\ldots,\ell$ as well as
\begin{equation}\label{pmj Conds}
\sum_{j=1}^r p_{m,j} = 1,\qquad 
\mathrm{and}\qquad\sum_{j = 1}^r p_{m,j}^{2m+1} = 0.
\end{equation}
Reality ensures that $\pi_{t,n}^{(m)}$ are automorphisms and the symmetric choice of coefficient on either side of $p_{m,\ell+1}$ is for time reversal symmetry. Clearly, there is no non-trivial positive solution of these equations, and (\ref{Higher order construction}) indeed has $p_{m,\ell+1} = -s_m(r-1)^{\frac{1}{2m+1}}<0$, as already pointed out. The appearance of such a negative time evolutions to improve the order of the approximation is reminiscent of the decomposition proposed in~\cite{haah2021quantum}.

\begin{rem}\label{rem:ScalingSmaller1For5}
If $r=3$, then $\vert s_m \vert, \vert 1-(r-1)s_m \vert>1$ with  $\lim_{m\to\infty} s_m =1, \lim_{m\to\infty} (1-(r-1)s_m ) =-1$ and so the individual time intervals in the product scale as $\frac{t}{n}$, independently of $m$ for large $m$. On the other hand, if $r=5,7,\ldots$, then $\vert s_m \vert, \vert 1-(r-1)s_m \vert<1$ with $\lim_{m\to\infty}s_m = \frac{1}{r-2}$ and  $\lim_{m\to\infty} (1-(r-1)s_m ) =-\frac{1}{r-2}$. Hence the individual time intervals in the product scale as $\frac{1}{(r-2)^m}\frac{t}{n}$. The inductive construction and the appearance of negative signs yield a fractal path in the time domain. This behaviour --- already noted in~\cite{suzuki1991general} ---  is exhibited in Figure~\ref{fig:Nice plots}. 

\begin{figure}
\begin{minipage}{0.45\linewidth}
\includegraphics[width=\textwidth]{Path_m9r3}
\end{minipage}\hfill
\begin{minipage}{0.45\linewidth}
\includegraphics[width=\textwidth]{Path_m9r5}
\end{minipage}
\caption{
The discrete time steps in $\sigma_\mu^{(m)}$ for $r=3$ on the left and $r=5$ on the right, for the same order $m=9$ and in the case $k=2$. There is an order of magnitude difference between the number of terms involved, with the ratio of the number of terms being $\frac{5^42+1}{3^42+1}\simeq 7.7$.}
\label{fig:Nice plots}
\end{figure}
\end{rem}


\section{Quantum simulation: Decomposition in commuting Hamiltonians} \label{sec:QuantumSimulation}

So far, the results are completely general in the sense that they do not require any assumption on the local Hamiltonians $K_{j,\Lambda}\in{\mathcal L}$ beyond their locality. In concrete applications however,
the choice of decomposition of $H$ is determined by the requirement that each $K_{j,\Lambda}$ is a sum of terms acting on spatially disjoint subsets of the lattice and hence mutually commuting. In the simple example of a one-dimensional lattice with nearest-neighbour interaction, namely $\Phi(X) = 0$ if $X\neq\{x,x+1\}$ for some $x\in{\mathbb Z}$, one would choose $\Phi_1,\Phi_2$ to be supported on pairs of neighbouring sites $\{2x,2x+1\}$, respectively $\{2x+1,2x+2\}$. 
Each dynamics $\tau_t^{\Phi_1}, \tau_t^{\Phi_2}$ is then strictly local and corresponds to the action of quantum gates in parallel, as illustrated in Figure~\ref{fig:FDQC}.
\begin{figure}
\begin{center}
\includegraphics[width=0.75\textwidth]{FDQC.pdf}
\end{center}
\caption{A cartoon of the action of a finite depth quantum circuit on a 1-dimensional spin chain, in blue. Since each layer is generated by a commuting Hamiltonian, it is a product of commuting unitaries. Propagation is induced by the alternating action of overlapping layers.}\label{fig:FDQC}
\end{figure}
In this setting where a product formula is referred to as a \emph{finite depth unitary quantum circuit}, Theorem~\ref{Thm: TRI-Trotter} provides a quantitative bound on the error in the approximation of the full dynamics $\tau_t$ by a circuit. The number of factors in the product automorphism is referred to as the \emph{depth} of the circuit. 
\begin{cor}
Let $\Phi$ be a finite range interaction, namely $\Phi(X) = 0$ if $D(X) >R$ for a fixed $R>0$. Let $Z\in{\mathcal P}_{\mathrm{fin}}(\Gamma)$ and $O\in{\mathcal L}(Z)$. For any $m\in{\mathbb N}$, there is a finite depth unitary quantum circuit of depth ${\mathcal O}(\epsilon^{-\frac{1}{m}})$ approximating $\tau_t^\Phi(O)$ within error $\epsilon$, as $\epsilon\to0$.
\end{cor}
\begin{proof}
The depth $h$ of the circuit $\pi_{t,n}^{(m)}$ is proportional to $n$. Hence, imposing that the bound~(\ref{main formula}) is less than $\epsilon$,
\begin{equation}\label{error scaling}
C \vvvert O\vvvert \vert Z \vert \frac{t^{\alpha+1}\ep{vt}}{h^{\alpha}}\leq \epsilon,
\end{equation}
yields the claim, since $\alpha\geq m$.
\end{proof}
\noindent We further point out that, as should be expected, the depth of the circuit diverges as $\vert Z\vert^{\frac{1}{m}}$ with the volume of the support $Z$. The depth is furthermore exponential in time, but with a rate $\frac{v}{m}$ that is smaller for a higher order product automorphism. The depth of the circuit proposed in~\cite{haah2021quantum} scales as ${\mathcal O}\left(t \: \mathrm{polylog}(N t\epsilon^{-1})\right)$, namely better in the total time. However, this is divergent in the size of the system $N$, and uses ${\mathcal O}\left(\mathrm{polylog}(Nt\epsilon^{-1})\right)$ additional ancilla qubits to achieve this.

In the context of finite range interactions and for a strictly local observable, one may wish to decompose the telescopic sum~(\ref{Telescopic}) rather as
\begin{equation}\label{FDQC:Telescopic}
 \tau^\Phi_t(O) - \pi_{t,n}^{(1)}(O)
= \sum _{j = 0}^{n-1} \left(\tau^\Phi_{\mu}\right)^j\left((\tau^\Phi_\mu -\sigma^{(1)}_\mu)\left( (\sigma^{(1)}_\mu)^{n-j-1}(O) \right)\right),
\end{equation}
since $(\sigma^{(1)}_\mu)^{n-j-1}(O)$ is strictly local. If the $K_{j,\Lambda}$'s are all commuting Hamiltonians, all factors in $(\sigma^{(1)}_\mu)^{n-j-1}$ induce no propagation beyond the range~$R$ of the interactions so that $(\sigma^{(1)}_\mu)^{n-j-1}(O)$ is strictly supported in $Z^{((n-j-1)(2k-1)R)}$ and of operator norm $O$. However, while the action of all derivations $\delta^\Phi,\delta^j$ is again strictly local, increasing the support by $R$, the bound $\Vert \delta^\Upsilon (\widetilde O)\Vert \leq C \vvvert \Upsilon\vvvert_{b,Y} \vert Y\vert \Vert \widetilde O\Vert$ valid for an observable $\widetilde O\in{\mathcal A}_Y$ yields an estimate
\begin{equation*}
C \vert Z \vert^3 n^{3d} \mu^3\Vert O\Vert,
\end{equation*}
where $d$ is the spatial dimension, for every term of~(\ref{FDQC:Telescopic}). One would therefore obtain $\Vert  \tau^\Phi_t(O) - \pi_{t,n}^{(1)}(O)\Vert = {\mathcal O}(n^{3d-2})$, emphasizing the need for a careful use of the Lieb-Robinson bound on the very short time intervals of width ${\mathcal O}(n^{-1})$. This also shows that the physical propagation is in fact much slower than what can be read off from the diagram in Figure~\ref{fig:FDQC}. 

\subsection{Long range interactions} While long range interactions pose no issue for our main theorem, Theorem~\ref{Thm: TRI-Trotter}, they cannot be decomposed as a finite depth unitary quantum circuit as just described. A necessary intermediate step is a truncation to finite range interaction. We now comment the error induced by neglecting the interactions between far enough lattice points.

Let $\Phi\in{\mathcal B}_b$ and $\Phi_R\in{\mathcal B}_b$ be defined by
\begin{equation*}
\Phi_R(X) = \begin{cases}
\Phi(X) & \text{if }D(X)\leq R \\ 0 & \text{otherwise}
\end{cases}.
\end{equation*}
Then, for any $b'<b$,
\begin{equation*}
\sum_{X\ni x}\frac{\Vert \Phi_R(X) - \Phi(X) \Vert}{\xi_{b'}(D(X))}
= \sum_{X\ni x: D(X)>R}\frac{\Vert \Phi(X) \Vert}{\xi_{b}(D(X))}\xi_{b-b'}(D(X))
\end{equation*}
so that $\vvvert \Phi_R - \Phi \vvvert_{b'} \leq \xi_{b-b'}(R+1) \vvvert \Phi\vvvert_b$. Moreover, 
\begin{equation*}
\tau_t^{\Phi_R}(O) - \tau_t^{\Phi}(O) = \int_0^t \tau_s^{\Phi_R}\left(\delta^{\Phi_R} - \delta^{\Phi}\right)(\tau^{\Phi}_{t-s}(O))ds
\end{equation*}
so that if $O\in{\mathcal L}(Z)$, 
\begin{equation*}
\left \vvvert \tau_t^{\Phi_R}(O) - \tau_t^{\Phi}(O)\right\vvvert
\leq C t \vvvert \Phi_R - \Phi \vvvert \vvvert O\vvvert
\leq C t\vvvert \Phi\vvvert\Inorm O\vvvert \xi_{\tilde b}(R+1)
\end{equation*}
for a $\tilde b>0$, since $\delta^{\Phi_R} - \delta^{\Phi} = \delta^{\Phi_R - \Phi}$. In other words, the error associated with the truncation of the interaction is superpolynomially small in the range $R$. In order to achieve an error ${\mathcal O}(\epsilon)$, the range must be chosen as $R = {\mathcal O}(\xi_{\tilde b}^{-1}(\epsilon)) = {\mathcal O}((\log{\epsilon^{-1}})^{\frac{1}{p}})$. A commuting decomposition for an interaction of range $R$ requires of the order of $k = R^d$ terms in $d$ dimensions, and since the constant in the Trotter error~(\ref{main formula}) is exponentially large in the number of terms, (\ref{error scaling}) and the lower bound $\exp\left((\log{\epsilon^{-1}})^{\frac{d}{p}}\right)\geq \epsilon^{-1}$ yield a circuit depth $h = {\mathcal O}(\epsilon^{-\frac{2}{m}})$.

\section*{Acknowledgements}

\noindent The authors are supported by NSERC of Canada. S.B.~would like to thank Marcel Schaub for introducing the specific interaction norm used in this work.

\bibliographystyle{unsrt}
