\section{Introduction}

The fractional calculus (FC) is a very important tool associated with several problems which appear in physics, engineering an other sciences \cite{podlubny,gorenflo,Debnath,Yusufoglu}. As an important example we mention  diffusion processes, particularly  to the fractional partial Burgers' equation (FPBE) appearing in the traffic flow and gas  dynamics  \cite{Wang}. On the order hand, the variational iteration method (VIM) is a relatively  new approaches to provide an analytical approximation to linear and nonlinear problems \cite{Odibat,Bal,Momami}. Those authors  consider the VIM  applied to the time-fractional partial Burgers' equation.\\ Here we are interested in  the VIM associated with the so-called generalized FPBE

\begin{equation}
{}_{c}{\sf D}_t^{\alpha} [u]= u_{xx} +Au^p u_{x},\quad u=u(x,t),
\end{equation}
\\
with ${\sf D}=\frac{\partial}{\partial t}, 0<\alpha\leq 1$ and  $1<\alpha\leq 2, A\in\mathbb{R}$ and $p>0$.\\

\noindent
Some particular results appearing in the literature, are recovered.

\noindent
The paper is organized as follows: in section 2, preliminaries  will be presented, particularly, a short review involving the FC, specifically the Caputo derivative  and its respective  Laplace transform, after of VIM it is analysed, in particular for nonlinear fractional partial differential equation. In section 3 we will present the VIM and the Burgers' equation and provide a lemma, with the proof involving the general case and showing some particular cases of the parameters  presented  as examples. In the section 4 we will discuss the  approximate  solutions for the generalized FPBE, recovering the results involving the classical Burgers' equation, presenting several examples with the respective graphics. Concluding remarks close the paper.

\section{Preliminaries}

In this section we will present  definitions and results that we use in the paper, a short review of FC and the VIM for the linear and nonlinear equations.

\subsection{Fractional Calculus}
First of all, we introduce the Riemann-Liouville (RL) fractional integral, considered
in the left, only \cite{Diethel}. Let $n$ be a positive integer and $\alpha \in \mathbb{C}$ such that
Re$(\alpha) > 0$, we define the RL fractional integral by means of

\begin{equation}
{\sf I}_t^{\alpha} [f(t)] := \frac{1}{\Gamma(n - \alpha)} \int_0^t
f(\tau){(t-\tau)^{\alpha + 1 - n}}{\mbox{d}}\tau\ ,\quad n-1 < \alpha < n.
\end{equation}

