\documentclass[fontsize=12pt,paper=a4,open=any,
,twoside=false,toc=listof,toc=bibliography,
captions=nooneline,captions=tableabove,english,DIV=16,numbers=noenddot,final]{scrartcl}

\usepackage{../hirostyle}
\addbibresource{../references.bib}
\synctex=1
\title{Calculating heat flows with HOPS}
\author{Valentin Boettcher}
\date{\today}
\begin{document}
\maketitle
\tableofcontents

\section{One Bath}

\subsection{Linear NMQSD, Zero Temperature}
As in~\cite{Hartmann2017Dec} we choose
\begin{equation}
  \label{eq:totalH}
  H = H_S + \underbrace{LB^\dagger + L^\dagger B}_{H_I} + H_B
\end{equation}
with the system hamiltonian \(H_S\), the bath hamiltonian
\begin{equation}
  \label{eq:bathh}
  H_B = \sum_\lambda \omega_\lambda a^\dag a,
\end{equation}
the bath coupling system operator \(L\) and the bath coupling bath
operator
\begin{equation}
  \label{eq:bop}
  B=\sum_{\lambda} g_{\lambda} a_{\lambda}
\end{equation}
which define the interaction hamiltonian \(H_I\).

We define the heat flow out of the system as in~\cite{Kato2015Aug}
through
\begin{equation}
  \label{eq:heatflowdef}
  J = - \dv{\ev{H_B}}{t}.
\end{equation}
Working, for now, in the Schr\"odinger picture the Ehrenfest theorem
can be employed to find
\begin{equation}
  \label{eq:ehrenfest}
  \i\partial_t\ev{H_B} = \ev{[H_B,H]} = \ev{[H_B,H_I]}.
\end{equation}
Thus, we need to calculate
\begin{eqnarray}
  \label{eq:calccomm}
  \begin{aligned}
    [H_B,H_I] &= [H_B, LB^\dag + L^\dag B] \\
    &= L[H_B, B^\dag ] + L^\dag [H_B, B] \\
    &= L[H_B, B^\dag ] - \hc.
  \end{aligned}
\end{eqnarray}
This checks out as the commutator has to be anti-hermitian due to
\cref{eq:ehrenfest}.
Using \([H_B, B^\dag ]=\sum_\lambda \omega_\lambda g^\ast_\lambda
a^\dag_\lambda\) it follows that
\begin{equation}
  \label{eq:expcomm}
  \begin{aligned}
    \ev{[H_B,H_I]} &= \sum_\lambda \omega_\lambda g^\ast_\lambda
    \ev{La^\dag_\lambda} - \cc
    = \sum_\lambda \omega_\lambda g^\ast_\lambda
    \ev{La^\dag_\lambda \eu^{\i \omega t}}_I - \cc\\
    &= \frac{1}{\i}\ev{L\partial_t{\sum_\lambda
        g^\ast_\lambda a^\dag_\lambda \eu^{\i \omega t}}}_I - \cc
    =\frac{1}{\i}\qty(\ev{L\dot{B}^\dag}_I  + \cc)
  \end{aligned}
\end{equation}
where we switched to the interaction picture with respect to \(H_B\)
in keeping with the standard NMQSD formalism.
In essence this is just the Heisenberg equation for \(H_I\). The
expression for \(J\) follows
\begin{equation}
  \label{eq:final_flow}
  J(t) = \ev{L^\dag\partial_t B(t) + L\partial_t B^\dag(t)}_I.
\end{equation}

From this point on, we will assume the interaction picture and drop the
\(I\) subscript.

The two summands yield different expressions in terms of the NMQSD.
For use with HOPS with the final goal of utilizing the auxiliary
states the expression \(\ev{L^\dag\partial_t B(t)}\) should be
evaluated.

We calculate
\begin{equation}
  \label{eq:interactev}
  \ev{L^\dag\partial_t B(t)}=\ev{L^\dag\partial_t B(t)}{\psi(t)} =
  \int \braket{\psi(t)}{z}\mel{z}{L^\dag\partial_tB(t)}{\psi(t)}\frac{\dd[2]{z}}{\pi^N},
