\section{The Hamiltonian}\label{The Hamiltonian}
In this section we define the Hamiltonians we will be considering in this article\footnote{This
will be extended to include multiple copies of the bosonic bath Hamiltonian in section
\ref{Multiple baths extension}.}. We consider a quantum system coupled to a bosonic bath. The
Hilbert space $\mathcal{S}$ of the system carries a ``free'' Hamiltonian $\hat H_{\mathcal{S}}$,
while the bosonic bath is described by the Fock space $\mathcal{B}:=\mathrm{\Gamma}(\mathfrak{h})$ over
the mode Hilbert space $\mathfrak{h}$, with the free Hamiltonian $d\hat{\mathrm{\Gamma}}(G)$ arising from
the second quantization of the positive self-adjoint operator $G$ on $\mathfrak{h}$. The coupling 
is via one mode $h\in\mathfrak{h}$ and a bounded in operator norm self-adjoint system operator
$\hat A_\mathcal{S}$, so that altogether we have on $\mathcal{S}\otimes\mathcal{B}$ the Hamiltonian
\begin{equation}\label{the bathtub}
    \hat H=\hat H_{\mathcal{S}}\otimes{\mathbbm{1}}_\mathcal{B}+{\mathbbm{1}}_{\mathcal{S}}\otimes d\hat{\mathrm{\Gamma}}(G)+\hat A_\mathcal{S}\otimes\hat\Phi(h),
\end{equation}
where $\hat\Phi$ denotes the usual field operator. $G$ is a function $g$ of momentum variables
$k$ (known as the \textit{dispersion relation}), and can be expressed in terms of bath creation and
annihilation operators as
\begin{equation}\label{eq:bath def}
    d\hat{\mathrm{\Gamma}}(G)=\int d{k}\ g({k})a_{{k}}^\dag a_{{k}}
\end{equation}
and
\begin{equation}\label{eq:int temr}
    \hat\Phi(h)=\int d{k}\ h({k})\bigl(a_{{k}}^\dag+a_{{k}}\bigr).