\noindent
The Caputo derivatives has been used by many authors in several physical applications \cite{Mainar,Silva,Oliveira,Figueiredo,Diethelm,Liu,Wazwaz}.
One reason for this choice is the fact that the initial conditions associated with the fractional differential equation
are usually expressed in terms of integer order derivatives. Let $n$ be a positive
integer and $\alpha \in \mathbb{C}$ such that Re$(\alpha) > 0$. We introduce the fractional
derivative of order $\alpha$ in the Caputo sense, denoted by ${}_{c}{\sf D}_t^{\alpha}[f(t)]$,
by means of the integral
\begin{equation}
{}_{c}{\sf D}_t^{\alpha} [f(t)] = \left\{
\begin{array}{l}
\displaystyle \frac{1}{\Gamma(n - \alpha)} \int_0^t \frac{f^{(n)}(\tau)}{(t-\tau)^{\alpha+ 1 - n}}{\mbox{d}}\tau\,,
\quad n-1 < \alpha < n,\\
\\
\displaystyle {\sf D}^n f(t) \equiv \frac{{\mbox{d}}^n}{{\mbox{d}}t^n} f(t)\,, \quad \alpha = n.
\end{array} \right.
\end{equation}
The relation involving the Caputo derivatives and the RL fractional integral is given by
\begin{equation}
{}_{c}{\sf D}_t^{\alpha}[f(t)] = {\sf I}_t^{n-\alpha} {\cdot} _{c}{\sf D}_t^n f(t)
\end{equation}
with the $n - 1 < \alpha < n$.

\noindent
As we have already said, the Laplace transform methodology is an efficient tool to discuss a fractional
differential equation. Then, we introduce the Laplace transform of the derivative in the Caputo sense. Denoting
by $\mathscr{L}$ the Laplace integral operator, we can write the Laplace transform of the
Caputo derivatives as follows
\begin{equation}
(\mathscr{L}[{}_{c}{\sf D}_t^{\alpha} f(t)])(s) = s^{\alpha}(\mathscr{L}[f(t)])(s) - \sum_{k=0}^{n-1}
s^{\alpha - k - 1} ({\sf D}^k f)(0^{+})\,,
\end{equation}
with $s$ is the parameter of the Laplace transform.

\noindent
As we have mentioned above, this expression shows that the Laplace transform of the fractional
derivative in the Caputo sense involves only the derivative of integer order evaluated in $t=0^{+}$,
conversely the corresponding RL derivative. In our particular problems, as we will be
seen in the sequence, we take the parameter $\alpha$ as a real number such that $0 < \alpha \leq 1$ in
problems involving (anomalous) diffusion and $1 < \alpha \leq 2$ in problems associated with
wave propagation.

\noindent
To close this subsection, as an example, we consider the fractional integral and the fractional
derivative of a power function $t^{\lambda}$, with $\lambda$  a real parameter. For the fractional
integral we have
\begin{equation}
{\sf I}_t^{\alpha} [t^{\lambda}] = \frac{\Gamma(\lambda + 1)}{\Gamma(\lambda + 1 + \alpha)} t^{\lambda + \alpha}
\end{equation}
with $\alpha \geq 0$, $\lambda > -1$ and $t > 0$. On the other hand, for the fractional derivative we have
\begin{equation}
{}_{c}{\sf D}_t^{\alpha}  [t^{\lambda}] = \frac{\Gamma(\lambda + 1)}{\Gamma(\lambda + 1 - \alpha)} t^{\lambda - \alpha}
\end{equation}
with $\alpha > 0$, $\lambda > -1$ and $t > 0$.

\subsection{Variational iteration method }

The VIM \cite{He1,He2,He3} was extended to fractional differential equations and has been one of  the methods frequently used. Classical and fractional differential equations are studied using  VIM \cite{He2,He3}. On the order hand classical and  fractional partial differential equations are studied in \cite{Mainardi,Waz1,Waz2,Ozer,Sheng,Ruiz,Hayat,Shah,Chen} in particular nonlinear dynamics for local fractional Burgers' equation arising in fractal flow is discussed in \cite{Tenreiro}.

\noindent
Here we consider a more general fractional differential equation
\begin{center}
\quad ${}_{c}{\sf D}_t^{\alpha} [u]+R[u]+N[u]=f(t)$,
\end{center}
\noindent
where ${}_{c}{\sf D}_t^{\alpha} [u]$ is the Caputo derivative, $R[u]$ is a linear term, $N[u]$ is a nonlinear one and $f(t)$ is a function associated with the non homogeneous term. Odibat and Momami in \cite{odibat} applied the VIM to the above equation and suggested a variational iteration formula

\begin{equation*}
\left\{
\begin{array}{l}
\displaystyle u_{n+1}=u_n+\int_0^t\lambda(t,\tau)(_{c}{\sf D}_t^{\alpha}u_n+R[u]+N[u]-f(\tau)){\mbox{d}}\tau,\qquad 0<t\\
       \\
\displaystyle \lambda(t,\tau)=-1,\qquad 0< \alpha\leq 1\\
\\
\displaystyle \lambda(t,\tau)=\tau-t,\qquad 1< \alpha\leq 2.
\end{array} \right.
\end{equation*}

\noindent
$\lambda(t,\tau)$ are known as the Lagrange multipliers associated with the variational iteration formula, this multipliers are evaluated with the general theory of Lagrange multipliers \cite{Inokuti}.

\section{VIM and fractional Burgers' equation}
In this section, we will present the VIM applied to the Burgers' equation:
\begin{equation}
{}_{c}{\sf D}_t^{\alpha} [u]= u_{xx} +Au^p u_{x},\quad u=u(x,t)
\end{equation}
with ${\sf D}=\dfrac{\partial}{\partial t}, 0<\alpha\leq 1$ and  $1<\alpha\leq 2, \quad A\in\mathbb{R}$ and $p>0$.\\\\

\noindent
The so-called correction functional  for Eq.(8) is

\begin{equation}
u_{k+1}(x,t)= u_k(x,t)+\frac{1}{\Gamma (\beta)}\int_{0}^{t}(t-\tau)^{\beta-1}\quad\lambda(t)\left(\dfrac{\partial^{\alpha}}{\partial t^{\alpha}}u_k-\dfrac{\partial^{2}}{\partial x^{2}}\bar{u}_{k}-A\bar{u}_{k}^{p}\dfrac{\partial}{\partial x}\bar{u}_k \right) {\mbox{d}}\tau,
\end{equation}

\noindent
where $\beta=\alpha+1-m$ and $m-1<\alpha\leq{m}$.\\
\noindent
The Eq.(9)  can be approximately expressed by means of
\begin{equation}
u_{k+1}(x,t)=u_k(x,t)+\int_{0}^{t}\lambda(t) \left(\dfrac{\partial^{m}}{\partial \tau^{m}}u_k-\dfrac{\partial^{2}}{\partial x^{2}}\bar{u}_{k}-A\bar{u}_{k}^{p}\dfrac{\partial}{\partial x}\bar{u}_{k}  \right){\mbox{d}}\tau.
\end{equation}
If  we use integration by parts and remembering that the stationary term in the functional  $\delta\bar{u}_{k}=0$, we get three cases as follow:\\

\noindent
a) For $m=1$ we have $0 < \alpha \leq 1$ and the correction functional  can be approximately by means of\\
\begin{equation*}
\delta u_{k+1}(x,t)=\delta u_{k}(x,t)+\delta \int_{0}^{t}\lambda(t)\left(\dfrac{\partial}{\partial t}u_k-\dfrac{\partial^{2}}{\partial x^{2}}\bar{u}_{k}-A\bar{u}_{k}^{p}\dfrac{\partial}{\partial x}\bar{u}_k\right){\mbox{d}}\tau.
\end{equation*}

\noindent
Thus,  we have the system

\begin{equation}
\begin{split}
1+\lambda(\tau)&=0\\
\lambda^{\prime}(\tau)&=0,\\
\end{split}
\end{equation}

\noindent
whose solution is $\lambda(\tau)=-1$. Then,  we obtain the following formula (note that, here we have $\beta=\alpha$)

\begin{equation}
u_{k+1}(x,t)=u_{k}(x,t)-\frac{1}{\Gamma(\alpha)}\int_{0}^{t}(t-\tau)^{\alpha-1}\left(\dfrac{\partial^{\alpha}}{\partial t^{\alpha}}u_k-\dfrac{\partial^2}{\partial x^2} u_{k}-A u_{k}^{p}\dfrac{\partial}{\partial x}{u}_{k} \right){\mbox{d}}\tau.
\end{equation}

\noindent
b) For $ m =2$ we have $1 <\alpha\leq 2$ and the correction functional  can be approximately expressed by means of\\
\begin{equation*}
\delta u_{k+1}(x,t)=\delta u_{k}(x,t)+\delta \int_{0}^{t}\lambda(t)\left(\dfrac{\partial^2}{\partial t^2}u_k-\dfrac{\partial^2}{\partial x^2}\bar{u}_{k}-A\bar{u}_{k}^{p}\dfrac{\partial}{\partial x}\bar{u}_{k}\right) {\mbox{d}}\tau.
\end{equation*}\\

\noindent
Thus, we get the system:
\begin{equation}
\begin{split}
\lambda^{\prime\prime}(\tau)=0,\\
\lambda(\tau)_{\tau=t}=0,\\
1-\lambda^{\prime}(\tau)_{\tau=t}=0,
\end{split}
\end{equation}

\noindent
whose solution can be written as $\lambda(\tau)=\tau-t$. Then,  for $ m = 2 $ and $1 <\alpha\leq 2$, we obtain the following interaction formula with $\beta=\alpha-1$

\begin{equation}
u_{k+1}(x,t)=u_k(x,t)-\frac{\alpha-1}{\Gamma(\alpha)}\int_{0}^{t}(t-\tau)^{\alpha-1}\left(\dfrac{\partial^{\alpha}}{\partial t^{\alpha}}u_k-\dfrac{\partial^2}{\partial x^2} u_{k}-A u_{k}^{p}\dfrac{\partial}{\partial x}{u}_{k} \right){\mbox{d}}\tau.
\end{equation}

\noindent
c) News Lagrange multipliers with $0<t$ and $0<\alpha$. In this case, our main result,  we propose a lemma with its proof, recover the two precedent results and present an example. \\