\end{equation}
where \(N\) is the total number of environment oscillators and
\(z=\qty(z_{\lambda_1}, z_{\lambda_2}, \ldots)\).
To that end,
\begin{equation}
  \label{eq:nmqsdficate}
  \begin{aligned}
    \mel{z}{\partial_tB(t)}{\psi(t)} &= \sum_\lambda g_\lambda
  \qty(\partial_t \eu^{-\i\omega_\lambda
    t})\partial_{z^\ast_\lambda}\ket{\psi(z^\ast,t)} \\
  &= \int_0^t \sum_\lambda g_\lambda
  \qty(\partial_t \eu^{-\i\omega_\lambda
    t})\pdv{\eta_s^\ast}{z^\ast_\lambda}\fdv{\ket{\psi(z^\ast,t)}}{\eta^\ast_s}\dd{s}\\
  &= -\i\int_0^t\dot{\alpha}(t-s)\fdv{\ket{\psi(z^\ast,t)}}{\eta^\ast_s}\dd{s},
  \end{aligned}
\end{equation}
where \(\eta^\ast_t\equiv -\i \sum_\lambda g^\ast_\lambda
z^\ast_\lambda \eu^{\i\omega_\lambda t}\).
With this we can write
\begin{equation}
  \label{eq:steptoproc}
  \ev{L^\dag\partial_t B(t)} = -\i \mathcal{M}_{\eta^\ast}\bra{\psi(\eta,
    t)}L^\dag\int_0^t\dd{s} \dot{\alpha}(t-s)\fdv{\eta^\ast_s} \ket{\psi(\eta^\ast,t)}.
\end{equation}
Defining
\begin{equation}
  \label{eq:defdop}
D_t = \int_0^t\dd{s} \alpha(t-s)\fdv{\eta^\ast_s}
\end{equation}
as in~\cite{Suess2014Oct} we can write
\begin{equation}
  \label{eq:final_flow_nmqsd}
  J(t) = -\i \mathcal{M}_{\eta^\ast}\bra{\psi(\eta,
    t)}L^\dag\dot{D}_t\ket{\psi(\eta^\ast,t)} + \cc,
\end{equation}
where we've used that the integral in \(D_t\) can be expanded over the
whole real axis. If we assume \(\alpha = \exp(-w t)\) then
\begin{equation}
  \label{eq:hopsj}
  J(t) = \i \mathcal{M}_{\eta^\ast}\bra{\psi^{(0)}(\eta,
    t)}wL^\dag\ket{\psi^{(1)}(\eta^\ast,t)} + \cc.,
\end{equation}
where \(\ket{\psi^{(1)}(\eta^\ast,t)}\) is the first HOPS hierarchy
state. This can be generalized to any BCF that is a sum of exponentials.

Interestingly one finds that
\begin{equation}
  \label{eq:alternative}
  \ev{L\partial_t B^\dag(t)} = \i\int\frac{\dd[2]{z}}{\pi^N}
  \dot{\eta}_t^\ast \mel{\psi(\eta,t)}{L}{\psi(\eta^\ast,t)}.
\end{equation}
However, this seems numerically problematic because \(\eta^\ast\) is not
differentiable in general. The previous expression has the advantage
that we utilize the first hierarchy states that are already being
calculated as a byproduct.

In the language of~\cite{Hartmann2021Aug} we can generalize to
\(\alpha(t) = \sum_i G_i \eu^{-W_i t}\) and thus
\begin{equation}
  \label{eq:hopsflowrich}
  J(t) = \sum_\mu\frac{G_\mu W_\mu}{\bar{g}_\mu} \i\mathcal{M}_{\eta^\ast}\bra{\psi^{(0)}(\eta,
    t)}L^\dag\ket{\psi^{\vb{e}_\mu}(\eta^\ast,t)} + \cc,