\end{equation}
The form Eq. (\ref{the bathtub}) is often referred to as a Pauli-Fierz Hamiltonian \cite{dere99},
specialized in our case by allowing one interaction term. We note that although $\hat A_\mathcal{S}$
is bounded, there are no constraints on the system Hamiltonian $\hat H_\mathcal{S}$, it can be
unbounded or otherwise and comprise of bosons, fermions, spins etc. Definining the \textit{free} Hamiltonian as $\hat H_0=\hat H_{\mathcal{S}}\otimes{\mathbbm{1}}_\mathcal{B}+{\mathbbm{1}}_{\mathcal{S}}\otimes d\hat{\mathrm{\Gamma}}(G) $, Eq. (\ref{the bathtub}) is well
defined on $\mathcal{D}(\hat H)=\mathcal{D}(\hat H_0)$ if \cite{Reed1975
\begin{equation}
\int\frac{h^2(k)}{g(k)}dk<\infty.
\end{equation}
It includes the spin-boson-Model achieved by letting $\mathcal{S}=\mathbb{C}^2$, $\hat H_\mathcal{S}=\mathcal{\alpha}\hat\sigma_z$, and $\hat A_\mathcal{S}=\hat\sigma_x$, where
$\mathcal{\alpha}$ is a positive constant and $\hat\sigma_x,$ $\hat\sigma_z$ are the Pauli
matrices.
Let $E(dx)$ be the projection valued spectral measure of $G$. Then we form the scalar measure
\begin{equation}\label{mu}
    \mu_0(d\omega)=\langle h,E(d\omega)h\rangle\ .
\end{equation}
The measure is absolutely continuous with respect to the Lebesgue measure, and we can write
$\mu_0(dx)= \pi^{-1} J(x)dx$ with the \textit{spectral density}\footnote{here the factor of
$1/\pi$ is for conformity with \cite{Weiss2001}} $J$.
When $g$ is monotone, it is defined as
\begin{equation}\label{SD_eq}
    J(\omega)=\pi h\left(g^{-1}(\omega)\right)^2\ \left|\frac{dg^{-1}(\omega)}{d\omega}\right|,
\end{equation}
where $g^{-1}$ is the inverse of $g$. If $g$ is not monotone, we would additionally have a sum
over inverse images $g^{-1}(\{\omega\}),$ and when the momentum variable has more dimensions, we
would also have an integral over the inverse image. The minimally closed interval containing its
support is $[\omega_{min},\omega_{max}]$ with $\omega_{min}:=\textup{inf}\; g\geq 0,$ $\omega_{max}:=\textup{sup}\; g$.
The case $\omega_{min}=0$ is called \textit{massless} where as $\omega_{min}>0$ is known as \textit{massive}.
We will be dealing with the case $\omega_{max}<\infty$, hence $\mu_0(dx)$ is determinate with all
moments finite (see e.g. \cite{chiara}). For later purposes, it will be convenient to recall 
that Hamiltonian Eq. \eqref{the bathtub} is isomorphic to the so-called \textit{standard form} defined 
on page 5 of \cite{misc14}. We summarise this standard form here, for the  convenience of the reader.
We observe that the bath in Hamiltonian Eq. \eqref{the bathtub} is fully specified by the triple 
$(\mathfrak{h},G,h)$. Suppose now we have another environment system $(\tilde{\mathfrak{h}},\tilde G,\tilde h)$ 
and a unitary operator $U : \mathfrak{h} \rightarrow \tilde{\mathfrak{h}}$
such that $U h = \tilde h$ and $U G U^\dag = \tilde G$. Under such an isomorphism all the details of 
system-environment are mapped into each other. This is formally done by the unitary operator 
$\hat{\mathrm{\Gamma}} (U) : \mathrm{\Gamma}(\tilde{\mathfrak h}) \rightarrow \mathrm{\Gamma}(\tilde{\mathfrak h})$. 
Note that the scalar measure Eq. \eqref{mu} is invariant under this isomorphisum $U$. In fact, 
it completely determines the triple $(\mathfrak{h},G,h)$ up to an isomorphisum. Indeed, we are permitted to set
\begin{align}\label{eq:std form}
\begin{split}
\tilde{\mathfrak{h}}&=L^2({\mathbbm{R}}^+,\mu),\\
(\tilde G\psi)(x)&=x \psi(x)\;\; \forall\, \psi\in \tilde{\mathfrak{h}},\\
\tilde{h}(x)&=1.
\end{split}
\end{align}
The unitary operator defining the isomorphism is given by $\left( U \mathrm{e}^{\mathrm{i} t G}h\right)(x)=\mathrm{e}^{\mathrm{i} t x}.$ The triple $(\tilde{\mathfrak{h}},\tilde G,\tilde h)$ defined by Eq. \eqref{eq:std form}, is referred to as the Standard form of $(\mathfrak{h}, G, h)$. We observe that when Eqs. \eqref{eq:bath def} and \eqref{eq:int temr} are written in standard form, $k$ is a scalar. We will use the standard form in theorems \ref{theorem dis 1} and \ref{theorem dis 2}.

\section{Discretisation of the continuous bath}\label{Discretisation of the continuous bath}
Let $P_n(x)$ be the real orthonormal polynomial of order $n$ with respect to the measure
$\mu_0(dx)=\pi^{-1}J(x)dx$ so that
\begin{equation}
    \int \mu_0(dx) P_n(x)P_m(x)=\delta_{m,n},\quad n,m\in {\mathbbm{N}}_0.
\end{equation}
The existence and uniqueness (up to a real phase) of these orthonormal polynomials is well
established \cite{gaut10}. We can use them to define a discretised Hamiltonian
\begin{equation}\label{discrete 1}
    \hat H_L=\hat H_\mathcal{S}\otimes{\mathbbm{1}}_\mathcal{B}+\hat A_\mathcal{S}\otimes\sum_{n=1}^L h_n^{(L)}(c_n^\dagger+c_n)+{\mathbbm{1}}_\mathcal{S}\otimes\sum_{n=1}^L\omega_n^{(L)}c_n^\dagger c_n,
    \quad L\in{\mathbbm{N}}_+
\end{equation}
where $\omega_n^{(L)}$ are the zeros of $P_L(x)$ and
\begin{equation}
    h_n^{(L)}=\frac{1}{\sqrt{\sum_{k=0}^{L-1} P^2_k(\omega_n^{(L)})}}
\end{equation}
The $c_n$ $(c_m^\dagger)$ on $\mathcal{B}$ are bosonic annihilation (creation) operators (and can
be expressed in terms of the field operators of Hamiltonian Eq. (\ref{the bathtub}) as detailed
by Eq. \eqref{eq:c n} in the proof). The discretised Hamiltonian, Eq. (\ref{discrete 1}) has a
nice interpretation in terms of Gauss quadrature rules: $\{\omega_n^{(L)} \}_{n=1}^L$ are the
Gauss knots and $\{(h^{(L)}_n)^2\}_{n=1}^L$ are the Gauss weights for the weight function $J(\cdot)/\pi$
(see section 1.4.2 for an introduction to Gauss quadrature and last paragraph of section 1.4.1
and Eq. (3.1.7) for the nodes and weights \cite{gaut10}). What is more, the theory of orthogonal
polynomials has established that $\omega_{min}\leq \omega_n^{(L)} \leq \omega_{max}$ for all
$L,n$ and $\{\omega_n^{(L)}\}_{n=1}^{L}$ interlace $\{\omega_n^{(L+1)}\}_{n=1}^{L+1}:$
\begin{equation}\label{interlace 1}
    \omega_{L+1}^{(L+1)}<\omega_{L}^{(L)}<\omega_{L}^{(L+1)}<\omega_{L-1}^{(L)}<\ldots<\omega_{1}^{(L)}<\omega_{1}^{(L+1)},
\end{equation}
thus discretisation in terms of the zeros of orthogonal polynomials corresponds to a natural
way of discretising a continuum. We note that the spectral density of Hamiltonian Eq. (\ref{discrete 1})
forms a pure point measure. This is in stark contrast with the spectral density of Hamiltonian
Eq. (\ref{the bathtub}), which is the weight function of an absolutely continuous measure with
respect to the Lebesgue measure.

This discretisation will lead to errors in time evolved observables on the system degrees
of freedom that can be bounded as determined by the following theorem. For self-adjoint
operators of the form $\hat O=\hat O_\mathcal{S}\otimes{\mathbbm{1}}_\mathcal{B}$, $\hat O_\mathcal{S}\in\mathcal{S}$
and initial normalised quantum states $\hat\varrho_0$ on $\mathcal{S}\otimes\mathcal{B}$ we find the following result.
\begin{theorem}\label{theorem dis 1}
    For any system observable $\hat O,$ $\| \hat O\|<\infty$ the error introduced on its expectation value at any time $t\geq 0$ when it's dynamics are approximated by the discretised Hamiltonian $\hat H_L,$ is bounded by
    \begin{align}\label{error H_L}
        \begin{split}
        &\bigl|\tr\bigl[\hat{O}\mathrm{e}^{-\mathrm{i} t\hat{H}}\hat{\varrho}_0\mathrm{e}^{\mathrm{i} t\hat{H}}\bigr]-\tr\bigl[\hat{O}\mathrm{e}^{-\mathrm{i} t\hat{H}_L}\hat{\varrho}_0\mathrm{e}^{\mathrm{i} t\hat{H}_L}\bigr]\bigr|^2 \\
        &\leq 8\eta_0 \|\hat O\|^2 \frac{\|\hat A_\mathcal{S}\|}{\omega_{max}}\frac{(\omega_{max} t)^{L+1} }{(L+1)!}  \left( \mathrm{e}^{\omega_{max} t}+1\right)\left( \|\vec
        \gamma_0\|^{1/2}+\eta_0\|\hat A_\mathcal{S}\|t \right),
     \end{split}
    \end{align}
for $L=1,2,3,\ldots$ where $\eta_0=\sqrt{\frac{2}{\pi}\int dx J(x)}$, and the operator $\vec \gamma_0$
encodes information regarding 2-point correlation functions in $\mathcal{B}$ of the initial state. More
precisely,
\begin{equation}\label{eq:gamma mat thorem 1}
\vec \gamma_0=
\left(\begin{array}{cc}
\vec \gamma_{xx} & \vec \gamma_{xp}\\
\vec \gamma_{px} & \vec \gamma_{pp}
\end{array}\right)
\end{equation}

where $\left[\vec\gamma_{xx}\right]_{n,m}=\tr[\hat x_n \hat x_m\hat\varrho_0]$, $\left[\vec\gamma_{xp}\right]_{n,m}=\tr[\hat x_n \hat p_m\hat\varrho_0]$, $\left[\vec\gamma_{px}\right]_{n,m}=\tr[\hat p_n \hat x_m\hat\varrho_0]$, $\left[\vec\gamma_{pp}\right]_{n,m}=\tr[\hat p_n \hat p_m\hat\varrho_0]$, $n,m\in{\mathbbm{N}}_+$ with
\begin{align}
    \begin{split}\label{eq:x p expansion}
       \hat x_n&=\frac{1}{\sqrt{2}}\int \mu_0^{(1/2)}(dx) P_{n-1}(x) \left(a_x^\dag+a_x\right),\\
       \hat p_n&=\frac{\mathrm{i}}{\sqrt{2}} \int \mu_0^{(1/2)}(dx) P_{n-1}(x) \left(a_x^\dag-a_x\right), \quad n\in {\mathbbm{N}}_+
    \end{split}
\end{align}
on $\mathcal{B}$ where $\mu_0^{(1/2)}(dx):=\sqrt{J(x)/\pi}\,dx$. Here $\hat x_n$ and $\hat p_n$ have been written in standard form (see section \ref{The Hamiltonian}).
\end{theorem}
See section \ref{LR bunds proof} for the proof of this theorem and appendix \ref{sec:alternative expressions for gamma 0} for alternative expressions for $\|\vec \gamma_0\|$ and examples including when $\tr_\mathcal{S}[\hat \varrho_0]$ is the vacuum state.

Theorem \ref{theorem dis 1} provides, to the best of our knowledge, for the first time, a bound on the
error that results from approximating a bath with absolutely continuous spectrum by a bath with pure point
spectrum in the form of a sum over discrete modes and vice versa.
The Lieb-Robinson light cone is achieved by choosing a tangent surface to the r.h.s. of Eq. \eqref{error H_L} bounding this non-linear function from above.
Prior to this work discretisations of continuum baths had a long history of being probed numerically
(see \cite{burk84,Kazansky1997,Shen08,iyer11,Ines15} and references here in). A fundamental insight was provided
by Burkey and Cantrell when they numerically observed that using the evenly spaced knots specified by the
trapezoid rule (which is referred to as \textit{Rice discretization}) seems not to be the most efficient
way to discretise a continuum \cite{burk84}. As pointed out in \cite{Shen08}, their choice of Gauss quadrature
rules to discretise the continuum is for the weight function $J(\cdot)/\pi,$ and thus theorem \ref{theorem dis 1}
applies to it. We will thus refer to the particular choice of Gauss weights and knots in Eq. \eqref{discrete 1}
as \textit{Burkey-Cantrell} discretisation.

However, Eq. (\ref{discrete 1}) is not the unique way to use Gauss quadrature rules to discretise
the continuum. The next theorem will present a sharper bound compared to that of Eq. \eqref{error H_L}
by discretising the bath according to the weight function $J(\sqrt{\cdot})/\pi$ instead. Via a trivial
change of variable followed by using the properties of $\mu_0(dx)$, one verifies that the measure
\begin{equation}\label{eq:phonon measure}
    \mu_1(dx):=\pi^{-1}J(\sqrt{x})dx,
\end{equation}
is determinate, with all moments finite. Again, we define $P_n'(x)$ as the real orthonormal
polynomial of order $n$ with respect to the measure $\mu_1(dx)$:
\begin{equation}\label{eq:ortho poly 2nd measure}
    \int \mu_1(dx) P_n'(x)P_m'(x)=\delta_{m,n},\quad n,m\in {\mathbbm{N}}_0.
\end{equation}
We use these orthogonal polynomials to define a discretised Hamiltonian:
\begin{equation}\label{discrete 2}
    \hat H_L'=\hat H_\mathcal{S}\otimes{\mathbbm{1}}_\mathcal{B}+\hat A_\mathcal{S}\otimes\sum_{n=1}^L \frac{h_n^{\prime(L)}}{\sqrt{2\bar{\omega}_n^{(L)}}}(c_n^{\prime\dagger}+c_n')+{\mathbbm{1}}_\mathcal{S}
    \otimes\sum_{n=1}^L\bar{\omega}_n^{(L)} c_n^{\prime\dagger} c_n', \quad L\in {\mathbbm{N}}_+
\end{equation}
where $\bar{\omega}_n^{(L)}=\sqrt{\omega_n^{\prime(L)}},$ $\omega_n^{\prime(L)}$ are the zeros of
$P_L'(x)$ and
\begin{equation}
    h_n^{\prime(L)}=\frac{1}{\sqrt{\sum_{k=0}^{L-1} P^{\prime2}_k(\omega_n^{\prime(L)})}}