\noindent
\begin{lem} If the correction functional of the{\rm{ Eq.(9)}} is given by the Riemann integration

\begin{equation}
u_{n+1}=u_{n}+\int_{0}^{t}\lambda(t,\tau)\left[ _c{\sf D}_t^{\alpha}(u_{n})-(\bar{u}_{n})_{xx}-A(\bar{u}_{n})^{p}(\bar{u}_{n})_{x}\right]{\mbox{d}}\tau,
\end{equation}

\noindent
with $0 <t$,  $0<\alpha$ and $(\bar{u}_{n})_{xx}$, $A(\bar{u}_{n})^p(\bar{u}_{n})_x$  are restrictions of the variations of the functional associated with {\rm{ Eq.(15)}}, then  the Lagrange multiplier is

\begin{equation*}
\lambda(t,\tau)=\frac{(-1)^{\alpha}(\tau-t)^{\alpha-1}}{\Gamma(\alpha)}.
\end{equation*}\\\\
\end{lem}

\textbf{Proof.} First of all, we transform Eq.(15) in its integral form and taking the Laplace  transform on both sides of  the new equation
\begin{equation}
\mathscr {L}[u_{n+1}]=\mathscr {L}[u_n+{\sf I}_t^{\alpha}\lambda_{RL}(t,\tau)(_{c}{\sf D}_{\tau}^{\alpha}u_n-(\bar{u}_n)_{x x}-A(\bar{u}_n)^{p}(\bar{u}_n)_x)]
\end{equation}
where $\lambda_{RL}(t,\tau)$ is the Lagrange multiplier of the integral form for  Eq.(15).  We consider
\begin{equation}
{\sf I}_t^{\alpha}\lambda_{RL}(t,\tau)_{c}{\sf D}_t^{\alpha}u_n=\frac{1}{\Gamma(\alpha)}\int_{0}^t(t-\tau)^{\alpha-1}\lambda(t,\tau)_c {\sf D}_t^{\alpha}u_n(\tau){\mbox{d}}\tau.
\end{equation}
Thus, if  $\lambda_{RL}(t,\tau)=\lambda(t-\tau)$, Eq.(17) is the convolution of the function
\begin{equation}
a(t)=\frac{\lambda(t)t^{\alpha-1}}{\Gamma(\alpha)}
\end{equation}
\noindent
and  $_c {\sf D}_t^{\alpha}u_n(t)$.
The terms $(\bar{u}_n)_{x x}$ and $A(\bar{u}_n)^p(\bar{u}_n)x$ which are considered restrictions on variations, implying $\delta(\bar{u}_n)_{x x}=0$ and $\delta A(\bar{u}_n)^p(\bar{u}_n)x=0$.

\noindent
Using the variational functional associated with Eq.(16) and the Laplace transform of the Caputo derivative, we obtain

\begin{equation*}
\delta\mathscr {L}[u_{n+1}]=\delta\mathscr {L}[u_{n}]+\delta\mathscr {L}[{\sf I}_t^{\alpha}\lambda_{RL}(_c {\sf D}_t^{\alpha}u_n-(\bar{u}_n)_{x x}-A(\bar{u}_n)^p(\bar{u}_n)x)],\\
\end{equation*}