\end{equation}
where \(\psi^{\vb{e}_\mu}\) is the \(\mu\)-th state of the first
hierarchy and \(\bar{g}_\mu\) is an arbitrary scaling introduced in
the definition of the hierarchy in~\cite{Hartmann2021Aug} to help with
the scaling of the norm.



\subsection{Nonlinear NMQSD, Zero Temperature}
\label{sec:nonlin}
In the spirit of the usual derivation of the nonlinear NMQSD we write
\begin{equation}
  \label{eq:newb}
  \begin{aligned}
  \ev{L^\dag\dot{B(t)}} &= \int \frac{\dd[2]{z}}{\pi^N} \eu^{-\abs{z}^2}
  \braket{\psi}{z}\!\braket{z}{\psi}
  \frac{\braket{\psi(t)}{z}\!\mel{z}{L^\dag\dot{B(t)}}{\psi(t)}}{\braket{\psi}{z}\!\braket{z}{\psi}}
  \\
  &= \int \frac{\dd[2]{z}}{\pi^N} \eu^{-\abs{z}^2}
  \frac{\mel{z(t)}{L^\dag\dot{B(t)}}{\psi(t)}}{\braket{\psi}{z(t)}\!\braket{z(t)}{\psi}},
  \end{aligned}
\end{equation}
where \(z_{\lambda}^{*}(t)=z_{\lambda}^{*}+\i g_{\lambda} \int_{0}^{t}
\dd{s} \eu^{-\i \omega_{\lambda} s}\ev{L^\dagger}_{s}\).
We find that next steps are the same as in \cref{sec:nonlin} by noting
\begin{equation}
  \label{eq:deriv_trick}
  \eval{\partial_{z^\ast_\lambda}}_{z^\ast=z_\lambda^\ast(t)} =
  \int_0^t\dd{s}\eval{\pdv{\eta^\ast}{z^\ast_\lambda}}_{z^\ast=z^\ast_\lambda(t)}
  \fdv{}{\eta^\ast_s(z^\ast=z^\ast(t))} =
  \int_0^t\dd{s}\eval{\pdv{\eta^\ast}{z^\ast_\lambda}}_{z^\ast=0}
  \fdv{}{\eta^\ast_s(z^\ast=z^\ast(t))},
\end{equation}
which does alter the definition of \(D_t\) but results in the same
HOPS equations.
The shifted process \(\tilde{\eta}^\ast=
\eta^\ast(z^\ast(t),t)=\eta^\ast(t) +
\int_0^t\dd{s}\alpha^\ast(t-s)\ev{L^\dag}_{\psi_s}\) appears directly
in the NMQSD equation but results only in a slight change in the
functional derivative. Note however that
\begin{equation}
  \label{eq:fdvclarification}
  \fdv{}{\eta^\ast_s(z^\ast=z^\ast(t))} \neq \fdv{}{\tilde{\eta}^\ast_s}
\end{equation}
which is not problematic as we have (implicit in~\cite{Diosi1998Mar})
\begin{equation}
  \label{eq:fdvhops}
  \fdv{}{\eta^\ast_s(z^\ast=z^\ast(t))} \ket{\psi(z^\ast)} = \fdv{}{\eta^\ast_s}\ket{\psi(z^\ast(t, z^\ast_0), t)}
\end{equation}
so that the usual HOPS hierarchy follows. Note \(z^\ast_0 = z^\ast(0)\).

Therefore,
\begin{equation}
  \label{eq:newbcontin}
  J(t) =
  -\i
  \mathcal{M}_{\tilde{\eta}^\ast}\frac{\mel{\psi(\tilde{\eta},t)}{L^\dag\dot{\tilde{D}}_t}{\psi(\tilde{\eta}^\ast,t)}}{\braket{\psi(\tilde{\eta},t)}{\psi(\tilde{\eta}^\ast,t)}}
  + \cc,
