add new chapters

This commit is contained in:
Valentin Boettcher 2022-06-21 10:58:53 +02:00
parent c130b890a6
commit f1f2ddf49b
2 changed files with 915 additions and 0 deletions

551
src/analytical_solution.tex Normal file
View file

@ -0,0 +1,551 @@
\chapter{Comparison with an Analytical Solution}
\section{One Oscillator, One Bath}
\label{sec:oneosc}
\subsection{Model}
\label{sec:model}
The model is given by the quadratic hamiltonian
\begin{equation}
\label{eq:hamiltonian}
H = \frac{Ω}{4}\qty(p^2+q^2) + \frac{1}{2} q
\sum_λ\qty(g_λ^\ast b_λ + g_λ
b^_λ)+\sum_λ\omega_λ b^_λ b_λ,
\end{equation}
where \(a,a^\) are the ladder operators of the harmonic
oscillator, \(q=a+a^\) and \(p=\frac{1}{\iu}\qty(a-a^)\) so
that \([q,p] = 2\iu\).
\subsection{Equations of Motion}
\label{sec:eqmot}
The Heisenberg equation yields
\begin{align}
\dot{q} &=Ω p \label{eq:qdot}\\
\dot{p} &= -Ω q - \int_0^t \Im[α_0(t-s)] q(s)\dd{s} + W(t) \label{eq:pdot}
\\
\dot{b}_λ &= -\iu g_λ \frac{q}{2} - \iu\omega_λ b_λ
\end{align}
with the operator noise
\(W(t)=-\sum_λ \qty(g_λ^\ast b_λ(0)
\eu^{-\iu\omega_λ t } + g_λ b_λ^†(0)
\eu^{\iu\omega_λ t })\),
\(\ev{W(t)W(s)}=α(t-s)\) and \(α_0 \equiv \eval{α}_{T=0}\).
The equation \cref{eq:pdot} arises from
\begin{equation}
\label{eq:bsol}
b_λ(t) = b_λ(0) \eu^{-\iu ω_λ t} - \frac{\iu g_λ}{2}_0^t
q(s) \eu^{-\iu ω_λ (t-s)}\dd{s}.
\end{equation}
The equations for \cref{eq:qdot} and \cref{eq:pdot} can be solved by
finding a matrix \(G(t)\) with \(G(0)=\id\) and
\begin{equation}
\label{eq:eqmotprop}
\dot{G}(t) = A G(t) - \int_0^t K(t-s) G(s)\dd{s},\quad A=\mqty(0 &
Ω \\& 0), \quad K(t)=\mqty(0 & 0\\ \Im[α_0(t)] & 0).
\end{equation}
Then
\begin{equation}
\label{eq:qpsol}
\mqty(q(t)\\ p(t)) = G(t)\mqty(q(0)\\ p(0)) + \int_0^tG(t-s)
\mqty(0\\ W(s))\dd{s}.
\end{equation}
Because we are only interested in solutions for \(t\geq 0\) and the
shape of the convolution in \cref{eq:eqmotprop} the solution may be
found by virtue of the Laplace transform.
Setting
\begin{equation}
\label{eq:laplprop}
\mathcal{L}\{G\}(z) = \int_0^\infty \eu^{-z\cdot t} G(t)
\end{equation}
leads to an algebraic formula
\begin{equation}
\label{eq:galgebr}
\mathcal{L}\{G\}(z) = \qty(z-A + \mathcal{L}\{K\}(z))^{-1}.
\end{equation}
\subsection{Solution}
\label{sec:solution}
We observe that
\begin{equation}
\label{eq:mdef}
M = z-A + \mathcal{L}\{K\}(z) = \mqty(z &\\ Ω +
\mathcal{L}\{\Im[α_0]\}(z) & z)
\end{equation}
and therefore
\begin{equation}
\label{eq:minv}
M^{-1} = \frac{1}{Ω^2 + Ω\mathcal{L}\{\Im[α_0]\}(z) + z^2}
\mqty(z & Ω \\ -(Ω + \mathcal{L}\{\Im[α_0]\}(z)) & z).
\end{equation}
From this we can conclude that \(G_{11}=G_{22}\).
Because \(\ev{W(s)}=0\) for thermal initial states of the bath we have
\begin{equation}
\label{eq:meanvals}
\mqty(\ev{q(t)}\\ \ev{p(t)}) = G(t)\mqty(\ev{q(0)}\\ \ev{p(0)}).
\end{equation}
Knowing this, we can deduce from \(\ev{\dot{q}}= Ω \ev{p}\) that
\begin{align}
\label{eq:onlyoneneeded}
G_{11} &= \frac{\dot{G}_{12}}{Ω} & G_{21} &=\frac{\ddot{G}_{12}}{Ω^2}.
\end{align}
These relations are true independent of the initial state of the
system. It therefore suffices if we concern ourselves with
\(G_{12}\). We nevertheless continue in full generality.
Assume that \(α_0(t)=\sum_{n=1}^N G_n \eu^{-W_n t - \i
\varphi_n}\) with \(W_n=\gamma_n + \i\delta_n\) and \(G_n,
\varphi_n, \gamma_n,\delta_n\in\RR\) for \(t\geq 0\).
This leads to a mathematically simple expression for the Laplace
transform
\begin{equation}
\label{eq:laplace_alpha}
\mathcal{L}\qty{\Im[α_0]}(z) = -\sum_n G_n\qty[\frac{(z+\gamma_n)\sin\varphi_n+\delta_n\cos\varphi_n}{\delta_n^2+(\gamma_n+z)^2}].
\end{equation}
Because \(\mathcal{L}\{\Im[α_0]\}\) appears in the denominator of
\cref{eq:minv} it is desirable to write \cref{eq:laplace_alpha} with a
common denominator. Introducing \(s_n = \sin\varphi_n,\, c_n =
\cos\varphi_n\) and \(z_n= -W_k\) we arrive
at
\begin{equation}
\label{eq:laplace_alpha_better}
\begin{aligned}
\mathcal{L}\qty{\Im[α_0]}(z) &= - \sum_n
G_n\frac{(z_n+\gamma_n)s_n+ \delta_nc_n}{(z-z_n)(z-z_n^\ast)} \\
&= -\frac{\sum_n G_n \qty((z_n+\gamma_n)s_n+
\delta_nc_n)\prod_{k\neq
n}(z-z_k)(z-z_k^\ast)}{\prod_{k}(z-z_k)(z-z_k^\ast)} \\
&= \frac{\sum_n f_n(z)\prod_{k\neq n}(z-z_k)(z-z_k^\ast)}{\prod_{k}(z-z_k)(z-z_k^\ast)}
\end{aligned}
\end{equation}
with the polynomials of first degree
\(f_n(z)=-G_n \qty((z_n+\gamma_n)s_n+\delta_nc_n)\). Because the
above expression is a rational function, the components of
\cref{eq:minv} are rational functions for which the Laplace transform
is particularly simple to invert using the residue theorem. With this
in mind we now calculate
\begin{equation}
\label{eq:prefactorrational}
\frac{1}{Ω^2 + Ω\mathcal{L}\{\Im[α_0]\}(z) + z^2}
% =\frac{\prod_{k}(z-z_k)(z-z_k^\ast)}{\qty[(z+\iΩ)(z-\iΩ)]\prod_{k}(z-z_k)(z-z_k^\ast)
% + \sum_nΩ f_n(z)\prod_{k\neq n}(z-z_k)(z-z_k^\ast)}\\
=\frac{f_0(z)}{p_1(z) + \sum_n q_n(z)}
=
\frac{f_0(z)}{\mu\prod_{n=1}^{N+1}(z-\tilde{z}_l)(z-\tilde{z}^\ast_l)}
= \frac{f_0(z)}{p(z)}
\end{equation}
where
\begin{align}
f_0(z) &= \prod_{k}(z-z_k)(z-z_k^\ast) \\
p_1(z) &= \qty[(z+\iΩ)(z-\iΩ)]\prod_{k}(z-z_k)(z-z_k^\ast) \\
q_n(z) &= Ω f_n(z)\prod_{k\neq n}(z-z_k)(z-z_k^\ast)
\end{align}
and \(\mu\in\RR\). The \(\tilde{z}_l\) are the roots of the real
polynomial
\begin{equation}
\label{eq:polyp}
p(z) = p_1(z) + \sum_{n=1}^{N}q_n(z)
\end{equation}
of degree \(2(N+1)\) where we \textbf{assume that there are
no roots with multiplicity greater than one}.
With this we can now calculate the inverse laplace transform of
expressions of the form \(\frac{f_0(z)g(z)}{p(z)}\) where \(g(z)\) is
any holonome function so that \(\frac{f_0(z)g(z)}{p(z)} \eu^{z\cdot
t}\) falls off fast enough for \(t\geq 0\),
\(\Re(z)>\max_l{\Re(\tilde{z}_l)}=\Delta\) and \(\Re(z) \rightarrow
-\infty\). With this we can close the contour of the inverse Laplace
transform
\begin{equation}
\label{eq:invlap}
\mathcal{L}^{-1}\qty{\frac{f_0(z)g(z)}{p(z)}}(t) =
\frac{1}{2\pi\i}\int_{\Delta - \i\infty}^{\Delta + \i\infty} \frac{f_0(z)g(z)}{p(z)} \eu^{z\cdot
t}\dd{z}
\end{equation}
to the left to obtain
\begin{equation}
\label{eq:simpleinvtrans}
\mathcal{L}^{-1}\qty{\frac{f_0(z)g(z)}{p(z)}}(t)
=
\sum_{l=1}^{N+1}\qty[\frac{f_0(\tilde{z}_l)g(\tilde{z}_l)}{p'(\tilde{z}_l)}
\eu^{\tilde{z}_l \cdot t} + \cc]
\end{equation}
where we assumed that \(g(z)^\ast=g(z^\ast)\) which is the case for
all our purposes. For completeness we give
\begin{equation}
\label{eq:pderiv}
p'(z) = 2\mu\sum_{k=1}^{N+1}\qty[(z-\Re(\tilde{z}_k))\prod_{\substack{n=1\\
n\neq k}}^{N+1}(z-\tilde{z}_n)(z-\tilde{z}^\ast_n)].
\end{equation}
We can immediately conclude that all elements of \(G\) are sums of
exponentials, just like the BCF. In particular
\begin{equation}
\label{eq:gfinal}
G(t) = \sum_{l=1}^{N+1}\qty[R_l \mqty(\tilde{z}_l & Ω \\ \frac{\tilde{z}_l^2}{Ω} & \tilde{z}_l)\eu^{\tilde{z}_l \cdot
t} + \cc]
\end{equation}
with \(R_l={f_0(\tilde{z}_l)}/{p'(\tilde{z}_l)}\).
\subsubsection{Negative Times}
The solution detailed above is only valid for positive times. Because
we strive to employ the same formalism again for negative times, we
will concern ourselves with the transformed quantities
\(\bar{X}(τ) = X(t(τ)) = X(-τ)\) so that \(τ ≥ 0\). It follows that
\(_τ \bar{X}(τ) = \bar{X}'(τ) = -\dot{X}(t(τ))\) so that
\begin{align}
\bar{q}' &= -Ω \bar{p} \label{eq:qtag}\\
\bar{p}' &= Ω \bar{q} + ∫_0^τ \Im[α_0(τ-s)] \bar{q}(s)\dd{s} - \bar{W}(τ) \label{eq:ptag}
\\
\bar{b}'_λ &= \iu g_λ \frac{q'}{2} + \iu\omega_λ b'_λ.
\end{align}
This leads to an equation for \(\bar{G}(τ)\), namely
\begin{equation}
\label{eq:eqmotpropbar}
\bar{G}'(τ) = -A \bar{G}(τ) + \int_0^τ K(τ-s) \bar{G}(s)\dd{s}.
\end{equation}
The solution is obtained from the \(t\geq 0\) case by substituting
\(A\rightarrow -A\) and \(K\rightarrow -K\).
We obtain for \(t\leq 0\)
\begin{equation}
\label{eq:gfinalbar}
\bar{G}(τ) = G(-τ) = G(t) = \sum_{l=1}^{N+1}\qty[R_l \mqty(\tilde{z}_l &\\ -\frac{\tilde{z}_l^2}{Ω} & \tilde{z}_l)\eu^{-\tilde{z}_l \cdot
t} + \cc]
\end{equation}
and
\begin{equation}
\label{eq:qpsolneg}
\mqty(q(t)\\ p(t)) = G(t)\mqty(q(0)\\ p(0)) - \int_t^0 G(t-s)
\mqty(0\\ W(s))\dd{s}.
\end{equation}
\subsection{Applications}
\label{sec:applications}
Knowing \(G\) and \(α\), we can calculate all observables of the
system. Simple closed form expressions of sums of exponentials can be
obtained by using an exponential expansions of \(α\). Throughout,
we assume a thermal bath initial state so that \(\ev{W(t)}=0\).
\subsubsection{Correlation Functions}
\label{sec:correl}
We proceed to calculate \(\ev{q(t)q(s)}\). For brevity we set
\(A=G_{11}\), \(B=G_{12}\), \(p_0=p(0)\) and \(q_0=q(0)\). Then,
\begin{equation}
\label{eq:qcorrel}
\ev{q(t)q(s)} =
\begin{aligned}[t]
& A(t)A(s) \ev{q_0^2} + B(t)B(s) \ev{p_0^2} +
A(t)B(s)\ev{q_0p_0} + B(t)A(s)\ev{p_0q_0} \\
& +\underbrace{_0^t\dd{l}_0^s\dd{r} B(t-l)B(s-r)α(l-r)}_{\equiv Λ(t, s)}.
\end{aligned}
\end{equation}
For a pure harmonic oscillator initial state \(\ket{n}\) we have
\begin{equation}
\label{eq:hoexp}
\begin{aligned}
\ev{q^2_0} &= \ev{p^2_0} = 2n+1 & \ev{q_0p_0} &= \iu.
\end{aligned}
\end{equation}
Note that \(p\) and \(p\) differ from the usual definition by a factor
of \(2\).
\subsubsection{Bath Enery Derivative}
\label{sec:bathflow}
With \cref{eq:qcorrel} we can calculate the time derivative of the
bath energy expectation value
\begin{equation}
\label{eq:bathderiv_1}
\begin{aligned}
\ev{\dot{H}_B} &= ∑_λ ω_λ \qty(\ev{b_λ^\dot{b}_λ} + \cc) \\
&=\frac{1}{4}_0^t\dd{s}\qty[\ev{q(s)q(t)}_λ\abs{g_λ}^2 \eu^{\i
ω_λ(t-s)} + \cc] +
\i_0^t\dd{s} G_{12}(s)∑_λ\abs{g_λ}\bar{n}_λ\qty[\eu^{\i ω_λ s}+\cc]\\
&= -\frac{1}{2}\Im\qty[∫_0^t\dd{s}\ev{q(t)q(s)}\dot{α}_0(t-s)] +
\frac{1}{2}_0^t\dd{s} G_{12}(s)\partial_s\qty[α(s)-α_0(s)] \\
&=
\begin{aligned}[t]
-\frac{1}{2}\Im&\qty[∫_0^t\dd{s}\ev{q(t)q(s)}\dot{α}_0(t-s)] \\&+
\frac12 G_{12}(t)[α(t)-α_0(t)]
-\frac{Ω}{2}_0^t\dd{s} G_{11}(s)\qty[α(s)-α_0(s)],
\end{aligned}
\end{aligned}
\end{equation}
where we've used \(\ev{b_λ(0)}=\ev{b_λ^0}=0\),
\begin{equation}
\label{eq:blambdadotexp}
\begin{aligned}
\ev{b_λ^\dot{b}_λ}= -\i\ev{b_λ^\qty(\frac{g_λ}{2}q + ω_λb_λ)} =
-\i\qty[\frac{g_λ}{2}\ev{b^_λ(t)q(t)} +
\underbrace{ω_λ\ev{b^_λ(t)b_λ(t)}}_{\in\RR\implies\text{cancels }\cc}].
\end{aligned}
\end{equation}
and
\begin{equation}
\label{eq:moreladot}
\begin{aligned}
\ev{b^_λ(t)q(t)} &= \ev{\qty(b(0)^{\dag}_λ\eu^{\i ω_λ t} +
\frac{\i}{2}_0^tg_λ^\ast q(s)\eu^{\i ω_λ (t-s)}\dd{s})q(t)} \\
&= \frac{\i g_λ^\ast}{2}_0^t\ev{q(s)q(t)\eu^{\i ω_λ(t-s)}}\dd{s}
- g_λ^\ast\bar{n}_λ∫_0^t G_{12}(s)\eu^{\i ω_λ s}\dd{s}
\end{aligned}
\end{equation}
with \(\bar{n}_λ=\ev{b^_λ(0)b_λ(0)}\).
For further evaluation of \cref{eq:bathderiv_1} we have to calculate
\begin{equation}
\label{eq:lambdafold}
\begin{aligned}
Λ(t)&=∫_0^t\dd{s}Λ(t,s)\dot{α}_0(t-s)
=∫_0^t\dd{s}_0^t\dd{l}_0^s\dd{r}
B(t-l)B(s-r)α(l-r)\dot{α}_0(t-s)\\
&=∫_0^t\dd{s}_0^t\dd{r}
\begin{aligned}[t]
Θ(s-r)&B(s-r)\dot{α}_0(t-s)\times\\
\biggl[
&_0^{t-r}\dd{u}B(t-r-u)α(u)+∫_0^{r}\dd{u}B(t-r+u)α^\ast(u)
\biggr]
\end{aligned}
\\
&= ∫_0^t\dd{r} g_1(t-r)\qty[g_2(t-r) + g_3(t,r)] = Λ_1(t) + Λ_2(t)
\end{aligned}
\end{equation}
This expression now only uses \(α(t)\) for \(t\geq 0\) so that we can
once again employ the exponential expansion for \(α\). In fact, all
other quantities in \cref{eq:lambdafold} have exponential expansion so
that we can now define\footnote{Note that this is inconsistent with
\cref{sec:solution}.}
\begin{equation}
\label{eq:expansions}
\begin{aligned}
α_0&=∑_k U_k\eu^{-Q_k t} & \dot{α}_0&=∑_k P_k\eu^{-L_k t} & α(t)
&= ∑_nG_n\eu^{-W_n t} \\
A(t) &= ∑_l A_l\eu^{-C_l t} & B(t) &= ∑_l B_l\eu^{-C_l t}.
\end{aligned}
\end{equation}
With this we can calculate,
\begin{align}
\label{eq:lambdaintegrals}
_r^t\dd{s}B(s-r)\dot{α}_0(t-s)
&=\sum_{m,k}\underbrace{\frac{B_mP_k}{L_k-C_m}}_{\equiv
Γ^1_{mk}}\qty[\eu^{-C_m(t-r)}-\eu^{-L_k(t-r)}]=g_1(t-r)\\
_0^{t-r}\dd{u}B(t-r-u)α(u)
&=\sum_{n,l}\underbrace{\frac{B_nG_l}{C_n-W_l}}_{\equiv
Γ^2_{nl}}\qty[\eu^{-W_l(t-r)}-\eu^{-C_n(t-r)}]=g_2(t-r)\\
_0^{r}\dd{u}B(t-r+u)α^\ast(u)
&=\sum_{n,l}\underbrace{\frac{B_nG_l^\ast}{C_n+W_l^\ast}}_{\equiv
Γ^3_{nl}}\qty[\eu^{-C_n(t-r)}-\eu^{-W_l^\ast r-C_n t}]=g_3(t,r)
\end{align}
and
\begin{align}
\label{eq:finalsummands}
Λ_1(t)&= ∑_{m,k,n,l}Γ^1_{mk}Γ^2_{nl}\qty[\frac{1-\eu^{-(C_m+W_l)t}}{C_m+W_l}
-
\frac{1-\eu^{-(C_m+C_n)t}}{C_m+C_n}-
\frac{1-\eu^{-(L_k+W_l)t}}{L_k+W_l}
+
\frac{1-\eu^{-(L_k+C_n)t}}{L_k+C_n}]\\
Λ_2(t)&=
\begin{aligned}[t]
_{m,k,n,l}Γ^1_{mk}Γ^3_{nl}\biggl[\frac{1-\eu^{-(C_m+C_n)t}}{C_m+C_n}
&-\frac{1-\eu^{-(L_k+C_n)t}}{L_k+C_n}
\\&-\frac{\eu^{-(C_n+W_l^\ast)t}-\eu^{-(C_m+C_n)t}}{C_m-W_l^\ast}
+\frac{\eu^{-(C_n+W_l^\ast)t}-\eu^{-(L_k+C_n)t}}{L_k-W^\ast_l}\biggr]
\end{aligned}
\end{align}
Also required for \cref{eq:bathderiv_1} are
\begin{align}
\label{eq:ABconv}
_0^t\dd{s}A(s)\dot{α}_0(t-s) &= ∑_{n,m}\underbrace{\frac{A_nP_m}{L_m-C_n}}_{\equiv
Γ^A_{nm}}\qty[\eu^{-C_n t}-\eu^{-L_m t}]\\
_0^t\dd{s}B(s)\dot{α}_0(t-s) &= ∑_{n,m}Γ^1_{nm}\qty[\eu^{-C_n t}-\eu^{-L_m t}]
\end{align}
and
\begin{multline}
\label{eq:nonzerotemplim}
_0^t\dd{s}A(s)\qty(α(s)-α_0(s)) =\\
_{m,n}\frac{A_nG_m}{C_n+W_m}\qty(1-\eu^{-(C_n+W_m)t}) - ∑_{m,n}\frac{A_nU_m}{C_n+Q_m}\qty(1-\eu^{-(C_n+Q_m)t}).
\end{multline}
This concludes the calculation. A possible measure of simplification
would be to write \cref{eq:bathderiv_1} as a sum of exponentials and
give explicit expressions for the coefficients and exponents. This is
not required for now. Code implementing this can be found under
\url{https://github.com/vale981/hopsflow}.
\section{Two Oscillators, Two Baths}%
\label{sec:twoosc}
The considerations of~\cref{sec:oneosc} can be straight forwardly
generalized to the case of two coupled oscillators coupled in turn to
a bath each.
We will not give explicit formulas for the results in terms of sums of
exponentials, as they are quite extensive and easily obtained via the
use of a computer algebra system or the aforementioned code.
\subsection{Model}
\label{sec:twomodel}
The model is again given by a quadratic hamiltonian
\begin{equation}
\label{eq:hamiltonian}
\begin{aligned}
H &= ∑_{i\in\qty{1,2}} \qty[H^{(i)}_O + q_iB^{(i)} + H_B^{(i)}] + \frac{γ}{4}(q_1-q_2)^2,
\end{aligned}
\end{equation}
where \(H_O^{(i)}= \frac{Ω_i}{4}\qty(p_i^2+q_i^2)\), \(B^{(i)}=\sum_λ\qty(g^{(i),\ast}_λb^{(i)}_λ + g^{(i)}_λ
b^{(i),†}_λ)\) and \(H_B^{(i)}=\sum_λ\omega_λ b^{(i),†}_λ b^{(i)}_λ\).
The \(b^{(i)}\) are the usual bosonic ladder operators of the baths.
The \(a_i^{(i)},a_i^{}\) are the ladder operators of the harmonic
oscillators and \(q_i=a_i+a_i^\) and \(p=\frac{1}{\iu}\qty(a_i-a_i^)\) so
that \([q_i,p_j] = 2\iuδ_{ij}\) and \([q_i,q_j] = [p_i,p_j] = 0\).
\subsection{Equations of Motion}
\label{sec:eqmot_two}
The Heisenberg equation yields
\begin{align}
\dot{q}_i &= Ω_i p_i \label{eq:qidot}\\
\dot{p}_i &= -(Ω_i+γ) q_i - \int_0^t \Im[α^{(i)}_0(t-s)] q_i(s)\dd{s} + W_i(t) \label{eq:pidot}
\\
\dot{b}^{(i)}_λ &= -\iu g^{(i)}_λ \frac{q_i}{2} - \iu\omega^{(i)}_λ b^{(i)}_λ
\end{align}
with the operator noise
\(W_i(t)=-\sum_λ \qty(g_λ^{(i),\ast} b^{(i)}_λ(0)
\eu^{-\iu\omega^{(i)}_λ t } + g_λ^{(i)} b_λ^{(i),†}(0)
\eu^{\iu\omega^{(i)}_λ t })\) satisfying \(\ev{W_i(s)}=0\) and
\(\ev{W_i(t)W_j(s)}=δ_{ij}α^{(i)}(t-s)\). We introduced \(α^{(i)}_0
\equiv \eval{α^{(i)}}_{T=0}\).
We have given most quantities an extra index and accounted for the
coupling between the two oscillators. Apart from this, the equations
of motion have the same structure as in \cref{seq:eqmot}.
Again, we obtain
\begin{equation}
\label{eq:bsoltwo}
b^{(i)}_λ(t) = b^{(i)}_λ(0) \eu^{-\iu ω^{(i)}_λ t} - \frac{\iu g^{(i)}_λ}{2}_0^t
q_i(s) \eu^{-\iu ω^{(i)}_λ (t-s)}\dd{s}.
\end{equation}
We can solve the equations for the \(q_i,\,p_i\)
by finding a matrix \(G(t)\) with \(G(0)=\id\) and
\begin{gather}
\label{eq:eqmotproptwo}
\dot{G}(t) = A G(t) - \int_0^t K(t-s) G(s)\dd{s}\\
A = \mqty(
0 & \Omega & 0 & 0 \\
-\gamma -\Omega & 0 & \gamma & 0 \\
0 & 0 & 0 & \Lambda \\
\gamma & 0 & -\gamma -\Lambda & 0),\;
K(τ) =
\mqty(0 & 0 & 0 & 0 \\
\Im[α^{(1)}_0(t-s)] & 0 & 0 & 0 \\
0 & 0 & 0 & 0 \\
0 & 0 & \Im[α^{(2)}_0(t-s)] & 0),
\end{gather}
where \(Ω=Ω_1\) and \(Λ=Ω_2\) for convenience.
Then
\begin{equation}
\label{eq:qpsol}
\mqty(q_1(t)\\ p_1(t)\\ q_2(t)\\ p_2(t)) = G(t)\mqty(q_1(0)\\ p_1(0) \\ q_2(0)\\ p_2(0)) + ∫_0^tG(t-s)
\mqty(0\\ W_1(s)\\ 0 \\ W_2(s))\dd{s}.
\end{equation}
With the Laplace transform find for \(t\geq 0\) the formula
\cref{eq:galgebr} for the Laplace transform of the solution, albeit
now with a more complicated matrix
\begin{equation}
\label{eq:mdeftwo}
M = z-A + \mathcal{L}\{K\}(z) = \mqty(z & -\Omega & 0 & 0 \\
\mathcal{L}\{\Im[α^{(1)}_0]\}(z)+\gamma +\Omega & z & -\gamma & 0 \\
0 & 0 & z & -\Lambda \\
-\gamma & 0 & \mathcal{L}\{\Im[α^{(2)}_0]\}(z)+\gamma +\Lambda & z)
\end{equation}
that we have to invert.
This can be done easily\footnote{We have use a computer algebra
system. There is probably a pattern to the inverse matrix which can
be found so that the solution for \(N>2\) oscillators may be found.}
and yields
\begin{equation}
\label{eq:minvtwo}
M^{-1}(z) = \frac{1}{\det[M](z)} \tilde{M}(z)
\end{equation}
where \(\tilde{M}\) is a matrix containing only polynomials of \(z\)
and of the Laplace transforms of the bath correlation functions.
The numerator
\begin{equation}
\label{eq:numerator}
\begin{aligned}
\det[M](z)=a& b \Lambda \Omega +a \left(\gamma
\Lambda \Omega +\Lambda ^2 \Omega +\Omega z^2\right)
+b
\left(\gamma \Lambda \Omega +\Lambda \Omega ^2+\Lambda z^2\right)\\
&+\gamma \Lambda ^2 \Omega +\gamma \Lambda
\Omega ^2+\Lambda ^2 \Omega ^2+\gamma \Lambda z^2+\gamma \Omega z^2+\Lambda ^2 z^2+\Omega ^2 z^2+z^4
\end{aligned}
\end{equation}
where \(a=\mathcal{L}\{\Im[α^{(1)}_0]\}\) and \(b=\mathcal{L}\{\Im[α^{(2)}_0]\}\).
Using the same approach as in \cref{sec:solution}, we arrive at an
expression similar to \cref{eq:prefactorrational} for
\((\det[M](z))^{-1}\). The polynomial \(p\) is now of degree
\(4 + 2 \qty(N^{(1)} + N^{(j)})\) where the \(N^{(i)}\) are the number of
terms in the expansions of bath correlation functions for each bath
and the function \(f_0\) now depends on both bath correlation
functions.
We ultimately find that \(G\) is a sum of
exponentials
\begin{equation}
\label{eq:gfinal}
G(t) = \sum_{l=1}^{2+N_1+N_2}\qty[R_l \tilde{M}(\tilde{z}_l)\eu^{\tilde{z}_l \cdot
t} + \cc]
\end{equation}
with \(R_l={f_0(\tilde{z}_l)}/{p'(\tilde{z}_l)}\).
\subsection{Applications}
\subsubsection{Correlation Functions}
\label{sec:correltwo}
We can now proceed to calculate the correlation functions
\(C(t,s) = \ev{x_i(t)x_j(s)}\) where the \(x_i\) are the phase space operators
of the two harmonic oscillators.
We find
\begin{equation}
\label{eq:generalcorr}
C_{ij}(t, s) = G_{ik}(t)G_{jl}(s) C(0,0)_{kl} +
\underbrace{_0^t\dd{l}_0^s\dd{r}G_{ik}(t-l)G_{jl}(s-r) \ev{W_k(l)W_l(r)}}_{_{ij}}.
\end{equation}
The matrix \(Θ_{ij}\) contains the bath-induced correlations and can
be calculated as in the single-oscillator case.
\subsubsection{Bath Enery Derivative}
\label{sec:bathflowtwo}
Similar to the calculations in \cref{sec:bathflow} we find
\begin{equation}
\label{eq:bathderivtwo}
\ev{\dot{H}_B^{(n)}}=-\frac12
\Im_0^tC_{2n-1, 2n-1}(t,s)\dot{α}_0^{(n)}(t-s)\dd{s} + \frac12 ∫_0^t
_{k=1,2}G_{2n-1,2k}_s\qty(α^{(k)}(s)-α_0^{(k)}(s)).
\end{equation}
This can be evaluated using the exponential expansions and yields
yet another sum of exponentials. The steady state flow can then be
found be setting all exponentials to zero.

364
src/hops_tweak.tex Normal file
View file

@ -0,0 +1,364 @@
\chapter{Some Notes on HOPS}
\section{Normalized HOPS}%
\label{sec:norm}
We introduce full HOPS vector \(Ψ = \qty(ψ, φ)\) which can be
decomposed into the zeroth hierarchy order state \(ψ\) and the
non-zero order states \(φ\).
The HOPS equations can then be written in an abstract manner as
\begin{equation}
\label{eq:HOPS}
\begin{aligned}
\dot{ψ} &= F(ψ, φ), & \dot{φ} &= G(ψ, φ),
\end{aligned}
\end{equation}
where \(c\cdot F(ψ, φ) = F(c\cdot ψ, c\cdot φ)\) and
\(c\cdot G(ψ, φ) = G(c\cdot ψ, c\cdot φ)\) for \(c\in\CC\)
The goal is to transform \(ψ \rightarrow \tilde{ψ}\) so that
\begin{equation}
\label{eq:goal}
\norm{\tilde{ψ}} = 1
\end{equation}
in a numerically stable manner.
Introducing the definitions \(\tilde{ψ} = \eu^{f(t)}ψ\) and
\(\tilde{φ} = \eu^{f(t)}φ\) with an
arbitrary\(f\colon \RR \rightarrow \CC\) we can begin to calculate
\begin{equation}
\label{eq:normdgl}
_t\norm{\tilde{ψ}}^2 = \tilde{ψ}^\qty(\dot{f}\tilde{ψ} +
F(\tilde{ψ}, \tilde{φ})) + \cc = \dot{f} \abs{\tilde{ψ}}^2 +
\tilde{ψ}^†F(\tilde{ψ}, \tilde{φ}) + \cc.
\end{equation}
We would now like to obtain \(_t\norm{\tilde{ψ}}^2 = 0\) as well as
\(\dot{f} > 0\) for \(\norm{\tilde{ψ}} < 0\) and \(\dot{f} < 0\) for
\(\norm{\tilde{ψ}} > 0\), so that \(\norm{\tilde{ψ}} = 1\) becomes a
stable fix-point.
Observing \cref{eq:normdgl}, we conclude that our goals can be
achieved by demanding
\begin{equation}
\label{eq:fdgl}
\dot{f} = \frac{\tilde{ψ}^†F(\tilde{ψ},
\tilde{φ})}{\norm{\tilde{ψ}}^2} + g\qty(\norm{ψ}^2).
\end{equation}
The first summand on its own would lead to \(_t\norm{\tilde{ψ}}^2 =
0\) norm conservation. The latter of our goals may be achieved by
choosing \(g(x) = \qty(1-x)\).
These choices lead to an altered HOPS equation
\begin{equation}
\label{eq:normedhops}
\dot{\tilde{Ψ}} = \qty[\frac{\tilde{ψ}^†F(\tilde{ψ},
\tilde{φ})}{\norm{\tilde{ψ}}^2}+\qty(1-\norm{\tilde{ψ}}^2)]\mqty(\tilde{ψ}\\
\tilde{φ}) + \mqty(F(\tilde{ψ},\tilde{φ}) \\ G(\tilde{ψ},\tilde{φ})).
\end{equation}
\section{Multiple Baths}
\label{sec:multibath}
We generalize the NMQSD and HOPS to \(N\) baths for Hamiltonians of
the form
\begin{equation}
\label{eq:multimodel}
H = H_\sys + ∑_{n=1}^N \qty[H_\bath\nth + \qty(L_n^†B_n + \hc)],
\end{equation}
where \(H_\sys\) is the (possibly time dependent) system Hamiltonian,
\(H_B\nth =_λω_λ\nth a_λ^{(n),†}a_λ\nth\),
\(B_n=_{λ} g_λ\nth a_λ\nth\) and the \(L_n={(\vb{L})}_n\) are
arbitrary operators in the system Hilbert space. This models a
situation where each bath couples with the system through exactly one
spectral density and is therefore not fully general.
\subsection{NMQSD}
\label{sec:nmqsd}
Following the usual derivation of the NMQSD \cite{Diosi1998Mar}, we
switch to an interaction picture with respect to the \(H_\bath\)
leading to
\begin{equation}
\label{eq:multimodelint}
H(t) = H_\sys + ∑_{n=1}^N \qty[L_n^†B_n(t) + \hc],
\end{equation}
with \(B_n=_{λ}L_n^† g_λ\nth a_λ\nth\eu^{-\iu ω_λ\nth t}\).
We will discuss the zero temperature case. The finite temperature
methods generalize straight forwardly to multiple baths. Projecting
on a Bargmann (unnormalized) coherent state basis
\(\qty{\ket{\vb{z}^{(1)},\vb{z}^{(2)},\ldots}=
\ket{\underline{\vb{z}}}}\) of the baths
\begin{equation}
\label{eq:projected}
\ket{ψ(t)} = ∫∏_{n=1}^N{\qty(\frac{\dd{\vb{z}\nth}}{π^{N_n}}\eu^{-\abs{\vb{z}}^2})}\ket{ψ(t,\underline{\vb{z}}^\ast)}\ket{\underline{\vb{z}}},
\end{equation}
where \(N_n\) are the number of oscillators in each bath.
We define
\begin{equation}
\label{eq:processes}
η^\ast_n(t) = {\qty(\vb{η}^\ast_t)}_n= -\iu_λg_λ^{(n),\ast} z_λ^{(n),\ast}\eu^{\iu ω_λ\nth t}
\end{equation}
and using
\(\pdv{z_λ^{(n),\ast}}=\dd{s}\pdv{η^\ast_n(s)}{z_λ^{(n),\ast}}\fdv{η^\ast_n(s)}\)
we arrive at
\begin{equation}
\label{eq:multinmqsd}
__t(\vb{η}^\ast_t) = -\iu H ψ_t(\vb{η}^\ast_t) +
\vb{L}\cdot\vb{η}^\ast__t(\vb{η}^\ast_t) - ∑_{n=1}^N L_n^†∫_0^t\dd{s}α_n(t-s)\fdv{ψ_t(\vb{η}^\ast_t)}{η^\ast_n(s)},
\end{equation}
where \(α_n(t-s)= {\qty(\vb{α}(t-s))}_n=_λ\abs{g_λ\nth}^2\eu^{-\iu ω_λ\nth(t-s)}\) are the
zero temperature bath correlation functions. The equation
\cref{eq:multinmqsd} becomes the NMQSD by reinterpreting the
\(\vb{z}\nth\) as normal distributed complex random variables by
virtue of monte-carlo integration of \cref{eq:projected}. The
\(η^\ast_n(t)\) become homogeneous gaussian stochastic processes
defined through
\begin{equation}
\label{eq:processescorr}
\begin{aligned}
\mathcal{M}^\ast_n(t)) &=0, & \mathcal{M}_n(t)η_m(s)) &= 0,
& \mathcal{M}_n(t)η_m(s)^\ast) &= δ_{nm}α_n(t-s).
\end{aligned}
\end{equation}
\subsection{Nonlinear NMQSD}
\label{sec:nonlin}
For the derivation of the lonlinear theory, the characteristic
trajectories of the partial differential equation of motion of
the Husimi-function
\begin{equation}
\label{eq:husimi}
Q_t(\underline{\vb{z}}, \underline{\vb{z}}^\ast) =
\frac{\eu^{-\abs{{\underline{\vb{z}}}}^2}}{π^{_n N_n}}
\braket{ψ(t, {\underline{\vb{z}}})}{ψ(t, {\underline{\vb{z}}}^\ast)}
\end{equation}
have to be determined.
Using \(_{\underline{\vb{z}}}\ket{ψ(t, {\underline{\vb{z}}}^\ast)} =
0\) and \(_{\underline{\vb{z}}^\ast}\bra{ψ(t, {\underline{\vb{z}}})} =
0\) because \(\ket{ψ(t, {\underline{\vb{z}}}^\ast)}\) is holomorphic
we derive
\begin{equation}
\label{eq:husimimotion}
_tQ_t(\underline{\vb{z}}, \underline{\vb{z}}^\ast) = -i
_{n=1}^N\qty[∂_{z_λ^{(n), \ast}}\eu^{-\iu ω_λ\nth
t}\ev{L^_n}_tQ_t(\underline{\vb{z}}, \underline{\vb{z}}^\ast) - \cc],
\end{equation}
where \(\ev{L^_n}_t = \mel{ψ(t, {\underline{\vb{z}}})}{L^_n}{ψ(t,
{\underline{\vb{z}}}^\ast)} / \braket{ψ(t, {\underline{\vb{z}}})}{ψ(t, {\underline{\vb{z}}}^\ast)}\).
The characteristics of \cref{eq:husimimotion} obey the equations of
motion
\begin{equation}
\label{eq:characteristics}
\dot{z}^{(n),\ast}_λ = \iu g_λ\nth \eu^{-\iu ω_λ\nth t} \ev{L^_n}_t
\end{equation}
for the stochastic state labels.
The microscopic dynamics can in-turn be gathered into a shift of the
stochastic processes
\begin{equation}
\label{eq:procshift}
\tilde{η}_n^\ast(t) = η_n^\ast(t) + ∫_0^t\dd{s}α_n^\ast(t-s)\ev{L^_n}_s
\end{equation}
and we obtain the nonlinear NMQSD equation
\begin{multline}
\label{eq:multinmqsdnonlin}
__t(\tilde{\vb{η}}^\ast_t) = -\iu H ψ_t(\tilde{\vb{η}}^\ast_t) +
\vb{L}\cdot\tilde{\vb{η}}^\ast__t(\tilde{\vb{η}}^\ast_t) \\-
_{n=1}^N
\qty(L_n^†-\ev{L^_n}_t)∫_0^t\dd{s}α_n(t-s)\eval{\fdv{ψ_t(\tilde{\vb{η}}^\ast_t)}{η^\ast_n(s)}}_{\vb{η}^\ast(s)
= \vb{η}(\underline{\vb{z}}^\ast(t), s)}.
\end{multline}
The notation
\({\vb{η}^\ast(s) = \vb{η}(\underline{\vb{z}}^\ast(t), s)}\) means
that we replace the microscopic \(z_λ^{(n),\ast}\) in
\cref{eq:processes} with the shifted ones obeying
\cref{eq:characteristics} and evaluate the resulting function at \(s\).
This awkward construction can be remedied by the convolutionless
formulation. It plays no great role in the HOPS formalism.
\subsection{Multi Bath HOPS in Fock-Space Formulation}
\label{sec:multihops}
Following the usual derivation~\cite{RichardDiss} (but with a
different normalization) and using an exponential expansion of the
BCFs \(α_n(τ)=_{\mu}^{M_n}=G_μ\nth\eu^{-W_μ\nth τ}\), we define
\begin{equation}
\label{eq:dops}
D_μ\nth(t) \equiv_0^t\dd{s}G_μ\nth\eu^{-W_μ\nth (t-s)}\fdv{η^\ast_n(s)}
\end{equation}
and
\begin{equation}
\label{eq:dops_full}
D^{\underline{\vb{k}}} \equiv
_{n=1}^N∏_{μ=1}^{M_n}
{\sqrt{\frac{\underline{\vb{k}}_{n,μ}!}{\qty(G\nth_μ)^{\underline{\vb{k}}_{n,μ}}}}
\frac{1}{i^{\underline{\vb{k}}_{n,μ}}}}\qty(D_μ\nth)^{\underline{\vb{k}}_{n,μ}},
\end{equation}
as well as
\begin{equation}
\label{eq:hierdef}
ψ^{\underline{\vb{k}}} \equiv D^{\underline{\vb{k}}}ψ.
\end{equation}
Using
\begin{equation}
\label{eq:commrelation}
[D^\kmat(t),η_n^\ast(t)] = \iu_{μ=1}^{M_n}
\sqrt{\kmat_{n,μ}G\nth_μ} D^{\kmat -
\mat{e}_{n,μ}}
\end{equation}
where \({\qty(\mat{e}_{n,μ})}_{ij}=δ_{ni}δ_{μj}\) we find after some algebra
\begin{multline}
\label{eq:multihops}
\dot{ψ}^\kmat = \qty[-\iu H_\sys + \vb{L}\cdot\vb{η}^\ast -
_{n=1}^N∑_{μ=1}^{M_n}\kmat_{n,μ}W\nth_μ]ψ^\kmat \\+
\iu_{n=1}^N∑_{μ=1}^{M_n}\sqrt{G\nth_μ}\qty[\sqrt{\kmat_{n,μ}} L_^{\kmat -
\mat{e}_{n,μ}} + \sqrt{\qty(\kmat_{n,μ} + 1)} L^_^{\kmat +
\mat{e}_{n,μ}} ].
\end{multline}
The HOPS equations \cref{eq:multihops} can also be rewritten in an
especially appealing form \cite{Gao2021Sep} if we embed the hierarchy
states into a larger Hilbert space using
\begin{equation}
\label{eq:fockpsi}
\ket{Ψ} = \sum_\kmat\ket{\psi^\kmat}\otimes \ket{\kmat}
\end{equation}
where
\(\ket{\kmat}=\bigotimes_{n=1}^N\bigotimes_{μ=1}^{N_n}\ket{\kmat_{n,μ}}\)
are bosonic Fock-states.
Now \cref{eq:multihops} becomes
\begin{equation}
\label{eq:fockhops}
\begin{aligned}
_t\ket{Ψ} &= \qty[-\iu H_\sys + \vb{L}\cdot\vb{η}^\ast -
_{n=1}^N∑_{μ=1}^{M_n}b_{n,μ}^\dag b_{n,μ} W\nth_μ +
\iu_{n=1}^N∑_{μ=1}^{M_n} \sqrt{G_{n,μ}} \qty(b^_{n,μ}L_n +
b_{n,μ}L^_n)] \ket{Ψ}\\
&= \tilde{H}\ket{Ψ}
\end{aligned}
\end{equation}
\section{Estimating the Norms of the Auxiliary States}
\label{sec:normest}
It is possible to find an (semi-rigorous) upper bound to the norms of
the auxiliary states. We will limit ourselves to one bath. The
generalization to multiple baths is straight forward.
Using \cref{eq:fockhops}, we can calculate
\begin{equation}
\label{eq:normdiff}
\begin{aligned}
\iu_t \norm{ψ^{\vb{k}}}^2
&= \bra{Ψ}\ket{k}\bra{k}\tilde{H}\ket{Ψ} - \cc\\
&= \qty^{\vb{k}})^\bra{k}
\qty[-\iu L η^\ast -\iu_{μ=1}^{M}b_{μ}^\dag b_{μ} W_μ
+∑_{μ=1}^{M} \sqrt{G_{μ}} \qty(b^_{μ}L +
b_μ L^†)]\ket{Ψ}- \cc\\
&= \Bigg[-\iu \qty^{\vb{k}})^†L η^\ast ψ^{\vb{k}}
-\iu_{μ=1}^{M}k_μ W_μ \norm{ψ^{\vb{k}}}^2\\
&\phantom{=}\quad -∑_{μ=1}^{M}\qty[\qty^{\vb{k}})^\sqrt{G_{μ}k_μ}^{\vb{k}-\vb{e}_μ} +
\qty^{\vb{k}})^\sqrt{G_{μ}(k_μ+1)}^{\vb{k}+\vb{e}_μ} ]\Bigg] - \cc.
\end{aligned}
\end{equation}
We can now further treat the this expression to find the steady state
norms of the states.
Assuming generically that the term containing the stochastic process
\(η\) vanishes in the time average (as is the case for the steady
state) we will drop it in the following.
Terms of the form \(\Im(ψ^† O φ)\) may be estimated as follows
\begin{equation}
\label{eq:genericest}
\abs{\Im^† O φ)} \leq \norm{ψ} \norm{O φ} \leq \norm{ψ}\norm{O}\norm{φ},
\end{equation}
where the norm on the operator is the standard linear operator norm
\(\norm{O} = \max_{x\in \mathcal{H}}\frac{\ev{O}{x}}{\braket{x}}\).
We now endeavor to find from \cref{eq:normdiff} an estimate of the
steady state norm of \(ψ^{\vb{k}}\). To this end we assume that the
coupling to higher hierarchy states generically lowers the norm and is
therefore neglected. Using \cref{eq:genericest} we can estimate the
influence of the coupling to lower states, choosing the sign
so that the contribution to the norm is positive.
With this we obtain
\begin{equation}
\label{eq:finalest}
_t \norm{ψ^{\vb{k}}}^2 = 0 = -∑_{μ=1}^{M}k_μ \Re[W_μ]
\norm{ψ^{\vb{k}}}^2 +
_{μ=1}^{M}\abs{\sqrt{G_{μ}k_μ}}\norm{ψ^{\vb{k}}}\norm{ψ^{\vb{k}-\vb{e}_μ}}\norm{L}
\end{equation}
and therefore
\begin{equation}
\label{eq:steadynorm}
\norm{ψ^{\vb{k}}} =
\frac{_{μ=1}^{M}\abs{\sqrt{G_{μ}k_μ}}\norm{ψ^{\vb{k}-\vb{e}_μ}}\norm{L}}{_{μ=1}^{M}k_μ \Re[W_μ]}.
\end{equation}
For the nonlinear method, the stochastic process obtains a shift whose
magnitude can be estimated as follows
\begin{equation}
\label{eq:shiftestimate}
\abs{η_{\mathrm{sh}}} \leq \norm{L}_0^\dd{s}\abs{α^\ast(t-s)} \leq
\norm{L} \sum_{μ=1}^M \frac{\abs{G_μ}}{\Re[W_μ]}.
\end{equation}
Assuming the contribution of the shift is norm enhancing, we arrive at
the expression
\begin{equation}
\label{eq:steadynorm_nonlin}
\norm{ψ^{\vb{k}}} =
\frac{_{μ=1}^{M}\abs{\sqrt{G_{μ}k_μ}}\norm{ψ^{\vb{k}-\vb{e}_μ}}\norm{L}}{_{μ=1}^{M}k_μ
\Re[W_μ] - \norm{L}^2 \sum_{μ=1}^M \frac{\abs{G_μ}}{\Re[W_μ]}}.
\end{equation}
This indicates, that we trade hierarchy depth for sample count in the
nonlinear method. Interestingly the divisor in
\cref{eq:steadynorm_nonlin} may vanish, leading to a breakdown of the
estimate. An interpretation would be, that for very strong coupling or
long bath correlation times the interaction \fixme{reference} with the
bath diverges and the method fails. On the other hand, the estimate
may simply be wrong and should account for the coupling to the lower
orders as well.
The relations \cref{eq:steadynorm,eq:steadynorm_nonlin} are recursive
and break off at \(ψ^0\), the norm of which can be assumed to be \(1\)
in the nonlinear method.
These ideas remain to be verified. Especially the assumptions should
be checked. For time dependent coupling, one may maximize the estimate
over all \(L(t)\).
\subsection{Truncation Scheme}
\label{sec:truncsch}
The norm of the \(\vb{k}\)th hierarchy state scales like
\({1} / {\sqrt{\max_μk_μ}}\). This fact in itself, however, is not
too meaningful as the magnitude of the coupling to the lower hierarchy
states is
\begin{equation}
\label{eq:couplingmag}
M_{\vb{k}} = \norm{L} \norm{ψ^{\vb{k}}} \max_μ \abs{\sqrt{G_μk_μ}},
\end{equation}
which balances out the scaling.
Calculating \(M_{\vb{k}}\) explicitly and demanding it to be small
(compared to some energy scale) nevertheless gives a convergent
truncation scheme below a certain coupling strength.
Some basic experimentation has shown, that the cutoff parameter has to
be tuned and is not universally valid.