\begin{equation*}
=\delta\mathscr {L}[u_{n}]+\delta[\mathscr {L}[a(s)]s^{\alpha}\mathscr {L}[u_n(s)]]-\delta\sum_{k=0}^{m-1}u^k(0^+)s^{\alpha-1-k},\\
\end{equation*}

\begin{equation}
= (1+ \mathscr {L}[a(s)] s^{\alpha})\delta\mathscr {L}[u_n],
\end{equation}\\

\noindent
where, using the Euler-Lagrange equation associated with Eq.(16), we have
\begin{equation*}
1+\mathscr{L}[a(s)]s^{\alpha}=0,
\end{equation*}

\noindent
resulting

\begin{equation}
\mathscr{L}[a(s)]=-\frac{1}{s^\alpha}.
\end{equation}
 Performing the inverse Laplace transform of Eq.(20) we get
\begin{equation}
a(t)=-\frac{t^{\alpha-1}}{\Gamma(\alpha)}.
\end{equation}

\noindent
Comparing Eq.(18) and Eq.(21) we have $\lambda_{RL}(t,\tau)=-1$ and using Eq.(16) we can write

\begin{equation*}
\displaystyle
u_{n+1}=u_n-{\sf I}_t^{\alpha}[_c {\sf D}_t^{\alpha}u_n-(\bar{u}_n)_{xx}-A(\bar{u}_n)^p(\bar{u}_n)_{x}],\\
\end{equation*}
\begin{equation*}
\displaystyle
{u_{n+1}=u_n}-{\displaystyle{\int_{0}^{t}}}\displaystyle{\frac{(t-\tau)^{\alpha-1}}{\Gamma(\alpha)}} {{\left[_{c}{\sf D}_t^{\alpha}u_n-(\bar{u}_n)_{xx}-A(\bar{u}_n)^p(\bar{u}_n)_{x} \right]}} {\mbox{d}}\tau,\\\\
\end{equation*}
\begin{equation*}
\displaystyle
u_{n+1}=u_{n}+\int_{0}^{t}\frac{(-1)^{\alpha}(\tau-t)^{\alpha-1}}{\Gamma(\alpha)} \left[ _{c}{\sf D}_{t}^{\alpha}(u_{n})-(\bar{u}_{n})_{xx}-A(\bar{u}_{n})^{p}(\bar{u}_{n})_{x} \right]\mbox {d}\tau.
\end{equation*}

\noindent
Finally, using Eq.(17) we obtain the following general expression
\begin{equation}
\lambda(t,\tau)=\frac{(-1)^{\alpha}(\tau-t)^{\alpha-1}}{\Gamma(\alpha)}.
\end{equation}

\noindent
In what follows, we recover the known results and discuss an example.

\begin{corol}
In {\rm{Eq.(22)}} if $\alpha=1$ then $\lambda(t,\tau)=-1$ and if $\alpha=2$ then $\lambda(t,\tau)=\tau-t$.\\
\end{corol}

\textbf{Proof.}
Direct substitution $\alpha=1$ ($\alpha=2$) in Eq.(22).\\\\