\end{equation}
where the dependence on \(\tilde{\eta}\) is symbolic and to be
understood in the context of \cref{eq:fdvhops}.

Again we express the result in the language of~\cite{Hartmann2021Aug}
to obtain
\begin{equation}
  \label{eq:hopsflowrich}
  J(t) = \sum_\mu\frac{G_\mu W_\mu}{\bar{g}_\mu}
  \i\mathcal{M}_{\eta^\ast}\frac{\bra{\psi^{(0)}(\eta,
      t)}L^\dag\ket{\psi^{\vb{e}_\mu}(\eta^\ast,t)}}{\bra{\psi^{(0)}(\eta,
      t)}\ket{\psi^{0}(\eta^\ast,t)}} + \cc.
\end{equation}

\subsection{Linear Theory, Finite Temperature}
The finite temperature case needs some additional considerations as
the previous sections dealt explicitly with mean values in a pure
state. The Ehrenfest theorem still holds in mixed states, but we would
like to recover the usual pure state zero temperature formalism. There
are multiple methods for dealing with a thermal initial such as the
thermofield method (see~\cite{Diosi1998Mar}), but because the results
discussed here are to be applied with the HOPS method we shall use the
method described in~\cite{Hartmann2017Dec}.

The shift operator
\begin{equation}
  \label{eq:shiftop}
  \vb{D}(y) = \bigotimes_\lambda \eu^{y_\lambda a_\lambda^\dag-y^\ast_\lambda a_\lambda}
\end{equation}
the ground state of the environment into an arbitrary
coherent state
\begin{equation}
  \label{eq:shiftwork}
  \vb{D}(y)\ket{0} = \ket{y}
\end{equation}
where \(y=(y_1,y_2,\ldots)\) as usual.

This allows us to write the density matrix of the system with a
thermal initial bath as
\begin{equation}
  \label{eq:shiftbath}
  \rho =
  \prod_\lambda\qty(\int\dd[2]{y_\lambda}
  \frac{\eu^{-\abs{y_\lambda}^2\bar{n}_\lambda}}{\pi\bar{n}_\lambda})
  U(t)\vb{D}(y)\ketbra{\psi}\otimes\ketbra{0}\vb{D}(y)^\dag U(t)^\dag.
\end{equation}
The usual step is now to insert \(\id =\vb{D}(y)\vb{D}^\dag(y)\) to
arrive at a new time translation operator
\begin{equation}
  \label{eq:utilde}
  \tilde{U}(t) = \vb{D}^\dag(y)U(t)\vb{D}(y)
\end{equation}
and to interpret the integral in \cref{eq:shiftbath} in a monte-carlo
sense which leads to a stochastic contribution to the system Hamiltonian
\begin{equation}
  \label{eq:thermalh}
  H_{\mathrm{sys}}^{\mathrm{shift}}=L \xi^{*}(t)+L^{\dagger} \xi(t)
\end{equation}
with the stochastic process
\begin{equation}
  \label{eq:xiproc}
  \xi(t):=\sum_{\lambda} g_{\lambda} y_{\lambda} \eu^{-\mathrm{i} \omega_{\lambda} t}
\end{equation}
with corresponding moments \(\mathcal{M}(\xi(t))=0=\mathcal{M}(\xi(t) \xi(s))\) and
\[
\mathcal{M}\left(\xi(t) \xi^{*}(s)\right)=\frac{1}{\pi} \int_{0}^{\infty} \mathrm{d} \omega \bar{n}(\beta \omega) J(\omega) e^{-\mathrm{i} \omega(t-s)}.
\]
Remember that we want to calculate
\begin{equation}
  \label{eq:whatreallymatters}
  \begin{aligned}
    \ev{L^\dag\dot{B}(t)} &= \tr[L^\dag\dot{B}(t)\rho(t)] \\
    &=\prod_\lambda\qty(\int\dd[2]{y_\lambda}
  \frac{\eu^{-\abs{y_\lambda}^2\bar{n}_\lambda}}{\pi\bar{n}_\lambda})\tr[L^\dag\dot{B}(t)
  U(t)\vb{D}(y)\ketbra{\psi}\otimes\ketbra{0}\vb{D}(y)^\dag U(t)^\dag] .
  \end{aligned}