\end{equation}
$c'_i$ $(c_j^{\prime\dagger})$ on $\mathcal{B}$ are bosonic annihilation (creation) operators (and
can be expressed in terms of the field operators of Hamiltonian Eq. (\ref{the bathtub}) as detailed
by Eq. \eqref{eq:d n} in the proof that is found in section 8.).

Here $\{\omega_n^{\prime(L)} \}_{n=1}^L$ are the Gauss knots and $\{(h^{\prime(L)}_n)^2\}_{n=1}^L$
are the Gauss weights for the weight function $J(\sqrt{\cdot})/\pi$. The knots satisfy
$\omega_{min}^2\leq \omega_n^{\prime(L)}\leq \omega_{max}^2$ for all $L,n$ and satisfy the same
interlacing properties as the Gauss knots $\{\omega_n^{(L)} \}_{n=1}^L$ observed in Eq. (\ref{interlace 1}).
To the best of our knowledge, this is a new discretisation, which, unlike the Burkey-Cantrell
discretisation, has not been probed numerically. We now state analogous theorems to that of
Eq. (\ref{error H_L}), but with a smaller r.h.s. for otherwise unchanged parameters. We will
make the distinction between when Hamiltonian Eq. (\ref{the bathtub}) is massless and massive.
\begin{theorem}\label{theorem dis 2}
    For any system observable $\hat O,$ $\| \hat O\|<\infty$ the error introduced on it's expectation value at any time $t\geq 0$ when it's dynamics are approximated by the discretised Hamiltonian $\hat H_L'$, is bounded by
    \begin{itemize}
    \item [] \textup{1)} Massive case \textup{($\omega_{min}>0$)}\\
    \begin{align}\label{error H_L'}
        \begin{split}
        &\bigl|\tr\bigl[\hat{O}\mathrm{e}^{-\mathrm{i} t\hat{H}}\hat{\varrho}_0\mathrm{e}^{\mathrm{i} t\hat{H}}\bigr]-\tr\bigl[\hat{O}\mathrm{e}^{-\mathrm{i} t\hat{H}_L'}\hat{\varrho}_0\mathrm{e}^{\mathrm{i} t\hat{H}_L'}\bigr]\bigr|^2 \\
        &\leq D_1\frac{(\omega_{max} t)^{2L+1} }{(2L+1)!}  \left( \mathrm{e}^{\omega_{max} t}+1\right)\left( \|\vec \gamma_0'\|^{1/2}+\eta_1\|\hat A_\mathcal{S}\| t \right),
        \end{split}
    \end{align}
    \item [] \textup{2)} Massless case \textup{($\omega_{min}=0$)}\\
    \begin{align}\label{error H_L'2}
        \begin{split}
        &\bigl|\tr\bigl[\hat{O}\mathrm{e}^{-\mathrm{i} t\hat{H}}\hat{\varrho}_0\mathrm{e}^{\mathrm{i} t\hat{H}}\bigr]
        -\tr\bigl[\hat{O}\mathrm{e}^{-\mathrm{i} t\hat{H}_L'}\hat{\varrho}_0\mathrm{e}^{\mathrm{i} t\hat{H}_L'}\bigr]\bigr|^2 \\
        &\leq D_1 \frac{(\omega_{max} t)^{2L+1}}{(2L+1)!}\left( \mathrm{e}^{ \omega_{max} t}+1\right) \left(
        \|\vec \gamma_0'\|^{1/2}+\eta_1\|\hat A_\mathcal{S}\|\frac{\mathrm{e}^{\omega_{max}t}-1}{\omega_{max}}
        \right) \mathrm{e}^{\omega_{max} t}
        \end{split}
    \end{align}
\end{itemize}
for $L=1,2,3,\ldots$ where $D_1=4\eta_1\|\hat O\|^2 \frac{\|\hat A_\mathcal{S}\|}{\omega_{max}}$,
$\eta_1=\sqrt{\frac{1}{\pi \omega_{max}}\int dx J(\sqrt{x})}$, and the operator $\vec \gamma_0'$ encodes information regarding 2-point correlation functions in $\mathcal{B}$ of the initial state. More precisely,
\begin{equation}\label{eq:gamma mat thorem 2}
    \vec \gamma_0'=
    \left(\begin{array}{cc}
    \vec \gamma_{xx}' & \vec \gamma_{xp}'\\
    \vec \gamma_{px}' & \vec \gamma_{pp}'
    \end{array}\right),
\end{equation}
where $\left[\vec\gamma_{xx}'\right]_{n,m}=\tr[\hat x_n' \hat x_m'\hat\varrho_0]$, $\left[\vec\gamma_{xp}'\right]_{n,m}=\tr[\hat x_n' \hat p_m'\hat\varrho_0]$,
$\left[\vec\gamma_{px}'\right]_{n,m}=\tr[\hat p_n' \hat x_m'\hat\varrho_0]$,
$\left[\vec\gamma_{pp}'\right]_{n,m}=\tr[\hat p_n' \hat p_m'\hat\varrho_0]$,
$n,m\in{\mathbbm{N}}_+$ with
\begin{align}
    \begin{split}\label{eq:x p expansion 2}
    \hat x_n'&=\sqrt{\omega_{max}}\int \mu_0^{(1/2)}(dx) P'_{n-1}(x^2)
    \left(a_x^\dag+a_x\right),\\
    \hat p_n'&=\frac{\mathrm{i}}{\sqrt{\omega_{max}}} \int \mu_0^{(1/2)}(dx) P'_{n-1}(x^2)\,x
    \left(a_x^\dag-a_x\right), \quad n\in {\mathbbm{N}}_+
    \end{split}
\end{align}
on $\mathcal{B}$. Here $\hat x_n'$ and $\hat p_n'$ have been written in standard form (see section \ref{The Hamiltonian}).
\end{theorem}
See section \ref{LR bunds proof} for proof. See section \ref{sec:Alternative expressions for vec gamma_0'}
for alternative expressions for $\|\vec \gamma_0'\|$\\

{\em Observations} The bound that has been obtained in theorem \ref{theorem dis 2} up to constant factors achieves
the same error as the Burkey-Cantrell discretisation but with only half the number of knot points.
This observation suggests an improved discretisation method using Gauss quadrature rules and may also
find applications in analytical and numerical work.

Interestingly, for theorems (\ref{theorem dis 1}) and (\ref{theorem dis 2}) we used a Lieb-Robinson bound to derive a relation between two non-local Hamitonians (in Eq. (\ref{the bathtub}) and its discretised counterparts Eqs. (\ref{discrete 1}) and (\ref{discrete 2}), every harmonic oscillator is coupled \textit{directly} to the system via $\hat A_\mathcal{S}$). But now, in these non-local Hamiltonians, the notion of a distance that is  normally associated with a lattice in the Lieb-Robinson bound plays the role of the number of Gauss knot points $L$.

Note also that for $\omega_{max}=\infty$, the upper bound in Eqs. (\ref{error H_L}) and (\ref{error H_L'})
diverge. Under the minimal assumptions that have been made in this manuscript regarding the initial
state $\hat \varrho_0$, this is to be expected since for $\omega_{max}=\infty$ we are sampling an
unbounded interval with a finite number of sample points\footnote{This observation, gives an intuitive
explanation to why certain types of Hamiltonians can be proven \textit{not} to admit a Lieb-Robinson
bound \cite{eise09}, since when the lattice is diagonalised via a unitary transformation, one should
expect to observe similar divergences}. The theory of orthogonal polynomials also tells us that if the
support of the spectral density is gapped, i.e. the spectral density vanishes strictly in an interval
$[\omega_i,\omega_f]$ with $\omega_{min}<\omega_i<\omega_f<\omega_{max}$, a situation known to occur in physical systems such as photonic crystals \cite{prior13},
then the discretised bath will have at most one discrete mode in the gap $[\omega_i,\omega_f]$.


\section{Multiple baths extension}\label{Multiple baths extension}
Often in physical settings, the quantum system on $\mathcal{S}$ is coupled to multiple continuous
baths via different interaction terms. In such circumstances, the Hamiltonian is on $\mathcal{S}\otimes
\mathcal{B}^{\otimes N}$ and reads
\begin{equation}
    \label{the bathtub multi}
    \hat H_\text{mul}=\hat H_\mathcal{S}\otimes{\mathbbm{1}}_\mathcal{B}^{\otimes N} + \sum_{m=1}^N \hat H^{(m)},
    \quad N\in {\mathbbm{N}}_+
\end{equation}
where
\begin{align}
    \begin{split}
        \hat H^{(m)}=&\hat A_\mathcal{S}^{(m)}\otimes{\mathbbm{1}}_\mathcal{B}^{\otimes m-1}\otimes\left( \int dk h^{(m)}(k) (a_k^{(m)\dag}+a_k^{(m)}) \right)\otimes{\mathbbm{1}}_\mathcal{B}^{\otimes N-m}\\
        &+{\mathbbm{1}}_\mathcal{S}\otimes{\mathbbm{1}}_\mathcal{B}^{\otimes m-1}\otimes\left(\int dk g^{(m)}(k)a^{(m)\dag}_k a_k^{(m)}\right)\otimes{\mathbbm{1}}_\mathcal{B}^{\otimes N-m}.
    \end{split}
\end{align}
Here each individual bath and interaction term $\hat H^{(m)}$ is defined as in Eq. (\ref{the bathtub})
and will have its own spectral density $J^{(m)}.$ Since they are independent bosonic baths, we also
have $[a_x^{(n)},a_y^{(m)\dag}]=\delta_{n,m}\delta(x-y),$ $[a_x^{(n)},a_y^{(m)}]=0.$ In analogy with
Eqs. (\ref{discrete 1}) and (\ref{discrete 2}), we can define discretised versions of $\hat H^{(m)}$
according to the two discretisation schemes considered in this article. Similarly to Eq.
(\ref{discrete 1}), for Burkey-Cantrell discretisation, we define
\begin{align}\label{discrete 1 n}
    \begin{split}
        \hat H^{(0,m)}_L=&\hat A_\mathcal{S}^{(m)}\otimes{\mathbbm{1}}_\mathcal{B}^{\otimes m-1}\otimes\left(\sum_{n=1}^L h_n^{(m,L)}(c_n^{(m)\dag}+c_n^{(m)})\right)\otimes{\mathbbm{1}}_\mathcal{B}^{\otimes N-m}\\
        &+{\mathbbm{1}}_\mathcal{S}\otimes{\mathbbm{1}}_\mathcal{B}^{\otimes m-1}\otimes\left(\sum_{n=1}^L\omega_n^{(m,L)}c_n^{(m)\dag} c_n^{(m)}\right)\otimes{\mathbbm{1}}_\mathcal{B}^{\otimes N-m}, \quad L,m\in {\mathbbm{N}}_+
    \end{split}
\end{align}
where the Gauss knots $\omega_n^{(m,L)}$ and Gauss weights $(h_n^{(m,L)})^2$ are calculated
from the weight function $J^{(m)}(\cdot)/\pi,$ for which we denote its minimum and maximum frequencies
by $\omega_{min}^{(m)}$ and $\omega_{max}^{(m)}$ respectively (see paragraph below Eq. (\ref{SD_eq})).
Similarly to the second discretisation procedure we considered, Eq. (\ref{discrete 2}), we define
\begin{align}\label{discrete 2 n}
    \begin{split}
        \hat H_L^{(1,m)}=&\hat A_\mathcal{S}^{(m)}\otimes{\mathbbm{1}}_\mathcal{B}^{\otimes m-1}\otimes\left(\sum_{n=1}^L \frac{h_n^{\prime(m,L)}}{\sqrt{2\bar{\omega}_n^{(m,L)}}}(c_n^{\prime(m)\dagger}+c_n^{\prime(m)})\right)
        \otimes{\mathbbm{1}}_\mathcal{B}^{\otimes N-m}\\
        &+{\mathbbm{1}}_\mathcal{S}\otimes{\mathbbm{1}}_\mathcal{B}^{\otimes m-1}\otimes\left(\sum_{n=1}^L\bar{\omega}_n^{(m,L)} c_n^{\prime(m)\dagger} c_n^{\prime(m)}\right) \otimes{\mathbbm{1}}_\mathcal{B}^{\otimes N-m}, \quad L,m\in {\mathbbm{N}}_+.
    \end{split}
\end{align}
We can now define the discretised Hamiltonian
\begin{equation}
    \hat H^{\vec q}_{\vec L}=\hat H_\mathcal{S}\otimes{\mathbbm{1}}_\mathcal{B}^{\otimes N}
    + \sum_{m=1}^N \hat H_{L_m}^{(q_m,m)}, \quad N\in {\mathbbm{N}}_+
\end{equation}
where $\vec q=(q_1,q_2,\ldots,q_N)$ is a binary string and $\vec L=(L_1,L_2,\ldots,L_N).$
With the definitions
\begin{align}
    \Big(\Delta^{(0,m)}_{L}(t)\Big)^2 = D_0^{(m)}\frac{(\omega_{max}^{(m)} t)^{L+1} }{(L+1)!}
    \left( \mathrm{e}^{\omega_{max}^{(m)} t}+1\right)\left( \|\vec \gamma_0^{(m)}\|^{1/2} +
    \eta_0^{(m)}\|\hat A_\mathcal{S}\|t \right),
\end{align}
\begin{align}
    \begin{split}
        \Big(\Delta^{(1,m)}_{L}(t)\Big)^2&=
        \begin{cases}
            f_L^{(m)}(t)  &\mbox{if } \omega_{min}^{(m)}>0\\
        	g_L^{(m)}(t)  & \mbox{if } \omega_{min}^{(m)}=0,
        \end{cases}\\
        f_L^{(m)}(t)&=C_1^{(m)} \frac{(\omega_{max}^{(m)} t)^{2L+1} }{(2L+1)!} ( \mathrm{e}^{ \omega_{max}^{(m)} t}+1)
        \left( \|\vec\gamma_0^{(m)}\|^{1/2}+\eta_1^{(m)}\|\hat A_\mathcal{S}^{(m)}\|t \right)\\
        g_L^{(m)}(t)&=C_1^{(m)} \frac{(\omega_{max}^{(m)} t)^{2L+1}}{(2L+1)!}( \mathrm{e}^{ 2\omega_{max}^{(m)} t}+1)\left( \|\vec \gamma_0^{(m)}\|^{1/2}+\eta_1^{(m)}(\mathrm{e}^{\omega_{max}^{(m)}t}-1)\|\hat A_\mathcal{S}^{(m)}\| \right)\mathrm{e}^{\omega_{max}^{(m)} t}
    	\end{split}
	\end{align}	
	
\begin{align}
    \begin{split}
        D_0^{(m)}&=8\eta_0^{(m)} \|\hat O\|^2 \|\hat A_\mathcal{S}^{(m)}\|/\omega_{max}^{(m)},\quad
        \eta_0^{(m)}=\sqrt{\frac{2}{\pi}\int dx J^{(m)}(x)},\\
        C_1^{(m)}&=4 \eta_1^{(m)} \|\hat O\|^2 \|\hat A_\mathcal{S}^{(m)}\|/\omega_{max}^{(m)},\quad \eta_1^{(m)}=\sqrt{\frac{1}{\pi\omega_{max}^{(m)}}\int dx J^{(m)}(\sqrt{x})},
    \end{split}
\end{align}
where $\Delta^{(q,m)}_{L}\geq 0,$ we have the following bounds.
\begin{corollary}\label{multi cor}
    For any system observable $\hat O,$ $\| \hat O\|<\infty$ the error introduced on it's expectation value at any time $t\geq 0$ when it's dynamics are approximated by the discretised Hamiltonian $\hat H^{\vec q}_{\vec L}$, is bounded by
    \begin{equation}\label{multi eq}
    \left| \tr\bigl[\hat{O}\mathrm{e}^{-\mathrm{i} t\hat{H}_\textup{mul}}\hat{\varrho}_0\mathrm{e}^{\mathrm{i} t\hat{H}_\textup{mul}}\bigr]-\tr\bigl[\hat{O}\mathrm{e}^{-\mathrm{i} t \hat H^{\vec q}_{\vec L}}\hat{\varrho}_0\mathrm{e}^{\mathrm{i} t\hat H^{\vec q}_{\vec L}}\bigr] \right|\leq \sum_{m=1}^N \Delta^{(q_m,m)}_{L_m}(t).
    \end{equation}
\end{corollary}
See section \ref{multiple chains extension proof} for proof. We thus see that the error
incurred scales linearly in the number of discretised continua.

\section{Explicit examples}\label{Explicit examples}
For a range of spectral densities from the literature we will now present the explicit expressions
for the frequencies and Gauss weights.

Consider the semi-circle law spectral density:
\begin{equation}
    J(\omega)=C\sqrt{(\omega_{max}-\omega)(\omega-\omega_{min})},
\end{equation}
for some constant $C>0$. This is the weight function of the Chebyshev polynomials of the second
kind on the interval $[\omega_{min},\omega_{max}]$. Their corresponding orthogonal polynomials
and their zeros are known explicitly \cite{gaut10}. The zeros of the $n$th order Chebyshev
polynomial on $[-1,1]$ are
\begin{equation}
    x_k=\cos \left( \frac{k}{n+1}\pi \right),\quad k=1,\ldots,n
\end{equation}
and hence for the Burkey-Cantrell discretisation (Eq. (\ref{discrete 1})), the discrete frequencies are
\begin{equation}
    \omega_k^{(L)}=\frac{(\omega_{min}-\omega_{max})}{2}\cos \left( \frac{k}{L+1}\pi \right) +\frac{(\omega_{max}+\omega_{min})}{2}.
\end{equation}
For the spectral density of the \textit{Rubin model} \cite{Weiss2001} we find
\begin{equation}
    J(\omega)=C\sqrt{(\omega^2_{max}-\omega^2)(\omega^2-\omega^2_{min})},
\end{equation}
for some constant $C>0$. From Eq. (\ref{ph map dens}), we have that this corresponds to
a Chebyshev weight function of the second kind for the measure $J(\sqrt{\cdot})/\pi$ used
in Eq. \eqref{eq:ortho poly 2nd measure}. Thus we have for the 2nd discretisation method
(Eq. (\ref{discrete 2}))
\begin{equation}
    \bar{\omega}_k^{(L)}=\sqrt{\frac{(\omega^2_{min}-\omega^2_{max})}{2}\cos \left( \frac{k}{L+1}\pi \right) +\frac{(\omega^2_{max}+\omega^2_{min})}{2}.}
\end{equation}
Note that since the Chebyshev polynomials are known explicitly, we can also find explicitly the
$h_n^{(L)}$ coefficients in both the above examples.

Now we will consider the frequently considered power-law spectral densities
\begin{equation}
    J(\omega)=2\pi \alpha (\omega_{max}-\omega_{min})(\omega-\omega_{min})^s,\quad -1<s
\end{equation}
which include the sub-ohmic $s<1$, ohmic $s=1$ and the super-ohmic $s>1$ case. We can use
the Jacobi Polynomials to describe the measures $J(\cdot)/\pi$ and $J(\sqrt{\cdot})/\pi$.
For the Jacobi polynomials $P_n^{(\alpha,\beta)}(x)$ on $[-1,1]$, we have the Buell inequalities
for their zeros $x_i$ (\cite{gobo75}, pg 125):
\begin{equation}
    \frac{i+(\alpha+\beta-1)/2}{n+(\alpha+\beta+1)/2}\pi<\nu_i<\frac{i}{n+(\alpha+\beta+1)/2}\pi,\quad i=1,\ldots,n
\end{equation}
where $\nu_i=\arccos x_i$ and $-1/2\leq \alpha\leq 1/2$, $-1/2\leq \beta\leq 1/2$, excluding the case $\alpha^2=\beta^2=1/4$.
Hence for the Burkey-Cantrell discretisation (Eq. (\ref{discrete 1}))
\begin{align}\label{pa f J}
    \begin{split}
        &\left[ 1-\cos\left(\frac{k}{L+(s+1)/2}\pi\right)\right]\frac{\omega_{max}-\omega_{min}}{2}
        < \omega_k^{(L)}-\omega_{min} \\
        &<\left[1-\cos\left(\frac{k+(s-1)/2}{L+(s+1)/2}\pi\right)\right]
        \frac{\omega_{max}-\omega_{min}}{2}
    \end{split}
\end{align}
for $-1/2\leq s\leq 1/2$. For the massless ($\omega_{min}=0$) 2nd discretisation method (Eq. (\ref{discrete 2})),
\begin{align}\label{ph f J}
    \begin{split}
        &\omega_{max}\sqrt{\left[ 1-\cos\left(\frac{k}{L+(s/2+1)/2}\pi\right)\right]/2}< \bar{\omega}_k^{(L)}\\
        &<\omega_{max}\sqrt{\left[1-\cos\left(\frac{k+(s/2-1)/2}{L+(s/2+1)/2}\pi\right)\right]/2}
    \end{split}
\end{align}
for $-1< s\leq 1$. Note that other bounds for the zeros are known for the values of
$-1<s$ not covered here. See for example \cite{elbe94}. In all the examples in this
section, the corresponding orthogonal polynomials are known explicitly and thus one can achieve explicit expressions for the Gauss weights too.

\section{Conclusion}
We derive two error bounds for dynamical observables when discretising a
continuum of harmonic oscillators according to Gauss quadrature rules. For one of the bounds,
numerical studies have probed this discretisation numerically in the past. For the second case,
no prior numerical studies have been performed. The second bound achieves a sharper bound for the
same parameters, suggesting that the second Gauss quadrature discretisation method may be
more efficient.

To prove these results, we make use of a unitary mapping based on orthogonal polynomials 
and a Lieb-Robinson bound, providing yet another application to these powerful tools. Attempts 
to understand the errors endured by dynamical observables as a function of time when approximating 
a Hamiltonian with absolutely continuous spectrum with one of pure point spectrum, although 
receiving interest since 1929 by various authors, has up to now only consisted in numerical studies.

\textit{Acknowledgements -- } M.W. would like to thank Marcus Cramer for motivating discussions. He
also acknowledges support from MOE Tier 3A Grant MOE2012-T3-1. MBP ackowledges support
by an Alexander von Humboldt Professorship, the ERC grant BioQ and the EU projects PAPETS
and SIQS.

\section{Proofs}\label{Proofs}
\subsection{Proofs of theorems \ref{theorem dis 1} and \ref{theorem dis 2}}\label{LR bunds proof}
For convenience, we will prove both theorems in parallel. The proof proceeds along three main steps.
First, we will state unitary transformations of the bath modes which allow Eq. (\ref{the bathtub})
to be written as a Hamiltonian in which the system on $\mathcal{S}$ couples to the first
site of a nearest neighbour coupled harmonic chain. This chain is then truncated at finite length
and a Lieb-Robinson bond is deployed to achieve error bounds for the expectation
value of systems operators. Subsequently, the finite chain is transformed back to a non-local Hamiltonian
and it is demonstrated that this achieves the discretised Hamiltonians as formulated
in theorems \ref{theorem dis 1} and \ref{theorem dis 2}.

In \cite{misc14}, it was shown that if $J(\cdot)$ has all moments finite, the Hamiltonian Eq. (\ref{the bathtub})
can be written in the forms
\begin{align}\label{marticle mapping hamiltonian eq}
\begin{split}
\hat H=& \hat H_\mathcal{S}\otimes{\mathbbm{1}}_\mathcal{B}+\sqrt{\beta_0(0)} \hat A_\mathcal{S}\otimes(b_0(0)+b_0^\dagger(0))\\
&+{\mathbbm{1}}_\mathcal{S}\otimes\sum_{n=0}^\infty \alpha_n(0)b_n^\dagger(0) b_n(0) +\sqrt{\beta_{n+1}(0)}(b_{n+1}^\dagger(0) b_n(0) +h.c.)
\end{split}
\end{align}
and
\begin{align}\label{phonon chain hamiltonian}
\begin{split}
\hat H=& \hat H_{S}\otimes{\mathbbm{1}}_\mathcal{B}+\sqrt{\frac{\beta_0(1)}{\omega_{max}}}\hat A_\mathcal{S}\otimes\hat X_0\\
&+{\mathbbm{1}}_\mathcal{S}\otimes\sum_{n=0}^{\infty} \left(\frac{\sqrt{\beta_{n+1}(1)}}{\omega_{max}}\hat X_n\hat X_{n+1}+ \frac{\alpha_n(1)}{2\omega_{max}}\hat X^2_n +\frac{\omega_{max}}{2}\hat P_n^2\right)\;\\
&+K{\mathbbm{1}}_\mathcal{S}\otimes{\mathbbm{1}}_\mathcal{B}.
\end{split}
\end{align}
where
\begin{equation}\label{eq:phnon ops}
\hat X_n:=\sqrt{\omega_{max}}\,(b_n^\dag(1)+b_n(1)),\quad \hat P_n:=\mathrm{i} (b_n^\dag(1)-b_n(1))/(2\sqrt{\omega_{max}}),\quad n\in{\mathbbm{N}}_0,
\end{equation}
are position and momentum operators satisfying the canonical commutation relations. It was established in Theorem 33 of \cite{misc14} that the operators $b_n(q),b^\dag_n(q)$ are well defined on $\mathcal{D}(\mathfrak{n}^{1/2})$, where $\mathfrak{n}=d\hat{\mathrm{\Gamma}}({\mathbbm{1}}_\mathfrak{h})$ is the second quantised boson number operator. $K\in {\mathbbm{R}}$ thus w.l.o.g.
from here on it will be set to zero since it does not contribute to Eq. (\ref{error H_L'}).
Now we will state a recently derived locality bound \cite{ourLR}.

The Hamiltonian the bound applies to is system with Hamiltonian $\hat{H}_\mathcal{S}$ on $\mathcal{S}$ coupled to a semi-infinite nearest neighbour bosonic chain of the form
\begin{equation}\label{the bath}
\hat{H}_\mathcal{B}=\frac{1}{2}\sum_{i,j=0}^\infty\left(\hat{x}_iX_{i,j}\hat{x}_j+\hat{p}_iP_{i,j}\hat{p}_j\right),
\end{equation}
where $X_{i,j}=X_{j,i}\in{\mathbbm{R}}$,  $P_{i,j}=P_{j,i}\in{\mathbbm{R}}$ where $X_{i,j}=P_{i,j}=0$ for $|i-j|>1$.
The system-bath coupling is of the form
$\hat{h}\otimes\hat{x}_0,$ where $\hat{h}$ acts on the system and is bounded in operator norm.
The total Hamiltonian reads
\begin{equation}\label{tot ham}
\hat{H}=\hat{H}_\mathcal{S}\otimes{\mathbbm{1}}_\mathcal{B}+\hat{h}\otimes\hat{x}_0+{\mathbbm{1}}_\mathcal{S}\otimes\hat{H}_\mathcal{B}.
\end{equation}
We then define a spatially truncated Hamiltonian
\begin{equation}
\begin{split}
\label{eq:H_L}
\hat{H}_L&=\hat{H}_\mathcal{S}\otimes{\mathbbm{1}}_\mathcal{B}+\hat{h}\otimes\hat{x}_0+{\mathbbm{1}}_\mathcal{S}\otimes\hat{H}^L_\mathcal{B},\\
\hat{H}^L_\mathcal{B}&=\frac{1}{2}\sum_{i,j=0}^{L-1}\left(\hat{x}_iX_{i,j}\hat{x}_j+\hat{p}_iP_{i,j}\hat{p}_j\right)
\end{split}
\end{equation}
and the constant $c$ such that $\|X_LP_L\|^{1/2}\le c$. Theorem 1 in \cite{ourLR} for $X,P>0$ or $X=P$ gives us the bound
\begin{equation}\label{bound 1}
\frac{\bigl|\tr\bigl[\hat{O}\mathrm{e}^{-\mathrm{i}\hat{H}t}\hat{\varrho}_0\mathrm{e}^{\mathrm{i}\hat{H}t}\bigr]-\tr\bigl[\hat{O}\mathrm{e}^{-\mathrm{i}\hat{H}_Lt}\hat{\varrho}_0\mathrm{e}^{\mathrm{i}\hat{H}_Lt}\bigr]\bigr|^2}{4\|\hat{O}\|^2\|\hat{h}\|/c}
\le
\frac{C(\|\vec \gamma_0\|^{\frac{1}{2}}+\|\hat{h}\|t)(c t)^{L+1}(\mathrm{e}^{ct}+1)}{(L+1)!}
,
\end{equation}
where
\begin{equation}
\begin{split}
C&= \|P_L\||X_{L-1,L}|/c^2+|P_{L-1,L}|/c
\end{split}
\end{equation}
and
\begin{equation}\label{eq:gamma 0 general}
\vec \gamma_0=\left(\begin{array}{cc}
\vec\gamma_{xx} & \vec\gamma_{xp}\\
\vec\gamma_{px} & \vec\gamma_{pp}
\end{array}\right),\;\;\; [\vec \gamma_{ab}]_{i,j}=\tr[\hat{a}_i\hat{b}_j\hat{\varrho}_0],\;\;\; a,b=x,p,
\end{equation}
collects the two-point bath correlations in the initial state. Also note that $|X_{L-1,L}|\le \|X_L\|\leq \|X\|$ and $|P_{L-1,L}|\le \|P_L\|\leq \|P\|$.
If $X,P>0$ or $X=P$ are not satisfied, we use theorem 3 in \cite{ourLR} to achieve the bound. Let $c'$ such that $\{\|X\|,\|P\|\}\leq c',$ then for all $X,$ $P$
\begin{equation}
\begin{split}\label{bound2}
&\frac{\bigl|\tr\bigl[\hat{O}\mathrm{e}^{-\mathrm{i}\hat{H}t}\hat{\varrho}_0\mathrm{e}^{\mathrm{i}\hat{H}t}\bigr]-\tr\bigl[\hat{O}\mathrm{e}^{-\mathrm{i}\hat{H}_Lt}\hat{\varrho}_0\mathrm{e}^{\mathrm{i}\hat{H}_Lt}\bigr]\bigr|^2}{4\|\hat{O}\|^2\|\hat{h}\|/c}\\
&\hspace{4cm}\le
\frac{C\mathrm{e}^{c^\prime t}(\|\vec \gamma_0\|^{1/2}+\|\hat{h}\|
\frac{\mathrm{e}^{c^\prime t}-1}{c^\prime})(ct)^{L+1}(\mathrm{e}^{ct}+1)}{(L+1)!}.
\end{split}
\end{equation}
If $P\propto{\mathbbm{1}}$, we may replace the factor $(ct)^{L+1}/(L+1)!$ by $(ct)^{2L+1}/(2L+1)!$ in the R.H.S. of Eqs. \eqref{bound 1}, \eqref{bound2}. We can readily apply these bounds to Eqs. (\ref{marticle mapping hamiltonian eq}) and (\ref{phonon chain hamiltonian}). First define position and momentum operators
\begin{equation}\label{eq:particle pos momentum ops def}
\hat x_n(0):=(b^\dagger_n(0)+b_n(0))/\sqrt{2},\quad \hat p_n(0):=\mathrm{i}(b^\dagger_n(0)-b_n(0))/\sqrt{2},\quad n\in{\mathbbm{N}}_0.
\end{equation}
Comparing Eq. (\ref{the bath}) with Eqs. (\ref{marticle mapping hamiltonian eq}), (\ref{phonon chain hamiltonian}) and the definition of the Jacobi matrices $\mathcal{J}(d\lambda^q)$ (see Eq. (162) in \cite{misc14}), we find\footnote{Here we use the notation of \cite{misc14} for reference purposes. In this article, we have used the notation $\mu_0(dx)\equiv d\lambda^0(x)$ and $\mu_1(dx)\equiv d\lambda^1(x)$. We will continue to use the notation of \cite{misc14} throughout this proof for these two measures}\\
\begin{itemize}
\item[] For Eq. (\ref{marticle mapping hamiltonian eq}):
\begin{equation}
X=P=\mathcal{J}(d\lambda^0),\quad \hat h=\sqrt{2\beta_0(0)}\hat A_\mathcal{S},\quad d\lambda^0(x)=\pi^{-1}J(x)dx.
\end{equation}
\item[] For Eq. (\ref{phonon chain hamiltonian}):
\begin{equation}
X=\frac{\mathcal{J}(d\lambda^1)}{\omega_{max}},\quad P={\mathbbm{1}}\;\omega_{max}, \quad\hat h=\sqrt{\frac{\beta_0(1)}{\omega_{max}}}\hat A_\mathcal{S},\quad d\lambda^1(x)=\pi^{-1}J(\sqrt{x})dx.\label{ph map dens}
\end{equation}
\end{itemize}
From Eqs (15,156,160) in \cite{misc14},
\begin{equation}
\beta_0(0)=\int dxJ(x)/\pi,\quad \beta_0(1)=\int dx J(\sqrt{x})/\pi.
\end{equation}
Since the spectrum of a Jacobi matrix is equal to its minimally closed support interval \cite{Teschl2000}, we have for the Eqs. (\ref{marticle mapping hamiltonian eq}), (\ref{phonon chain hamiltonian}): $\|X\|=\|P\|=\sqrt{\|XP\|}=\omega_{max}$, and $X>0$ iff $\omega_{min}>0$.
The r.h.s. of Eqs. (\ref{error H_L}), (\ref{error H_L'}) and (\ref{error H_L'2}) are a direct consequence of these bounds. We will now proceed to apply another unitary transformation. This time, we will apply it to the above spatially truncated Hamiltonians to write them in terms of Gauss quadrature.\\
Both Eqs. (\ref{marticle mapping hamiltonian eq}) and (\ref{phonon chain hamiltonian}) can be written in the compact form (see Eq. (162) in \cite{misc14})
\begin{align}\label{the generalised mapping in terms of the jacobi mat eq}
\begin{split}
H=&\hat H_{S}\otimes{\mathbbm{1}}_\mathcal{B}+\sqrt{\beta_0(q)}\,\hat A_\mathcal{S}\otimes(b^\dagger(q)+b(q))\\
&+{\mathbbm{1}}_\mathcal{S}\otimes\frac{q}{2}\left[ \vec b^\top (q)\left(\mathcal{J}(d\lambda^q)-\frac{q}{4}{\mathbbm{1}}\right)\vec b(q) +h.c.\right]\\
&+{\mathbbm{1}}_\mathcal{S}\otimes\vec b^\dagger(q) \left(\mathcal{J}(d\lambda^q)+\frac{q}{4}{\mathbbm{1}}\right)\vec b(q)
\end{split}
\end{align}
where
\begin{align}
\begin{split}
\vec b^\dagger(q)&:=(b_0^\dagger(q),b_{1}^\dagger(q),b_{2}^\dagger(q),b_{3}^\dagger(q),\ldots),\\
\vec b(q)&:=(b_0(q),b_{1}(q),b_{2}(q),b_{3}(q),\ldots)^\top ,
\end{split}
\end{align}
and Eq. (\ref{marticle mapping hamiltonian eq}) is given when $q=0$ while Eq. (\ref{phonon chain hamiltonian})
when $q=1$. In the case of the truncated chains $\hat H_L$and $\hat H_L'$, we simply replace
$\mathcal{J}$ with $\mathcal{J}_{L}:=\mathcal{J}_{[1:L;1:L]},$ ${\mathbbm{1}}$ with ${\mathbbm{1}}_{[1:L;1:L]}$
and $\vec b$ with $\vec b_{[1:L]}$.
Since Jacobi matrices are real symmetric, they are diagonalisable via an orthogonal transformation. Thus
\begin{equation}
    \mathcal{J}_L(d\lambda^q)=O_{L}(d\lambda^q)D^{(L)}(d\lambda^q)O^{\top}_{L}(d\lambda^q)
\end{equation}
where $\textup{diag}D^{(L)}(d\lambda^q)=(\omega^{(L)}_1 (d\lambda^q),\ldots,\omega^{(L)}_L(d\lambda^q))$
are the eigenvalues of $\mathcal J (d\lambda^q)$ and $O_{L}(d\lambda^q)=(\vec v_1(d\lambda^q),\ldots,\vec v_L(d\lambda^q))$ with $\vec v_i(d\lambda^q)$ the normalised eigenvector for $\omega^{(L)}_i (d\lambda^q)$.
Let $\tilde P_n(d\lambda^q,x)$ be the real orthogonal polynomial of order $n$ corresponding to the
measure $d\lambda^q(x)$:
\begin{equation}
    \int d\lambda^q(x) \tilde P_n(d\lambda^q;x)\tilde P_m(d\lambda^q;x)=||\tilde P_n ||_{d\lambda^q}^2\,\delta_{m,n},\quad n,m\in {\mathbbm{N}}_0
\end{equation}
where $m=n$ defines $|| \tilde P_n ||_{d\lambda^q}$.
The set $\{\tilde P_n(x;d\lambda^q)\}_{n=0}^\infty$ exists and for a specific choice of
normalisation $\{|| \tilde P_n ||_{d\lambda^q}\}_{n=0}^\infty$, it is unique up to a real
phase \cite{gaut10}. It is known that  $\{\omega_i^{(j)}(d\lambda^q)\in {\mathbbm{R}}\}_{i=1}^j$ are
the zeros of $\tilde P_j(d\lambda^q;x)$ and
\begin{align}
\begin{split}
    &\vec v_i(d\lambda^q)=\\ &\vec
    (\tilde P_0(d\lambda^q,\omega_i^{(L)}(d\lambda^q)),\ldots,\tilde  P_{L-1}(d\lambda^q,\omega_i^{(L)}(d\lambda^q)))^\top \frac{1}{\sqrt{\sum_{k=0}^{L-1}\tilde{P}^2_k(d\lambda^q;\omega_i^{(L)}(d\lambda^q))}}
\end{split}
\end{align}  \cite{gaut10}.
Let us define
\begin{equation}
    h_j^{(L)}(d\lambda^q):=\frac{\sqrt{\beta_0(d\lambda^q)}\tilde P_0(d\lambda^q;\omega_j^{(L)}(d\lambda^q))}
    {\sqrt{\sum_{k=0}^{L-1}\tilde{P}^2_k(d\lambda^q;\omega_j^{(L)}(d\lambda^q))}},\quad j\in {\mathbbm{N}}_+.
\end{equation}
Hence for Eq. (\ref{marticle mapping hamiltonian eq}):
\begin{equation}
\hat H_L=\hat H_\mathcal{S}\otimes {\mathbbm{1}}_\mathcal{B}+\hat A_\mathcal{S}\otimes\sum_{n=1}^L h_n^{(L)}(d\lambda^0)(c_n^\dagger+c_n)+{\mathbbm{1}}_\mathcal{S}\otimes\sum_{n=1}^L\omega_n^{(L)}(d\lambda^0)c_n^\dagger c_n
\end{equation}
where $\vec c:=O^{\top}_L(d\lambda^0)\vec b(0)_{[1:L]},\;$ $\vec c^\dag:=\vec b^\dag(0)_{[1:L]} O_L(d\lambda^0)$.
We can easily verify that $[c_i,c_j^\dagger]=\delta_{i,j},$ $[c_i,c_j]=0$ and are hence bosonic creation annihilation operators. If we use orthonormal polynomials (i.e. $||\tilde P_n||_{d\lambda^q}=1$), then the zeroth order polynomial
is $\tilde P_0(d\lambda^q;x)=1/\sqrt{\beta_0(q)}$ and hence we find the results stated in theorem \ref{theorem dis 1}.
For  Eq. (\ref{phonon chain hamiltonian}), we have
\begin{equation}
\hat H_L=\hat H_\mathcal{S}\otimes {\mathbbm{1}}_\mathcal{B}+\hat A_\mathcal{S}\otimes\sum_{n=1}^L\frac{h_n^{(L)}(d\lambda^1)}{\sqrt{2\bar{\omega}_n^{(L)}}}(c_n^{\prime\dagger}+c'_n)+{\mathbbm{1}}_\mathcal{S}\otimes\sum_{n=1}^L\bar{\omega}_n^{(L)}c_n^{\prime\dagger} c'_n
\end{equation}
where $\bar{\omega}^{(L)}_n:=\sqrt{\omega^{(L)}_n(d\lambda^1)},$
\begin{align}
\begin{split}
 c'_n:=&\sqrt{\frac{\bar{\omega}_n^{(L)}}{2}}\Bigl[\left(1+\frac{1}{2\bar{\omega}_n^{(L)}}\right)\left[O^{\top}_L(d\lambda^1)\vec b(1)_{[1:L]}\right]_n \\
&+\left(1-\frac{1}{2\bar{\omega}_n^{(L)}}\right)\left[\vec b^\dag(1)_{[1:L]} O_L(d\lambda^1)\right]_n\Bigr]\end{split}\\
\begin{split}
c^{\prime\dagger}_n:=&\sqrt{\frac{\bar{\omega}_n^{(L)}}{2}}\Bigl[\left(1+\frac{1}{2\bar{\omega}_n^{(L)}}\right)\left[\vec b^\dag(1)_{[1:L]} O_L(d\lambda^1)\right]_n \\
&+\left(1-\frac{1}{2\bar{\omega}_n^{(L)}}\right)\left[O^{\top}_L(d\lambda^1)\vec b(1)_{[1:L]}\right]_n\Bigr]
\end{split}
\end{align}
with $[c'_i,c_j^{\prime\dagger}]=\delta_{i,j},$ $[c'_i,c'_j]=0$ thus achieving the results of theorem \ref{theorem dis 2}.
From the theory of orthogonal polynomials \cite{gaut10}, it is also known that $\omega_n^{(L)}(d\lambda^q)$ are contained in the support interval of $d\lambda^q$ and that the zeros of $P_n(d\lambda^q;x)$ alternate with those of $P_{n+1}(d\lambda^q;x),$ that is
\begin{align}
\begin{split}
&\omega_{n+1}^{(n+1)}(d\lambda^q)<\omega_{n}^{(n)}(d\lambda^q)<\omega_{n}^{(n+1)}(d\lambda^q)<\omega_{n-1}^{(n)}(d\lambda^q)<\ldots \\
&<\omega_{1}^{(n)}(d\lambda^q)<\omega_{1}^{(n+1)}(d\lambda^q)
\end{split}
\end{align}
where $\omega_{i}^{(n+1)},$ and $\omega_{j}^{(n)}$ are ordered in descending order.

Using Eqs. (158) and (161) in \cite{misc14}, we can establish the relation between the field operators $a(f)=\int dx f(x)a_x,$ $a^\dag(f)=\int dx f(x)a_x^\dag$ in Eq. (\ref{the bathtub}) and the operators $c_n,c_n^\dag,$  $c'_n, c^{\prime\dag}_n$.
\begin{align}
\begin{split}
c_n&=\sum_{j=0}^{L-1} \frac{\tilde P_{j}(d\lambda^0;\omega^{(L)}_n(d\lambda^0))}{\sqrt{\sum_{k=0}^{L-1}\tilde P^2_k(d\lambda^0;\omega_n^{(L)}(d\lambda^0))}}\;a(\vec \gamma^0_n), \label{eq:c n}\\
c_n^\dag&=\sum_{j=0}^{L-1} \frac{\tilde P_{j}(d\lambda^0;\omega^{(L)}_n(d\lambda^0))}{\sqrt{\sum_{k=0}^{L-1}\tilde P^2_k(d\lambda^0;\omega_n^{(L)}(d\lambda^0))}}\;a^\dag(\vec \gamma^0_n),
\end{split}\\
\begin{split}
c'_n&=\sqrt{\frac{\bar \omega_n^{(L)}}{2}}\left( \sum_{j=0}^{L-1}\frac{\tilde P_j(d\lambda^1;\omega_n^{(L)}(d\lambda^1))}{\sqrt{\sum_{k=0}^{L-1}\tilde P^2_k(d\lambda^1;\omega_n^{(L)}(d\lambda^1))}}\left( a^\dag(f_j^n)+a(g_j^n)\right) \right),\\
c_n^{\prime\dag}&=\sqrt{\frac{\bar \omega_n^{(L)}}{2}}\left( \sum_{j=0}^{L-1}\frac{\tilde P_j(d\lambda^1;\omega_n^{(L)}(d\lambda^1))}{\sqrt{\sum_{k=0}^{L-1}\tilde P^2_k(d\lambda^1;\omega_n^{(L)}(d\lambda^1))}}\left( a(f_j^n)+a^\dag(g_j^n)\right) \right),\label{eq:d n}
\end{split}
\end{align}
where $f_j^n(x):=\gamma_j^1(x)(1+g(x)/\bar \omega_n^{(L)}),\,$ $
g_j^n(x):=\gamma_j^1(x)(1-g(x)/\bar \omega_n^{(L)}).$
If $x\in {\mathbbm{R}}$,
\begin{equation}
\gamma_n^q(x):=P_n(d\lambda^q;g^{q+1}(x))\sqrt{\frac{J(g^{q+1}(x))}{\pi}\left| \frac{dg(x)}{dx} \right|},
\end{equation}
and can be generalised for higher dimensions by writing Hamiltonian Eq. \eqref{the bathtub} in  \textit{standard form} (see section \ref{The Hamiltonian}). We now turn our attention to $\vec \gamma_0$. We start with deriving its value for the case of theorem \ref{theorem dis 1}. From Eqs. \eqref{eq:gamma 0 general}, \eqref{marticle mapping hamiltonian eq} and recalling the definitions \eqref{eq:particle pos momentum ops def}, we find\footnote{Here the index labeling has been shifted by one to start from 0. This is for consistency with the notation in \cite{misc14}.} $[\vec \gamma_{xx}]_{i,j}=\tr[\hat x_{i-1}(0)\hat x_{j-1}(0)\hat\varrho_0],$ $[\vec \gamma_{xp}]_{i,j}=\tr[\hat x_{i-1}(0)\hat p_{j-1}(0)\hat\varrho_0],$ $[\vec \gamma_{px}]_{i,j}=\tr[\hat p_{i-1}(0)\hat x_{j-1}(0)\hat\varrho_0],$ $[\vec \gamma_{pp}]_{i,j}=\tr[\hat p_{i-1}(0)\hat p_{j-1}(0)\hat\varrho_0],$ $i,j\in{\mathbbm{N}}_+.$
From Eq. (160-161) in \cite{misc14} we have that
\begin{align}
\begin{split}\label{eq:x p expansion proof}
\hat x_n(0)&=\int \mu_0^{(1/2)}(dx) P_{n}(x) (a_x^\dag+a_x),\\
\hat p_n(0)&=\mathrm{i} \int \mu_0^{(1/2)}(dx) P_{n}(x) (a_x^\dag-a_x), \quad n\in {\mathbbm{N}}_0,
\end{split}
\end{align}
thus achieving Eqs. \eqref{eq:gamma mat thorem 1}, \eqref{eq:x p expansion} in theorem \ref{theorem dis 1}.
We achieve the expression for $\| \vec \gamma'_0\|$ in theorem \ref{theorem dis 2} similarly. Namely, from Eq. \eqref{eq:phnon ops} and Eq. (158) in \cite{misc14}, we find
\begin{align}
\begin{split}\label{eq:x p expansion proof 2}
\hat X_n'&=\sqrt{\omega_{max}}\int \mu_0^{(1/2)}(dx) P'_{n}(x^2) \left(a_x^\dag+a_x\right),\\
\hat P_n'&=\frac{\mathrm{i}}{\sqrt{\omega_{max}}} \int \mu_0^{(1/2)}(dx) P'_{n}(x^2)\,x \left(a_x^\dag-a_x\right), \quad n\in {\mathbbm{N}}_0
\end{split}
\end{align}
thus after the definitions $\hat x_n':=\hat X_{n-1}$, $\hat p_n':=\hat P_{n-1}$, $n\in{\mathbbm{N}}_+$ we achieve Eqs. \eqref{eq:gamma mat thorem 2}, \eqref{eq:x p expansion 2} in theorem \ref{theorem dis 2}.
 \qed
\subsection{Multiple chains extension proof}\label{multiple chains extension proof}
Here we will prove corollary \ref{multi cor}.
Let
\begin{equation}
G^{(L_1,L_2,\ldots,L_N)}=
\tr\bigl[\hat{O}\mathrm{e}^{-\mathrm{i} t \hat H^{\vec q}_{\vec L}}\hat{\varrho}_0\mathrm{e}^{\mathrm{i} t\hat H^{\vec q}_{\vec L}}\bigr]
\end{equation}
 where we will denote the scenario that the $nth$ bath has not been discretised, by replacing $L_n$ by $\infty$. We now add and subtract $G$ $N\!-\!1$ times to the r.h.s. of Eq. (\ref{multi eq}) each time discretising one more bath and starting with only one bath discretised. We then apply the triangle inequality and arrive at
\begin{align}
\begin{split}
&\left| \tr\bigl[\hat{O}\mathrm{e}^{-\mathrm{i} t\hat{H}_\textup{mul}}\hat{\varrho}_0\mathrm{e}^{\mathrm{i} t\hat{H}_\textup{mul}}\bigr]-\tr\bigl[\hat{O}\mathrm{e}^{-\mathrm{i} t \hat H^{\vec q}_{\vec L}}\hat{\varrho}_0\mathrm{e}^{\mathrm{i} t\hat H^{\vec q}_{\vec L}}\bigr] \right| \\ &\leq \left|
G^{(\infty,\ldots,\infty,\infty)}
-G^{(\infty,\ldots,\infty,L_N)}  \right|
 +\left| G^{(\infty,\ldots,\infty,L_N)}-G^{(\infty,\ldots,\infty,L_{N-1},L_N)}  \right|\\ &\quad+\ldots+\left| G^{(\infty,L_2,L_3,\ldots,L_N)}
G^{(L_1,L_2,L_3,\ldots,L_N)}
\right|.
\end{split}
\end{align}
In every pair on the r.h.s. of the inequality, there is always one bath which is discretised in one of the $G$ terms but is not discretised for the other $G$ term. We can thus define all the other baths to be part of the system Hamiltonian and then apply theorem \ref{theorem dis 1} or theorem \ref{theorem dis 2} to it. This gives us Eq. (\ref{multi eq}).\qed