\noindent
\textbf {Example 1:} Analogously to the results that can be obtained by [8], we have the approximated solution for the FPBE when $0 <\alpha\leq 1$ using the initial condition $u(x,0)= g(x)$. The problem to be consider is
\begin{equation}
\left\{
\begin{array}{l}
\displaystyle _c{\sf D}_t^{\alpha}u+u\frac{\partial u}{\partial x}-\frac{\partial^2 u}{\partial x^2}=0,\quad t>0,\quad 0<\alpha\leq{1}, \\\\
\displaystyle u(x,0)=g(x), \quad 0\leq x \leq 1.
\end{array} \right.
\end{equation}

\noindent
Using {\bf Corollary 1} we have $\lambda(t,\tau)=-1$ and  substitution  in the correction functional Eq.(15) we  obtain
 $$u_{n+1}(x,t)=u_n(x,t)-{\sf I}_t^{\alpha}\left[\dfrac{\partial^{\alpha}}{\partial t^{\alpha}}u_n-u_n(u_n)_x-(u_n)_{x x}\right].$$

\noindent
The correction functional can be written as
\begin{equation}
u_{n+1}(x,t)=u_n(x,t)-\frac{1}{\Gamma(\alpha)}\int_{0}^{t}(t-\tau)^{\alpha-1}[_c{\sf D}_t^\alpha u_n+u_n(u_n)_x-(u_n)_{x x}]{\mbox{d}}\tau.
\end{equation}

\noindent
Using Eq.(24) we obtain the following approximations

\begin{equation*}
 \displaystyle
 u_0= u(x,0)=g(x),
\end{equation*}

 \noindent
\begin{align*}
 \displaystyle
 u_1= g(x)-[g g^{\prime}-g^{\prime\prime}]{\displaystyle\frac{{t^{\alpha}}}{\Gamma(\alpha+1)}},
\end{align*}
\begin{equation*}
 \displaystyle
 u_2= u_{1}+(2g(g^{\prime})^2+g^2g^{\prime\prime}-2g g^{(3)}-4g^{\prime}g^{\prime\prime}+g^{(4)}){{\frac{t^{2\alpha}}{\Gamma(1+2\alpha)}}}
\end{equation*}
\noindent
\begin{equation*}
\qquad\qquad\qquad\qquad\qquad\qquad -(gg^{\prime}-g^{\prime\prime})[(g^{\prime})^2+g g^{\prime\prime}-g^{(3)}]{{\frac{\Gamma(1+2\alpha)}{\Gamma^2(1+\alpha)}\frac{t^{3\alpha}}{\Gamma(1+3\alpha)}}}.
\end{equation*}

\noindent
Similar calculus are using for obtain others approximations. To close the section we prove another corollary.

\noindent
\begin{corol}
With $\lambda(t,\tau)=\displaystyle{\frac{(-1)^\alpha(\tau-t)^{\alpha-1}}{\Gamma(\alpha)}}$ and correction functional as
\begin{equation}
u_{n+1}(x,t)=u_n(x,t)+\frac{1}{\Gamma(\beta)}\int_{0}^{t}(t-\tau)^{\beta-1}\lambda(t,\tau)[_c{\sf D}_t^{\alpha}-(u_n)_{x x}-A(u_n)^p(u_n)_{x}]{\mbox{d}}\tau
\end{equation}
where $\beta=\alpha-[\alpha]$ and $[\alpha]$ is integer part of $\alpha$, then the iteration equation for $0<\alpha\leq 1$ and $1<\alpha\leq 2$, is\\
\begin{equation}
u_{n+1}=u_n-{\sf I}_t^{\alpha}[_c{\sf D}_t^{\alpha}-(u_n)_{x x}-A(u_n)^p(u_n)_{x}],
\end{equation}\\
and
\begin{equation}
u_{n+1}=u_n-(\alpha-1){\sf I}_t^{\alpha}[_c{\sf D}_t^{\alpha}u_n-(u_n)_{x x}-A(u_n)^p(u_n)_{x}].
\end{equation}

\noindent
respectively.
\end{corol}

\noindent
\textbf{Proof.}
By the {\bf Corollary 1}, $\lambda(t,\tau)=-1$ for $0<\alpha\leq{1}$ therefore $\beta=\alpha$, which imply

\begin{equation}
u_{n+1}(x,t)=u_n(x,t)-\frac{1}{\Gamma(\alpha)}\int_{0}^{t}(t-\tau)^{\alpha-1}[_c{\sf D}_t^\alpha u_n-A(u_n)^p(u_n)_x-(u_n)_{x x}]{\mbox{d}}\tau,
\end{equation}

\noindent
or in the following form

\begin{equation}
u_{n+1}(x,t)  =  u_n-{\sf I}_t^{\alpha}[_c{\sf D}_t^{\alpha}u_n-(u_n)_{x x}-A(u_n)^p(u_n)_x].
\end{equation}

\noindent
Now for $1<\alpha\leq 2$ therefore $\beta=\alpha-1$, by the {\bf Corollary 1}, we have $\lambda(t,\tau)=\tau-t$ and using the relation

\begin{align*}
\frac{1}{\Gamma(\alpha-1)}=\frac{-(\alpha-1)}{\Gamma(\alpha)},
\end{align*}

\noindent
we obtain
\begin{equation*}
{\displaystyle{u_{n+1}(x,t)  =  u_n+\frac{1}{\Gamma(\alpha-1)}\int_{0}^{t}(t-\tau)^{\alpha-1}[_c{\sf D}_t^{\alpha}u_n-(u_n)_{x x}-A(u_n)^p(u_n)_x]{\mbox {d}}\tau}},
\end{equation*}

\noindent
or in the following form

\begin{equation}
\noindent
{\displaystyle{u_{n+1}(x,t) =  u_n-(\alpha-1){\sf I}_t^{\alpha}[_c{\sf D}_t^{\alpha}u_n-(u_n)_{x x}-A(u_n)^p(u_n)_x]}}.
\end{equation}

\section{Approximate solutions for FBPE}

\noindent
There are some researchers that consider fractional Burgers' equation to model the diffusion behaviour of the flow through porous medium. In this section we will consider three examples of the fractional Burgers' equation, two of them with $0<\alpha\leq{1}$ and another one for the case $1<\alpha\leq{2}$, with its respective initial conditions. Here $u$ is the flow of velocity, the viscosity coefficient is consider equal to $-1$, and without  loss of generality we take $A=-1$ and $p=1$ in the Eq.(1). Note that, the viscosity coefficient  corresponds the second derivative in the Eq.(1). \\

\noindent
For those three examples, we consider an approximation of the solution in the following  series form

\begin{center}
$u(x,t)= u_{0}(x,t)+u_{1}(x,t)+u_{2}(x,t)+\cdots$
\end{center}

\noindent
To make the graphics, we stop the series in the three term. \\

\noindent
We mention  that the effect of the fractional derivative recover the memory effect associated with physical phenomena for   $0<\alpha\leq{1}$,  $1<\alpha\leq{2}$ \cite{Deng}. In particular with suitable initial conditions, when $A=0$, $\alpha=1$ and $\alpha=2$ in Eq.(1), we  recover the memory effect associated with a heat equation and wave equation, respectively. Also, with $\alpha=1$ in Eq.(1), we  obtain  the classical  Burgers' equation.

\noindent
To close this section we will present some examples with graphics to see better the effect involving  the parameter associated with the derivative.\\\\

\noindent
\textbf{Example 2:}
Consider the problem

\begin{equation}
\left\{
\begin{array}{l}
\displaystyle _c{\sf D}_t^{\alpha}u+u\frac{\partial u}{\partial x}-\frac{\partial^2 u}{\partial x^2}=0,\quad t>0,\quad 0<\alpha \leq {1}, \\

\\
\displaystyle u(x,0)=g(x)=\sin(\pi x) , \quad 0<x\leq{1}.
\end{array} \right.
\end{equation}\\

\noindent
Using Eq.(25) with $\lambda(t,\tau)=-1$ we obtain

\begin{equation*}
{\displaystyle{u_0(x,t)= \sin(\pi x)}},
\end{equation*}

\begin{equation*}
{\displaystyle{u_1(x,t)= \sin(\pi x)-\pi \sin(\pi x)[\cos(\pi x)+(\pi)^2]\frac{t^2}{\Gamma(\alpha+1)}}},
\end{equation*}

\begin{equation*}
\displaystyle
u_2(x,t)= u_{1}(x,t)+(2\pi^2 \sin(\pi x)\cos^2(\pi x)-\pi^2 \sin^4(\pi x)+2\pi^3 \sin(\pi x) \cos(\pi x)
\end{equation*}

\begin{equation*}
\qquad\qquad+4\pi^3 \cos(\pi x)\sin(\pi x)+\pi^4\sin^4(\pi x))
{\displaystyle{\frac{t^{2\alpha}}{\Gamma(1+2\alpha)}}}
\end{equation*}

\begin{equation*}
\qquad\qquad\qquad\qquad\qquad -\pi \sin(\pi x)\left[\cos(\pi x)+(\pi)^{2}\right]\pi^2 \cos(\pi x)-\pi^2 \sin^3(\pi x)
\end{equation*}

\begin{equation*}
\qquad\qquad \quad\qquad\qquad\qquad\qquad\qquad\qquad-\pi^4 \sin(\pi x))\displaystyle{\left(\frac{\Gamma(1+2\alpha)t^{3\alpha}}{\Gamma^2(1+\alpha)\Gamma(1+3\alpha)}\right)}.
\end{equation*}

\noindent
The graphics in Figure 1, with $\alpha=0.2$ and Figure 2 with $\alpha=1$ elucidate the approximation of $u(x,t)$.

\vspace{1 cm}

\begin{figure}[H]
    \centering
    \includegraphics [scale=0.6] {figure1}
    \caption{$u(x,0)=\sin(\pi x)$, with  $\alpha=0.2$.}
    \label{figRotulo}
  \end{figure}

\begin{figure}[H]
    \centering
    \includegraphics [scale=0.6] {figure2}
    \caption{$u(x,0)=\sin(\pi x)$, with  $\alpha=1$.}
    \label{figRotulo}
  \end{figure}

\noindent
\textbf{Example 3:} Consider the  Example 2 with the initial condition given by $ u(x,0)=g(x)={\displaystyle{\frac{e^x}{1+e^x}}}$.
Using Eq.(25) with $\lambda(t,\tau)=-1$ we get

\begin{equation*}
\displaystyle{u_{0}(x,t)  =  \frac{e^x}{1+e^x}},
\end{equation*}

\begin{equation*}
{\displaystyle{u_1(x,t) =  \frac{e^x}{1+e^x}-\left(\frac{e^{x}(e^{x}+e^{2x}-1)}{(1+e^{x})^3}\right) \left(\frac{t^\alpha}{\Gamma(\alpha+1)}\right)}},
\end{equation*}
\\\\

\noindent
\begin{equation*}
{\displaystyle{u_2(x,t) =\frac{e^x}{1+e^x}- \left(\frac{e^{x}(e^{x}+e^{2x}-1)}{(1+e^{x})^3}\right)   \left(\frac{t^\alpha}{\Gamma(\alpha+1)}\right)}}
\end{equation*}

\begin{equation*}
\qquad\qquad\qquad +\left(\displaystyle{\frac{e^{x}(1-12e^{x}-29e^{2x}+4e^{3x}+37e^{4x}+50e^{5x}+11e^{6x}+2e^{7x})}{(1+e^{x})^8}} \right)
\end{equation*}

\begin{equation*}
\qquad\qquad\qquad\qquad\times\left(\frac{t^\alpha}{\Gamma(1+2\alpha)}\right)-
\left(\displaystyle{\frac{e^{x}(1-4e^{x}-10e^{2x}+6e^{3x}+2e^{4x}+e^{5x})}{(1+e^{x})^6} }\right)
\end{equation*}

\begin{equation*}
\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\quad\times\left(\displaystyle{\frac{\Gamma(1+2\alpha)}{\Gamma^2(1+\alpha)}\frac{t^{3\alpha}}{\Gamma(1+3\alpha)}}\right).
\end{equation*}

\noindent
Also here, the graphics in Figure 3, with $\alpha=0.2$  and Figure 4, with $\alpha=1$, elucidate the approximation of $u(x,t)$.

\begin{figure}[H]
    \centering
    \includegraphics [scale=0.5] {figure3}
    \caption{$u(x,0)=\displaystyle \frac{e^x}{1+e^x}$, with $\alpha=0.2$.}
    \label{figRotulo}
\end{figure}
\begin{figure}[H]
    \centering
    \includegraphics [scale=0.5] {figure4}
    \caption{$u(x,0)=\displaystyle \frac{e^x}{1+e^x}$, with $\alpha=1$.}
    \label{figRotulo}
\end{figure}

\noindent
\textbf{Example 4:}
Consider the problem

\begin{equation*}
\left\{
\begin{array}{l}
\displaystyle _c{\sf D}_t^{\alpha}u+u\frac{\partial u}{\partial x}-\frac{\partial^2 u}{\partial x^2}=0,\quad t>0,\quad 1< \alpha \leq 2, \\
       \\
\displaystyle u(x,0)=\sin(\pi x),\qquad u_t(x,0)=0,\quad 0\leq x \leq 1.
\end{array} \right.
\end{equation*}\\

\noindent
In this case, using Eq.(25) with $\lambda(t,\tau)=\tau-t$ we  obtain

\begin{equation*}
\noindent
{\displaystyle{u_0(x,t)=\sin(\pi x)}},
\end{equation*}
\begin{equation*}
\noindent
{\displaystyle{u_1(x,t)  =  \sin(\pi x)-\alpha(\alpha-1)[\pi \sin(\pi x)\cos(\pi x)-\pi^2 \sin(x)]\left[\frac{t^\alpha}{\Gamma(\alpha+1)}\right]}},
\noindent
\end{equation*}
\begin{equation*}
{\displaystyle{u_2(x,t)  =  \sin(\pi x)-\alpha(\alpha-1)[\pi \sin(\pi x)\cos(\pi x)-\pi^2 \sin(x)]\left[\frac{t^\alpha}{\Gamma(\alpha+1)}\right]}}
\end{equation*}
\begin{equation*}
\qquad\qquad+\alpha(\alpha-1)[\sin(\pi x)(2\pi^2 \cos^2(\pi x)-\pi^2 \sin^3(\pi x)+6\pi^3 \cos(\pi x)
\end{equation*}
\begin{equation*}
\qquad\qquad\qquad+\pi^4 \sin^3(\pi x)]\left[\frac{t^2\alpha}{\Gamma(1+2\alpha)}\right]-\alpha(\alpha-1)(\pi \sin(\pi x))(\cos(\pi x)+x^2)
\end{equation*}
\begin{equation*}
\qquad\qquad\qquad\quad\times(\pi^2 \cos(\pi x)-\pi^2 \sin^3(\pi x)-\pi^4 \sin(\pi x)\left[\frac{\Gamma(1+2\alpha)t^3\alpha}{\Gamma^2(1+\alpha)\Gamma(1+3\alpha)}\right].
\end{equation*}

\noindent
The graphics in Figure 5, with $\alpha=1.2$ and Figure 6, with $\alpha=2$ one can see the evolution of the approximation of $u(x,t)$.

\begin{figure}[H]
    \centering
    \includegraphics [scale=0.5] {figure5}
    \caption{$u(x,0)=\sin(\pi x),u_{t}(x,0)=0$, with $\alpha=1.2$.}
   \label{figRotulo}
  \end{figure}

 \begin{figure}[H]
    \centering
    \includegraphics [scale=0.5] {figure6}
    \caption{$u(x,0)=\sin(\pi x),u_{t}(x,0)=0$, with $\alpha=2$.}
   \label{figRotulo}
  \end{figure}

\section{Conclusions}

The FC is very useful in the recuperation of the memory of phenomena, by used of fractional derivative in  time variable. News Lagrange multipliers for FPBE are  identified by Laplace transform and particular cases are recovered. Using this multipliers and the VIM we obtained approximations of the solutions for FPBE taking three terms only. Then, we conclude that the VIM is a powerful and efficient technical to approximate solutions to FPBE.

\section*{References}


\section{Introduction}
The ability of describing images with sentences has numerous applications like helping the visually impaired independently browse the Internet. Recently, with competitions like Microsoft COCO's Image Captioning \cite{coco}, there has been an increased interest in the task of automatic image description generation \cite{google}. With this interest, there is a dire need for large scale datasets that can be used for training these sentence generation models. Datasets like COCO \cite{coco} and Flickr30M \cite{flickr} have been collected by crowdsourcing the description task to human workers on Amazon Mechanical Turk. Once the sentences are generated by one crowd worker, both \cite{coco, flickr} send their sentences to additional crowd workers to verify the accuracy of the sentences. The biggest bottleneck in growing these datasets to a much larger scale has been the cost of generating these sentences and verifying their accuracy.

Humans are strikingly proficient at "filling in the blanks" --- whether it be crosswords, hangman, or Wheel of Fortune. We enjoy partially filled puzzles because of the feeling of simultaneously knowing and not knowing the full answer \cite{crossword}. Previous research by von Ahn and Dabbish show gamification to be an effective vehicle for labeling images in datasets for free \cite{gwap1}. However, such games were only limited to single word annotations \cite{gwap1}. This paper explores another game mechanism in order to generate more complex, full-sentence annotations. We propose a gamification method to crowdsource sentence annotations of images by having players write sentences for other players to guess. Additionally, we show that there exists a direct correlation between the accuracy of the sentence description and a player's ability to guess the sentence.

Motivated to reduce the cost of collecting a large image captioning dataset, we present a game that achieves the following: \begin{enumerate}\item Generates sentence descriptions of images \item Verifies that these sentences are accurate \item Captures sentences of better quality than those collected by Amazon Mechanical Turk (AMT). \end{enumerate}
\begin{figure}[t]
\centering
\includegraphics[width=3.3in]{game3}
\caption{A screenshot of SentenceRacer's interface. The left side displays the image and the state of the verified words so far. The right side displays the chat interface and the scoreboard for guessing.}
\label{fig-example2}
\end{figure}
\section{System Description}
SentenceRacer is played with a minimum of three players. Each round of the game rotates a leader position. The leader sets a sentence describing the image for other players to guess. After eliminating stop words from the sentence, we allow all players to see guesses made by other players while blocking the leader's communication with the guessers. Players only have limited time to guess the words set by the leader. A correct guess rewards both the guesser and the leader with points and reveals the guess' position in the sentence. This reward system implicitly motivates the leader to write descriptive sentences about the image, as they will be easily guessed by the other players.

\section{Data and Analysis}
To gather the data, we took ten groups of four volunteers and ran each group on the same ten randomly sampled images from Microsoft's COCO dataset \cite{coco}.
\subsection{People Find SentenceRacer to be more fun}

Qualitative results suggest that SentenceRacer is more fun in comparison to the task of image captioning. Surveys comparing SentenceRacer and a standard AMT image captioning task show that players find SentenceRacer more fun and engaging particularly because of the social and fast-paced aspects of the game.

\subsection{SentenceRacer's Sentences are Confirmed by AMT}
Sentences collected by SentenceRacer were sent to AMT for verification by three crowd workers. A sentence is verified by AMT if at least two out of the three workers agree that the sentence accurately describes the image. A sentence is considered verified by SentenceRacer if all the words in the sentence can be guessed by the players. We found that 87.8\% of sentences verified by SentenceRacer were also verified by AMT, while only 54.9\% of the sentences not verified by SentenceRacer were verified by AMT. We also found that the sentences collected from SentenceRacer have a higher percentage of verified sentences (87.8\%) than those collected from AMT workers (85.5\%) on the same images.

We investigated the relation of the percentage of sentences verified with the number of remaining blanks left in the game. Table 2 shows that the percentage verified increases as the number of blanks decreases. The tail end of the blanks is sparse, causing the data to have high variance. However, we believe that this trend still shows that SentenceRacer's verification process adequately determines whether a sentence accurately describes an image. The number of blank spaces are directly correlated with how likely a sentence will be verified.

\begin{table}[h!]
\centering
\begin{tabular}{|c|c|c|c|c|}
\hline
Source & Total Blanks & \# Sentences & Verified (\%) \\ \hline
& 4 & 7 & 42.80 \\
& 3 & 12 & 50.00 \\
SentenceRacer & 2 & 9 & 33.30 \\
& 1 & 12 & 75.00 \\
& 0 & 49 & \textbf{87.80} \\ \hline \hline
AMT & - & 200 & 85.50 \\ \hline
\end{tabular}
\caption{Correlation between number of blanks and percentage of verified sentences. As the number of blanks decreases, the percentage of sentences verified by AMT increases. Also in comparison, sentences collected from AMT have a lower verification percentage than the sentences collected by SentenceRacer.}
\end{table}

\subsection{SentenceRacer has Higher Sentence Quality than AMT}

Figure 2 shows the quality of some sentences we received from both tasks on AMT and from playing SentenceRacer. We measure sentence quality by the amount of information we can extract from the sentence describing an image. The average number of objects, object attributes, and pairwise-object relationships per sentence is a basic indicator of sentence quality \cite{scenegraph}.
\begin{figure}[h!]
\centering
\includegraphics[width=3.3in]{table4}
\caption{Comparing verified sentences from AMT and SentenceRacer .}
\label{fig-example}
\end{figure}
Table 2 shows that SentenceRacer has statistically significant more objects and relationships and may suggest that SentenceRacer provides more attributes as well. We believe SentenceRacer's sentence quality stems from the rule that correct guesses reward both the guesser and the leader. Players are incentivized to write longer sentences, leading to higher averages of objects, relationships, and attributes, than AMT tasks, where this incentive is absent.

\begin{table}[h!]
\centering
\begin{tabular}{|C{2.5cm}|c|c|c|}
  \hline
& Objects & Relationships & Attributes\\ \hline
\textbf{AMT}\ \ \ \ \ \ \ \ \ \ \ (n=200) & 2.30 & 1.02 & 1.17 \\ \hline
\textbf{SentenceRacer} (n=49) & 2.98 & 1.88 & 1.45 \\ \hline
P-values & $<$ 0.01 & $<$ 0.001 & 0.1 \\ \hline
\end{tabular}
\caption{T-test showing that increased number of objects, relationships, and (potentially) attributes.}
\end{table}

\section{Conclusion}
In this paper, we demonstrate how SentenceRacer is able to collect sentences describing images, an expensive task for Computer Vision research. This system introduces the idea of collecting and verifying sentences through a game that uses contextual cues as a means of entertainment and verification. Our evaluations suggest that this game is more enjoyable than standard methods of collecting sentences. SentenceRacer can also simultaneously perform the collection and verification of sentences. Finally, we show that the sentences collected by SentenceRacer are of higher quality than those collected by AMT.

We hope to explore how a list of taboo words may increase the diversity of the sentences collected. We also hope to investigate creative ways of using the waiting period between rounds to possibly attempt crowdsourcing other tasks such as grounding objects in the sentence within the image itself.

\balance

\bibliographystyle{acm-sigchi}


\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