\end{equation}

To recover the zero temperature formulation of this expectation value we
again insert a \(\id\), but have to commute \(\vb{D}(y)^\dag\) past
\(\dot{B}(t)\). This leads to the expression
\begin{equation}
  \label{eq:pureagain}
  \begin{aligned}
    \ev{L^\dag\dot{B}(t)} &=\prod_\lambda\qty(\int\dd[2]{y_\lambda}
    \frac{\eu^{-\abs{y_\lambda}^2\bar{n}_\lambda}}{\pi\bar{n}_\lambda})\\
    &\qquad\times\tr[L^\dag(\dot{B}(t) + \dot{\xi}(t))
    \vb{D}^\dag(y) U(t)\vb{D}(y)\ketbra{\psi}\otimes\ketbra{0}\vb{D}^\dag(y)U(t)^\dag\vb{D}(y)] \\
    &=\prod_\lambda\qty(\int\dd[2]{y_\lambda}
    \frac{\eu^{-\abs{y_\lambda}^2\bar{n}_\lambda}}{\pi\bar{n}_\lambda})\tr[L^\dag\qty{\dot{B}(t) + \dot{\xi}(t)}
    \tilde{U}(t)\ketbra{\psi}\otimes\ketbra{0} \tilde{U}(t)^\dag].
  \end{aligned}
\end{equation}
which returns us to the zero temperature formalism with a transformed
Hamiltonian and the replacement
\begin{eqnarray}
  \label{eq:breplacement}
  B(t) \rightarrow B(t) + \xi(t)
\end{eqnarray}
which plausibly corresponds to the \(L^\dag\) part of \(H_I + H_{\mathrm{sys}}^{\mathrm{shift}}\).

The appearance of \(\dot{\xi}(t)\) may cause concern. However, for
twice differentiable \(\mathcal{M}(\xi(t)\xi^\ast(s))\) the sample
trajectories are smooth.

Alternatively we can calculate
\begin{equation}
  \label{eq:gettingarounddot}
  \begin{aligned}
    \ev{\dot{H}_{\mathrm{sys}}^{\mathrm{shift}}} &=
    \dv{\ev{H_{\mathrm{sys}}^{\mathrm{shift}}}}{t} -
    \frac{1}{\iu}\qty(\ev{H_{\mathrm{sys}}^{\mathrm{shift}}H} -\ev{H
      H_{\mathrm{sys}}^{\mathrm{shift}}}) \\
    &=\dv{\ev{H_{\mathrm{sys}}^{\mathrm{shift}}}}{t} -
    \frac{1}{\iu}\ev{[H_{\mathrm{sys}}^{\mathrm{shift}}, H]} \\
    &=\dv{\ev{H_{\mathrm{sys}}^{\mathrm{shift}}}}{t} -
    \frac{1}{\iu}\ev{[H_{\mathrm{sys}}^{\mathrm{shift}}, H_I]}.
  \end{aligned}
\end{equation}


Now,
\begin{equation}
  \label{eq:hshcomm}
  [H_{\mathrm{sys}}^{\mathrm{shift}}, H_I] = \xi(t) [L^\dag, L]
  B(t)^\dag + \xi^\ast(t) [L, L^\dag] B
\end{equation}
and therefore
\begin{equation}
  \label{eq:finalex}
  \ev{[H_{\mathrm{sys}}^{\mathrm{shift}}, H_I]} = -i \mathcal{M}_{\eta^\ast}\mel{\psi}{\xi(t)^\ast[L,L^\dag]D_t}{\psi}.
\end{equation}
This is an expression that we can easily evaluate with the HOPS
method.

\printbibliography{}
\end{document}




%%% Local Variables:
%%% mode: latex
%%% TeX-master: t
%%% End:
