Dephasing dynamics of an impurity coupled to an anharmonic environment

TL;DR

This paper investigates how anharmonic environments influence the dephasing and information back-flow of an impurity, revealing temperature-dependent effects and limitations of Gaussian approximations.

Contribution

It demonstrates the impact of anharmonicity on dephasing dynamics and non-Markovian behavior, highlighting differences from Gaussian environment assumptions.

Findings

High temperatures cause strong dephasing.

Low temperatures induce significant information back-flow.

Gaussian approximation poorly models low-temperature dynamics.

Abstract

We analyze the dephasing dynamics of an impurity coupled to an anharmonic environment. We show that a strong anharmonicity produces two different effects depending on the environment temperature: for high temperatures, the system suffers a strong dephasing, while for low temperatures there is a strong information back-flow (as measured by the Breuer-Laine-Piilo (BLP) non-Markovianity measure). Both dephasing and back-flow are particularly significant when the anharmonic potential allows environment states very close to the dissociation limit. In contrast, the information back-flow is suppressed when assuming the environment to be Gaussian. In this regard, we find that the Gaussian approximation is particularly poor at low temperatures.

Equations

$$H_S={\omega }_s{\sigma }_z,$$

$$H_k=\frac{p_k^2}{2m}+D_k\left(e^{-2{\alpha }_k{x}_k}-2e^{-{\alpha }_k{x}_k}\right),$$

$$H_I=S\otimes B,$$

$$B_k$$

$$E_{kn}:=-\frac{{\omega }_k\Lambda }{2}-\frac{{\omega }_k}{2\Lambda }{\left(n+\frac{1}{2}\right)}^2+{\omega }_k\left(n+\frac{1}{2}\right).$$

$$H_k$$

$$B_k$$

$${\rho }_s(t)={\text{Tr}}_E\{{\mathcal{U}}^{-1}(t,t_0){\rho }_s(0)\otimes {\rho }_B^{\beta }(0)\mathcal{U}(t,t_0)\},$$

$$\mathcal{U}(t,t_0)$$

$$+$$

$$C^l(0)={\text{Tr}}_E\{B^l{\rho }_E\}=\sum _{k_1,\dots ,k_l}{\text{Tr}}_E\{B_{k_1}\cdots B_{k_l}{\rho }_E\}.$$

$$\alpha (t-s)={\text{Tr}}_E\{\tilde{B}(t)\tilde{B}(s){\rho }_E^{\beta }\},$$

$$C^0$$

$$C(t)$$

$$\underset{\epsilon \to 0^{+}}{\lim }⟨m|x|m⟩$$

$$-$$

$$\underset{\epsilon \to 0^{+}}{\lim }$$

$$\times$$

$${\omega }_k$$

$$g_k$$

$$J(\omega ):=\Theta \left(2{\omega }_c-\omega \right)\eta \frac{\omega }{{\omega }_c}{e}^{-\frac{\omega }{{\omega }_c}},$$

$$H$$

$$e^{iHt}$$

$$\rho (t)$$

$$\displaystyle\rho_{s}(t)=\Phi(t)\left[\rho_{s}(0)\right]=\left(\begin{array}[]{cc}\rho_{s}^{11}(0)&\chi(t)\rho_{s}^{12}(0)\\[5.0pt] \chi^{*}(t)\rho_{s}^{21}(0)&\rho_{s}^{22}(0)\end{array}\right)$$

$$\chi (t)=e^{2i{\omega }_{s}t}\prod _k{\text{tr}}_k\left(e^{-i{H}_k^{-}t}{\rho }_k^{\beta }{e}^{i{H}_k^{+}t}\right),$$

$$\chi (t)=e^{2i{\omega }_{s}t-\Gamma (t)},$$

$$\Gamma (t)$$

$$=$$

$$H_k^{+}$$

$$H_k^{-}$$

$$e^{\pm i{H}_k^{+}t}$$

$$e^{\pm i{H}_k^{-}t}$$

$${\rho }_0=\frac{1}{2}\left(11+\frac{1}{2}{\sigma }_x\right).$$

Full text

I Introduction

An accurate description of hybrid systems such as molecular ensembles containing different types of degrees of freedom is still challenging. An often convenient approach is to consider some of these freedoms as an open system, while others that operate at a different time scale are treated as an environment and described by its statistical properties Breuer and Petruccione (2002); Brandes (2005); Rivas and Huelga (2011); Rivas et al. (2014); Breuer et al. (2016); de Vega and Alonso (2017). The open system may represent for instance an electronic component while the environment describes a set of nuclei. In other situations the system is an impurity coupled to a complex environment that can be either in condensed matter form, a liquid or a gas. One of the most successful approaches consists in approximating such environment as a set of harmonic oscillators, an idea that was first put forward by Feynman and Vernon Feynman and Vernon Jr. (1963). In this context, the best known model is the spin-boson Leggett et al. (1987), which is very extended to describe energy transfer between two or more molecules (or molecular states) in the presence of an environment Gilmore and McKenzie (2006); Nazir (2009). Examples of these systems are antenna molecules within photosynthetic complexes that are coupled to a protein environment Xu and Schulten (1994); Plenio and Huelga (2008); Ishizaki and Fleming (2009); Rebentrost et al. (2009); de Vega (2011); Mohseni et al. (2013), or electron-transfer reactions between electronic donor and acceptor states that are conditioned by the motion of nuclear degrees of freedom Tang (1993); Barbara et al. (1996); Mavros et al. (2016). In both cases, the protein and nuclear degrees of freedoms are represented by a harmonic bath producing dephasing. Other situations where electron transfer is affected by an environment include transport in polymers Brédas et al. (2004), dynamics of organic molecules within solar cells Zhao et al. (2012), or impurities in strongly condensed matter or liquids Tanimura and Mukamel (1993).

However, the emergence of new scenarios and the rapid development of experimental techniques such as time-resolved nonlinear spectroscopy are giving more information on the dynamics of complex molecular systems, therefore requiring more detailed models and analysis Okumura and Tanimura (1996). For instance, when the electronic dynamics occur in the presence of nonpolar liquids or low-frequency intramolecular modes (such as torsional motion) the harmonic approximation is known to fail Tang (1993); Evans (2000); Lockwood et al. (2002); Wang and Thoss (2004, 2007); Thoss and Wang (2006). Also, nuclear environments may have a complex structure where some modes are highly anharmonic, while the remaining ones can be treated as harmonic and considered to be linearly coupled to the anharmonic ones Garg et al. (1985); Hsieh and Kapral (2014). Another relevant type of environment that is present in condensed matter physics and quantum technological devices corresponds to a set of spins or spin bath Prokof’ev and Stamp (2000). For these spin environments it is well known that the statistics (i.e. the behaviour of different order correlation functions of the environment coupling operators), is highly non-Gaussian, and therefore their description in terms of a spectral density is highly inaccurate. In such context of quantum technological devices, where noise-optimized quantum control is required, an accurate characterization of the environment and its non-Gaussian features is essential Norris et al. (2016); Sung et al. (2019).

In this paper, we consider an impurity coupled to an environment consisting of a set of independent Morse oscillators López-López et al. (2011); Hsieh and Cao (2018) producing dephasing. Advantageously, the Morse oscillator Hamiltonian is analytically solvable Morse (1929), giving rise to a spectrum of discrete energy levels, or bound states, and a set of scattering states. Such scattering states form a continuum that produces an additional decoherence and is describable as a harmonic bath Burkey and Cantrell (1984); Kazansky (1997). Here, we are not interested in a complete description of the whole Morse spectrum, but rather we focus in its discrete part. This allows us to make a systematic study of the effects of the anharmonicity on the open system dynamics: by tuning the anharmonicity of the potential we will be able to vary from a regime where each environment oscillator has a large number of equally spaced bound states, thus a harmonic limit, to a highly anharmonic situation with only two bound states. In other words, our model allows us to describe the dephasing dynamics of an impurity coupled to a harmonic bath and to a spin bath, as well as all the intermediate anharmonic regimes. As shown in Fig. (1) the Hamiltonian describes a star configuration, with the impurity coupled to a set of independent oscillators, such that the harmonic and spin limits correspond to the pure dephasing version of the spin-boson and central spin models, respectively.

The paper is organized as follows. We first discuss general concepts of the model in Sec. II as well as the dynamical map that describes the reduced dynamics of the impurity. Then we divide the analysis in three parts: In Sec. III we discuss the changes in the impurity dephasing time, that is intimately related to the persistence of quantum mechanical properties such as coherence and entanglement along the evolution, when considering two impurities initially entangled Yu and Eberly (2002). In Sec. IV we consider the ratio between the flow of information inwards and outwards of the impurity Cosco and Maniscalco (2018), which determines the amount of information that is lost from the system into the environment. In turn, such ratio is strongly linked to the presence of BLP non-Markovianity Rivas et al. (2014); Breuer et al. (2016); Li et al. (2018), a concept that has gained increasing attention in various contexts including quantum information and quantum metrology, biological systems, ultra-cold atoms or quantum thermodynamics de Vega and Alonso (2017). In Sec. V we analyze the degree in which the impurity dynamics can be well described with a Gaussian map, such that it only depends on the environment second order moment or correlation function. This will allow us to explore when will the non-Gaussianity of the environment become more or less relevant depending on the environment temperature and the anharmonicity parameter. Finally, in Sec. VI we draw some conclusions and outlook.

II A spin coupled to a Morse oscillator environment

We consider a system with a total Hamiltonian $H=H_{S}+H_{E}+H_{I}$. This corresponds to an impurity with Hamiltonian ($\hbar=1$)

[TABLE]

where $\sigma_{z}$ is a spin ladder operator that can be written in terms of the impurity internal basis $|\pm\rangle$ as $\sigma_{z}=|+\rangle\langle+|-|-\rangle\langle-|$, coupled to an anharmonic environment composed of independent Harmonic oscillators with a Hamiltonian $H_{E}=\sum_{k}H_{k}$ ($k\in\{1,\dots,40\}$), where

[TABLE]

is describing each $k$-th oscillator in terms of a Morse potential. Here we have defined the position operator as $x_{k}=\sqrt{\alpha_{k}/2m\omega_{k}}(b_{k}+b_{k}^{\dagger})$, where $b_{k}$ ($b_{k}^{\dagger}$) the standard harmonic annihiliation (creation) operators. The depth and the width of the Morse potential is determined by two different parameters, $D_{k}$ and $\alpha_{k}$, respectively. Notice that this represents a specific case of the right panel of Fig.1 in which the potential part $V(x_{k})=D_{k}\left(e^{-2\alpha_{k}x_{k}}-2e^{-\alpha_{k}x_{k}}\right)$ is centered at the origin and equal for both impurity levels. Finally, the coupling is described through the Hamiltonian

[TABLE]

where $S:=\sigma_{z}$ and $B$ are system and environment coupling operators, respectively. Just like for the free Hamiltonian we can write $B=\sum_{k}B_{k}$, with

[TABLE]

with $g_{k}$ a constant that determines the coupling strength of the $k$-th oscillator to the system. Thus, we consider the coupling as linear and therefore proportional to the position operator (also called reaction coordinate in the context of electron transfer) of each Morse oscillator $B_{k}\sim x_{k}$ Okumura and Tanimura (1996); Wang and Thoss (2004); Thoss and Wang (2006); Okumura and Tanimura (1996); Wang and Thoss (2007); Mavros et al. (2016); Kananenka et al. (2018).

II.1 Properties of the Morse potential

Proposed in 1929 Morse (1929) the Morse potential was one of the first empirical models to describe anharmonicities in the vibration of diatomic molecules, and is now widely used to describe complex anharmonic effects López-López et al. (2011); Kananenka et al. (2018); Hsieh and Cao (2018). The Morse potential presents various interesting features. First of all it is analytically solvable, yielding a finite number of bound states with negative energies given by

[TABLE]

These energies depend on the depth and the width of the potential through a new parameter $\Lambda$, such that $D_{k}=\frac{\omega_{k}\Lambda}{2}$ and $\alpha_{k}=\sqrt{\frac{\omega_{k}}{\Lambda}}$. Thus, $\omega_{k}$ is the frequency of the harmonic part of the potential, while the parameter $\Lambda$ tunes the anharmonic component. Note that while we consider the same anharmonicity parameter for all environment oscillators, their frequencies are linearly distributed.

Besides the bound states, the Morse potential present a continuous spectrum of positive energy eigenstates, also known as scattering states, that will not be considered here for simplicity. Overall, they will correspond to a continuous bath of free particles, to which the impurity may be coupled too and which may produce additional dephasing. In the new eigenstate basis, with the $n$-th energy eigenstate of the $k$-th oscillator written as $|kn\rangle$ and the respective eigenenergies $E_{kn}$, we find that Eqs. (2) and (4) can be written as

[TABLE]

where $c^{k}_{nm}=\langle kn|B_{k}|km\rangle$ are the coefficients for the bath part of the interaction operator.

A second interesting property is that the number of bound states is uniquely given by the integer part of $\Lambda+\frac{1}{2}$, when $\Lambda+\frac{1}{2}$ is not an integer. When $\Lambda+\frac{1}{2}$ is an integer, the number of bound states is $\Lambda-\frac{1}{2}$. This allows us to, by simply varying $\Lambda$, tune our bath between a spin bath limit, corresponding to $\Lambda\in]3/2,5/2]$, to a harmonic limit, recovered for $\Lambda\rightarrow\infty$. Note that at 3/2 the Morse oscillator has only one bound state. Between these limits, we also find two different types of regions in $\Lambda$ which will determine very strongly the nature of the dynamics:

1. I)

The regions where all environment states are strongly bounded, corresponding to $\Lambda\equiv\Lambda^{n,-\epsilon}=n+1/2-\epsilon$, with $n$ integer and $0\leq|\epsilon|\ll 1$ a real number. 2. II)

The regions where a new bound state is formed, which moreover is weakly bounded, corresponding to $\Lambda\equiv\Lambda^{n,\epsilon}=n+1/2+\epsilon$, with $n$ integer and $0<\epsilon\ll 1$ a positive real number.

A typical situation in both regions can be observed in Fig. 2 for the case of $n=2$, when choosing $\epsilon=0$ (region (I)) and $\epsilon=0.1$ (region (II)).

II.2 Properties of the open system dynamics

We now consider an initially decorrelated state so that the total density operator is written as $\rho(t=0)=\rho_{0}\otimes\rho_{E}^{\beta}$, where $\rho_{0}$ is the initial state of the system, and the environment is in a thermal equilibrium $\rho_{E}^{\beta}:=e^{-\beta H_{E}}/\textmd{Tr}_{E}\{e^{-\beta H_{E}}\}$, with the inverse temperature $\beta=1/k_{B}T$ and $k_{B}$ is the Boltzmann constant. In the general case, the reduced density matrix of the system in the interaction picture can be written as

[TABLE]

where ${\mathcal{U}}(t,t_{0})$ is the evolution operator in interaction picture, which can be expanded with a Dyson series as

[TABLE]

Considering this expansion as well as the case $H_{I}=SB$, we find that the time evolution of $\rho_{s}(t)$ can be re-written in terms of an infinite series of $l$-th order correlation functions which at zero times have the general form

[TABLE]

Here, the sum covers all configurations of $\{k_{1},\dots,k_{l}\}$, and can be decomposed into sums over sectors where different numbers of those $k$ coincide. As noted above, our model describes a linear system-environment coupling, such that for each oscillator $B_{k}\sim x_{k}$, with $x_{k}$ the position operator. Three important consequences can be extracted from such linear coupling to an anharmonic environment:

1. (i)

The first order term ($l=1$) in the expansion (11) needed to recover Eq. (8) is non-zero. This is because the coupling is linear and the environment Hamiltonian describes an asymmetric potential, such that $\left\langle B\right\rangle_{E}^{\beta}:=\textmd{tr}_{B}\{B\rho_{E}^{\beta}\}\neq 0$. This fact can be remedied by considering a renormalized version of the Hamiltonian, such that $H=\tilde{H}_{S}+H_{E}+\tilde{B}\sigma_{z}$, where $\tilde{H_{S}}=\omega_{s}\sigma_{z}+\sigma_{z}\left\langle B\right\rangle_{E}^{\beta}$, and $\tilde{B}=B-\left\langle B\right\rangle_{E}^{\beta}$. Hence, the correlations involved in Eq. (8) are equal to the ones in Eq. (11) but replacing each $B_{k}$ by their renormalized counterpart $\tilde{B}_{k}$. 2. (ii)

After the first order term, the most important contribution, i.e. the second order term or correlation function

[TABLE]

has a real part that does not decay to zero at finite temperatures but to a constant positive value. In terms of the eigen-values and eigenvectors of the environment $\{E_{kn},|kn\rangle\}$, the correlation function (12) can be formally written as $\alpha(t)=C(t)+C^{0}$, where

[TABLE]

corresponds to an offset, while

[TABLE]

corresponds to a time-dependent component, where we have defined $\Delta_{np}^{k}=E_{kn}-E_{kp}$. When we have a sufficiently large number of interfering phases, such time dependent factor decays. However, the time independent contribution will produce a real offset that will only be cancelled in the zero temperature limit, i.e. when we have that $C^{0}=\sum_{k}|\langle k0|\tilde{B}_{k}|k0\rangle|^{2}=|\langle\tilde{B}\rangle^{\infty}_{E}|^{2}=0$. 3. (iii)

The offset is particularly important in regions of type (II), i.e. for values of $\Lambda$ for which a weakly bound state exists, since these correspond to highly asymmetric bound states leading to a large overlap $\langle km|\tilde{B}_{k}|km\rangle$, where we define $m$ as the index of such weakly bound state. The presence of the weakly bounded state has two important consequences: First, it holds that:

[TABLE]

where we have just inserted the matrix element given in (133) and took the limit. In consequence, the offset becomes very large nearby such limit, and when the temperature is high enough so as to have a large initial population in the upper energy level. The second consequence is that all the transitions to the bound state with highest energy are suppressed. To show this, we consider a bound state $|l\rangle$ having lower energy than $|m\rangle$, then

[TABLE]

Here we took (B), inserted $N=m+\varepsilon$ and used that the functions containing $\varepsilon$ are continous functions so we could pull the limit into them. The consequence is that the time-dependent term of the correlation function, given by Eq. (14) is not affected by the presence of the bound state, contrary to the offset part. 4. (iv)

So far we have focused on the first and the second order moments of the expansion (11). However, in the present anharmonic case the structure of such expansion is much richer than in the harmonic case. First, the terms in $\textmd{Tr}_{E}\{\tilde{B}_{k_{1}}\cdots\tilde{B}_{k_{l}}\rho_{E}\}$ having an odd number of identical $k_{j}$ do no longer vanish. That is, $\textmd{Tr}_{E}\{\tilde{B}_{k}^{2n+1}\cdots\tilde{B}_{k}\rho_{E}\}\neq 0$, when $n$ is integer. Moreover, the even order components can no longer be decomposed as products of second order components of the form (12). Naturally, these two properties are fulfilled in the harmonic bath case, leading to system dynamical equations that can be written solely in terms of the correlation function (12). Indeed, Wick’s theorem holds in the harmonic case, since there $H_{k}$ is quadratic and one can decompose $B_{k}(t)=e^{iH_{k}t}B_{k}e^{-iH_{k}t}=f_{k}(t)b_{k}(0)+f^{*}_{k}(t)b^{\dagger}_{k}(0)$ (with $f_{k}(t)$ a time-dependent function). Away from this limit, such a decomposition is no longer possible and Wick’s theorem can no longer be applied.

One way to diminish the relative weight of higher order terms with respect to the second order one (12) is to choose an appropriate scaling for the coupling strengths $g_{k}$ present in the environment coupling operators $B_{k}$. In detail, if we choose them to scale as $g_{k}\sim 1/\sqrt{K}$, where $K$ is the number of oscillators, it can be shown that higher order correlations will eventually vanish in the limit of large $K$, an idea that was originally proposed in Makri (1999) and that is further discussed in Bramberger et al. .

In our case, we consider a fixed number of oscillators $K=40$ and explore to which extent the higher order terms are relevant to the description of the system. In order to be consistent with the harmonic oscillator limit, we further consider the standard choice in this limit for the frequency distribution and coupling strengths,

[TABLE]

Here $J(\omega)$ is the spectral density, which we consider to be of ohmic type

[TABLE]

with a special hard cut and $\omega_{c}$ is a cut-off frequency, and $\eta$ is a parameter with which we modulate the strength of the system-environment coupling. For the harmonic case, the chosen linear discretisation $\omega_{k}=k\Delta\omega$, with $\Delta\omega=2\omega_{c}/K$, gives rise to a revival time (which is the time at which finite size effects of the environment will start to occur) that in our case is $T=\pi/\Delta\omega\approx 20$.

II.3 The dynamical map

Because we are considering pure dephasing, we have that $\left[H_{S},S\right]=0$, which implies that the Hamiltonian is in block diagonal form and can be written as

[TABLE]

where we have defined the projectors $P_{\pm}:=|\pm\rangle\langle\pm|$ in terms of the eigenstates $|\pm\rangle$ of $\sigma_{z}$ with eigenvalues $\pm 1$, and we have decomposed both $H_{E}$ and $B$ in terms of local operators, $H_{k}^{\pm}:=H_{k}\pm B_{k}$ with $H_{k}$ corresponding to the $k$-th oscillator. Thus, the time-evolution operator can be computed as

[TABLE]

To obtain this expression, we have used the fact that the two terms in (22) commute, as well as the locality of the terms in the exponential to factorize it. Thus, computing Eq. (23) requires only having to exponentiate local operators which can be done very efficiently. Indeed, the total Hilbert space dimension scales like $d^{K}$, where $d$ is the local dimension, and we only have to exponentiate matrices of size $d$. Thus, the total density operator $\rho(t)=e^{-iHt}\rho_{0}\otimes\rho_{E}^{\beta}e^{iHt}$ can be rewritten as

[TABLE]

where we have defined $H^{\pm}:=\sum_{k}H_{k}^{\pm}$. Taking the partial trace over the bath degrees of freedom yields the reduced density matrix of the system, which has the following form

[TABLE]

where we have defined the decaying factor as

[TABLE]

where $\rho_{k}^{\beta}$ is the thermal state of the $k$-th oscillator with respect to its free Hamiltonian. We note that despite of the presence of an offset (13) which gives rise to an ill defined weak coupling master equation, we have numerically found that the map (28) is invertible for all the parameter regimes we considered here (not shown), which suggests that a time-local master equation is still well-defined Hall et al. (2014).

II.3.1 Harmonic and spin limits

In the harmonic oscillator limit, the decay factor acquires the usual form

[TABLE]

where

[TABLE]

where $\alpha(t-s)=\textmd{Tr}_{E}\{\tilde{B}(t)\tilde{B}(s)\rho_{E}^{\beta}\}$. Moreover, note that in the harmonic limit $\tilde{B}=B$, since $\langle B\rangle^{\beta}=0$. In addition, the offset (13) of the correlation function vanishes.

In the opposite limit of a spin bath, i.e. when $1.5<\Lambda\leq 2.5$, the operators $H_{k}$ and $B_{k}$ can be written in terms of the Pauli matrices $\sigma_{j}^{k}$ which are elements of the vector $\vec{\sigma}_{k}=(1\!\!1_{k},\sigma^{k}_{x},\sigma^{k}_{y},\sigma^{k}_{z})$, where we have defined the first element $\sigma_{0}^{k}=1\!\!1_{k}$ in terms of the unit operator in the Hilbert space of the $k$-th oscillator. In this representation, we find that

[TABLE]

where each component of $\vec{c}_{k}$ and $\vec{d}_{k}$ is defined as $c^{k}_{j}=\frac{1}{2}\textmd{Tr}_{k}(H^{+}_{k}\sigma^{k}_{j})$ and $d^{k}_{j}=\frac{1}{2}\textmd{Tr}_{k}(H^{-}_{k}\sigma^{k}_{j})$, respectively. Thus, the exponentials appearing in Eq. (28) can be simplified as

[TABLE]

III Dephasing time

We first analyze the dephasing time of the system. i.e. the decaying of the off-diagonal elements of the reduced density matrix. For all the numerical results we choose frequency units of $\omega_{s}=1$. Moreover, we consider the initial state

[TABLE]

Although this choice is rather arbitrary, it ensures that there are initial coherences that allow us to analyze dephasing. In other words our analysis is independent from the initial condition as long as it contains coherences. We define the decay time $\tau_{d}$ as the time such that

[TABLE]

where the reduced density matrix is given by Eq. (27).

Figs. 3(a,b) display the decay time $\tau_{c}$ for different values of the anharmonicity parameter $\Lambda$ and different initial temperatures of the bath. For the following discussion it is important to remember that as $\Lambda$ increases, the anharmonicity decreases. In detail, in Fig. 3(a), we explore the values of $\Lambda\equiv\Lambda^{n,\epsilon=0}=n+1/2$, which correspond to type (I) regions where the highest energy bound state is as far as possible from the dissociation limit, i.e. as strongly bounded as possible. With this choice, we find that $\tau_{c}$ has a quite smooth behavior with the anharmonicity, displaying a general tendency to increase towards the harmonic limit. Such trend reverses towards higher anharmonicities, where the number of bound states becomes increasingly limited and the dephasing slows down again.

As shown in Fig. 3(b), such a smooth behavior is disrupted when including a finer grid of $\Lambda$ values, now taking closer points spaced by $\delta\Lambda=0.1$. One can clearly identify that the regions of type (I) end up with a peak at values $\Lambda=n+1/2$ (where the top of the peaks correspond to the points shown in Fig. 3(a)), followed by a dip at values $\Lambda^{n,\epsilon}$ with $\epsilon$ very small, i.e. when the regions (II) begin. Indeed, as discussed in Fig. 3(a), at high anharmonicities each oscillator has only a few energy levels, which in principle hinders the dephasing. However, when entering each region (II) a new bound state is formed, and the offset becomes so large that it compensates for this effect and leads to a strong dephasing.

Another interesting feature of the figure is that in the spin region the decay time goes to infinity at low temperatures and at low $\varepsilon$ at the beginning of region (II). Indeed, we have seen in Eq. (18) that all the transitions to the bound state with the highest energy are suppressed in regions (II). Thus, when we only have two bound states the limit $\varepsilon\rightarrow 0$ gives rise to bath operators $H^{\pm}_{k}$ that are diagonal. To see this, we find that for $\Lambda=1.5+\varepsilon$ and small $\varepsilon$,

[TABLE]

and similarly $\lim_{\varepsilon\rightarrow 0^{+}}d^{k}_{x}=0$ and $\lim_{\varepsilon\rightarrow 0^{+}}d^{k}_{y}=0$. In addition,

[TABLE]

and their limit $\varepsilon\rightarrow 0^{+}$ will diverge as predicted by the general case eq. (16). If we now consider the eqs. (LABEL:eq:spindeph2) for a very small $\varepsilon$, and ignore for this argument the phase that comes from the $c_{0}$ and the $d_{0}$ contribution because it does not affect the absolute value of the coherences, we find

[TABLE]

where we have defined the phases \phi_{k}^{+}=\lim_{\varepsilon\rightarrow 0^{+}}\big{\|}\vec{c}_{k}\big{\|} and \phi^{-}_{k}=\lim_{\varepsilon\rightarrow 0^{+}}\big{\|}\vec{d}_{k}\big{\|}. Strictly speaking the above equations are ill-defined in the limit $\lim_{\epsilon\rightarrow 0}$. For this reason, we look at small but not infinitesimal $\epsilon$. Also, we shall note that the sign between the two terms in the second equation is flipped as compared to the first. This is because the $\textmd{tr}_{k}(H_{k}\sigma_{z})$ term in equation (42) becomes insignificant compared to the $\textmd{tr}_{k}(B_{k}\sigma_{z})$ term. The phases $\phi_{k}^{\pm}$ are going to infinity towards the limit $\varepsilon\rightarrow 0^{+}$. However, if we stay at a very small, yet not tiny value of $\varepsilon$ we can consider that $\phi_{k}^{\pm}\approx\Gamma^{\pm}$, where $\Gamma^{\pm}$ are two large phases, and define $\Gamma_{T}=\Gamma^{-}+\Gamma^{+}$, such that in Eq. (28) we now have

[TABLE]

where $\textmd{tr}_{k}(\sigma_{z}\rho_{k})=\tanh\left(\frac{\beta\Delta E_{k}}{2}\right)$, with $\Delta E_{k}$ the energy difference of the $k$-th spin. Thus at zero temperature we simply find

[TABLE]

The reason why the dephasing time does not increase to infinity at the beginning of other type (II) regions in the plot, is because as soon as there are more bound states some off diagonal transitions to intermediate states are allowed, which gives rise to dynamics even in the zero temperature case.

The large offset at the beginning of regions (II) can be observed in Fig. 3(c), and as argued before in point (iii) it is particularly relevant at high temperatures. From this plot, one may be tempted to believe that the offset is only important in such regions of type (II). However, when representing in Fig.(3)(d) one can see that relative value of the offset with respect to the maximal value of the time dependent part, i.e. $C^{0}/C(t=0)$ is non-negligible in all regions of $\Lambda$, including the region (I) where there is no weakly bounded state. Thus, on the one hand the offset explains the short dephasing time at high temperatures ($\beta=1$) with respect to lower temperatures, and on the other hand, for all temperatures it explains the tendency for a shorter dephasing towards $\Lambda$ smaller.

IV Outflow and backflow of information

We now analyze how the outflow and backflow of information is affected by the degree of anharmonicity of the environment. To this aim, we consider that the back-flow of information, as obtained from the BLP measure can be written as Cosco and Maniscalco (2018)

[TABLE]

where $[t_{1n},t_{2n}]$ are the time intervals over which $|\chi(t))|$ increases, and $\chi(t)$ is given by Eq. (28). In a similar way as in Cosco and Maniscalco (2018) we consider the ratio between the information backflow, given by Eq. (49) and the analogous quantity corresponding to the information flow to the environment, ${\mathcal{N}}_{+}$,

[TABLE]

The upper panel of Fig. (4)(a) represents the ratio $R$ for different values of anharmonicity $\Lambda$ and for the same coupling value ($\eta=2$) considered in the dephasing time analysis. We find that the only region where there is a non-negligible back-flow of information is the highly anharmonic one, particularly the spin-bath region. Fig. (4)(b) represents the same quantity but now for weaker coupling, $\eta=0.01$, showing that in this case the back-flow of information is also significant for higher values of $\Lambda$, particularly nearby the transition regions where a new environment bound state has been formed. Interestingly, the presence of back-flow is related to the following features:

• •

At high temperatures there is almost no back-flow. As seen before, the presence of a large offset gives rise to a fast dephasing time which in turn eliminates any possibility of back-flow and therefore non-Markovianity.

• •

At lower temperatures, the offset is not as important as to eliminate all the structure in the dephasing dynamics, and the maximal back-flow is observed precisely at the same values of $\Lambda$ where a new bound state has just been created (regions (II)).

The dynamics of $\chi(t)$ can be observed in Fig. (4)(c) (strong coupling) and Fig. (4)(d) (weak coupling) for $\beta=7$. In detail, one can observe that while the real and imaginary parts of $\chi(t)$ are oscillatory, its absolute value is a monotonically decreasing function for $\Lambda=n+0.5=2.5$ and becomes non-monotonically decreasing (i.e. giving rise to back-flow) for $\Lambda=2.6$, i.e. for a Morse potential value having a weak bound state.

An important comment is here in order. Since our environment is finite, a border effect is to be expected at a certain time. For $K=40$ oscillators, the revival time can be estimated to be $T=20$ in the harmonic limit (see de Vega et al. (2015) and references therein for details), but it is much harder to compute in the anharmonic case. The reason is that, while in the harmonic case the revival time is the time at which the correlation function, having decayed, starts to grow again, in the anharmonic case we have seen that such function is not the only one that comes into play. For $\eta=2$ the system off-diagonal elements have already decayed to zero at all considered values, which allows to consider our analysis of dephasing time. However, this is not the case for smaller values, like $\eta=0.01$, where at $T=20$ there are curves that have not yet decayed. This means that the present analysis of forward and back-flow has been performed within a time-frame up to $T=20$ where neither we can consider that the system has always completely decayed, nor we can exclude border effects in the anharmonic limit. A rigorous analysis of such finite size effects produced by non-Harmonic environments is out of the scope here, but would be very interesting. In particular, it is interesting to ask whether the measured back-flow is due to such finite-size effects, or to the nature of the environment itself, as considered in the harmonic case.

Finally, we computed the non-Markovianity by considering the Gaussian map $\chi^{\textmd{Gauss}}(t)$ given by Eq. (29) and the correlation function (12), finding that is zero for all parameter regimes considered. It remains to be further analyzed if Gaussian evolutions (or their related properties such as non-Markovianity) obey some extremality property like that of Gaussian states Wolf et al. (2006).

V Non-Gaussian nature of the bath

We have seen in the last sections that nearby the limit where a weakly bounded upper state is present in the environment oscillators there is a strong information back-flow, particularly at low temperatures. In this section we will further explore the nature of such non-Markovianity and show that it is linked to a strong non-Gaussianity of the bath. To this aim, we will compare the map (27) computed with the exact non-Gaussian $\chi(t)$ corresponding to the Morse environment, Eq. (28), with the Gaussian $\chi^{\textmd{Gauss}}(t)$ (again as given by Eq. (29) with the correlation function (12)). Doing so, we are comparing the exact dynamics with that obtained by assuming that the map is Gaussian, and therefore fully determined by the second order moment.

The density plot (5) shows the error $E(t)=|\chi(t)-\chi^{\textmd{Gauss}}(t)|$ for low (left panel) and high (right panel) temperatures, for different times and anharmonicities. It can be observed in general that the Gaussian map becomes a particularly bad approximation in the regions nearby the formation of a new bound state (but not only) and at low temperatures. This behavior can be best observed in Fig. (6), which represents the time-averaged error

[TABLE]

where $\rho_{s}(t)$ represents the exact evolution, $\rho_{s}^{\textmd{Gauss}}(t)$ the Gaussian-approximated evolution as given by $\chi^{\textmd{Gauss}}(t)$, and $D(A,B)=(1/2)\textmd{Tr}\{\sqrt{(A-B)^{2}}\}$ represents the trace distance between the two density matrices $A$ and $B$. It can be observed that the error is large for all temperatures explored except for very high values ($\beta=1$). We note that in the scale that we show it seems to increase slightly with $\Lambda$. This is a local effect which disappears, as expected, when going further to the harmonic limit.

The lower panels show how the emergence of a new bound state (for $\Lambda=2.6$) gives rise to a strong departure of the Gaussian approximation with respect to the exact case.

VI Conclusions

We have analyzed the dephasing dynamics of an impurity coupled to an environment of oscillators containing a varying degree of anharmonicity. We first analyze the variation of the dephasing time with the anharmonic parameter of the environment, $\Lambda$, and find the following:

• •

On a coarse-grained scale in $\Lambda$, the decoherence time increases towards both the harmonic and anharmonic limits. In the anharmonic limit the decoherence is slow due to the few bound states in the environment. When the anharmonicity grows this effect is compensated by the fact that the environment correlation function does not decay to zero, but to a real value, the offset. Such offset occurs due to the presence of bound states $|p\rangle$ that are highly asymmetric, i.e. for which $\langle p|x_{k}|p\rangle>0$ where $x_{k}$ is the displacement operator with respect to the potential minimum. Hence, the offset becomes negligible towards the harmonic limit, which explains the slowing down of the dephasing.

• •

On a finer scale, for a given integer number $n$ there are two different regions in $\Lambda$: type (I) regions where $\Lambda=n+1/2-\epsilon$, $\epsilon$ small, where we have an integer number of bound states $N=n$ that are all strongly bounded, i.e energetically far from the dissociation limit, and type (II) regions for $\Lambda=n+1/2+\epsilon$, $\epsilon$ small, where we have $N=n+1$ bound states and the highest energy one, $|m\rangle$, is weakly bounded. The displacement of such weakly bounded state with respect to the center of the potential is large, leading to an overlap $\langle n|x_{k}|n\rangle\gg 1$ that produces a particularly large offset and thus a very strong dephasing in those regions.

Secondly, we analyze the non-Markovianity or back-flow of information, and find that it only occurs at high anharmonicities and that it is concentrated in regions (II). Moreover, the backflow is only relevant at low temperatures, which suggests that while it is indeed linked to a large anharmonicity it is nevertheless hindered by the presence of the offset. Indeed, at high temperatures such offset becomes dominant and gives rise to a dephasing dynamics that is too fast to allow for backflow.

A similar situation occurs when considering how relevant or necessary it is to account for the non-Gaussianity of the bath. To check this we consider the correlation function of the bath and build the corresponding Gaussian map to compare it with the exact non-Gaussian one. As it turns out, for high temperatures the effect of the non-Gaussianity becomes negligible and both exact and Gaussian versions give approximately the same dynamics. Indeed, the offset is also present in the Gaussian version of the map and is so large that it dominates the decay. The non-Gaussianity, and therefore the effect of the higher order moments, starts to be important for lower temperatures, when the offset is not so relevant but still a large anharmonicity is present. In addition, the Gaussian analogue map predicts no-backflow in any regime, which reinforces the idea that the backflow is related to the anharmonicity and non-gaussian character of the bath.

Overall, this analysis unveils the presence of a variety of dynamical regimes for an impurity coupled to an anharmonic environment, and which depend not only on the number of bound states in the environment, $|m\rangle$, but also on their location within the anharmonic potential and therefore on their degree of asymmetry. While the relevance of each of these regimes and situations remains to be determined and analyzed in each physical situation and specific model, this study is perhaps one of the first to systematically discuss how rich the dynamics of an open system may become when smoothly departing beyond the standard harmonic bath situation. In this context, concepts such as environment size effects, origin and nature of the information backflow Breuer et al. (2009) and other related concepts such as divisibility A et al. (2009), existence of a dynamical equation (i.e. invertibility of the map) Li et al. (2018); Hall et al. (2014) or validity of the weak coupling approximation may need to be revisited and further analyzed in the future.

Acknowledgements.

The authors acknowledge D. Wierichs, C. Parra, A. Recati, C. Hubig, Christiane Koch, N.O. Linden and U. Schollwöck for interesting discussions. This research was financially supported by DFG-grant GZ: VE 993/1-1 and supported in part by the National Science Foundation under Grant No. NSF PHY-1748958.

Appendix A Bound states of Morse oscillators

We discuss in this section how to derive the bound states of the Morse potential. To this aim we will follow the derivations in Morse (1929) but describe in more detail some of the steps that in the original paper were only sketched. In addition we will put more emphasis on the criteria to decide whether a solution is physical or not. We start with the Hamiltonian:

[TABLE]

We will now do a few transformations in order to solve the stationary Schrödinger equation related to this Hamiltonian. First we define $x^{\prime}:=\alpha x$. Leading to:

[TABLE]

Now we rescale this to make the Hamiltonian dimensionless: $H_{0}=\frac{2m}{\hbar^{2}\alpha^{2}}H^{\prime}$. Then the Hamiltonian becomes:

[TABLE]

With $D^{\prime}=\frac{2m}{\hbar^{2}\alpha^{2}}D:=(N+\frac{1}{2})^{2}$. From now on we will omit the primes. The bound states $\psi(x)$ for the energy $E<0$ must fulfill the following:

[TABLE]

In order to make this a linear differential equation we substitute $z=(2N+1)\exp(-x)$. Which leads us to:

[TABLE]

Now we substitute $\psi(z)=z^{b/2}\exp(-az)F(z)$. After a rather long but straight forward calculation this results in:

[TABLE]

In order to simplify this equation we choose $a=1/2$ and $b^{2}=-4E$. Because with this choice we obtain:

[TABLE]

This is the Laguerre equation:

[TABLE]

As discussed in Morse (1929), one set of solutions for this equation are the Laguerre polynomials and indeed these are the only polynomial solutions. In the original paper other possible solutions are immediately discarded as unphysical. In the following we will give proof that this is the case, which justifies the approach in Morse (1929). Hence we start with the series ansatz:

[TABLE]

Inserting this into the Laguerre equation gives us:

[TABLE]

Shifting the indices in the first two sums:

[TABLE]

Thus for $n=0$ we obtain the following relation between the first two coefficients:

[TABLE]

For $n>0$ we then obtain the recurrence relation:

[TABLE]

Iterating this equation gives us:

[TABLE]

Where $(x)_{n}$ is the Pochhammer symbol defined as

[TABLE]

which is equivalent to: $(x)_{n}=\prod_{l=0}^{n-1}(x+l)$. Please note that because of the definition of the Pochhammer symbols the series terminates if $\lambda$ is a non-negative integer. This leaves us with two cases: One is that $\lambda$ is a non-negative integer and the other one is when this is not the case. Lets first investigate the case where $\lambda$ is not an non-negative integer. We now show that in this case we do not obtain physical solutions of the Schrödinger equation. The argument is that if the series does not terminate, then the wavefunction will not be normalisable, and will therefore not represent a physical bound state. One can see this by splitting the series into a polynomial and a non-polynomial part and then approximate the behaviour of the non-polynomial part. The solution of the Laguerre equation is:

[TABLE]

First of all we note that with out loss of generality we can say that there exists an $N$ such that for all $n>N:(-\lambda+N+n)>0$. Now we split the sum as

[TABLE]

We can restrict ourselves to $b>0$, since if $b<0$ then we would have a singularity in the wave function due to the substitution we made earlier. This would make it unphysical. With this the coefficients in the second sum all have the same sign. The first sum is just a polynomial and can not give us any problems, because in our substitution from earlier we have a term that scales as $\exp(-z/2)$. Thats why we abandon the polynomial part in this consideration. Now if $-\lambda>(1+b)$ then we are already done. Because then the following holds:

[TABLE]

The first equality holds because we are only interested in the solution for $z$ in $[0,\infty)$. This means that in the case where $-\lambda>(1+b)$ the solution scales even stronger in $z$ then $\exp(z)$ and thus it is unphysical, because neither our factor of $\exp(-z/2)$ nor any polynomial can compete with that. So now lets see how this goes for $-\lambda<(1+b)$. In this case this is a bit more difficult. First we need to note that $\lim_{n\rightarrow\infty}(-\lambda+n)/(1+b+n)=1$. This means for all $\epsilon>0$ there exists an $N^{\prime}$ such that for all $n\geq N^{\prime}$:

[TABLE]

From which follows that for $n\geq N^{\prime}:\frac{-\lambda+n}{1+b+n}>1-\epsilon$. With this we can investigate the non-polynomial part further:

[TABLE]

Now lets just look at the non-polynomial part again:

[TABLE]

Here $c_{N^{\prime}}:=\frac{(-\lambda)_{N^{\prime}}}{(1+b)_{N^{\prime}}}$ and $\epsilon$ can be chosen to be any positive real number. Thus if one chooses for example $\epsilon=1/4$ one can see that the wave function is not normalisable. Again the second term does not make it normalisable again, since it is just a polynomial. With this extensive proof we showed that the only physical solutions are obtained for $\lambda$ being a non-negative integer. With this we obtain an equation for $b$ and for the energy since $b^{2}/4=-E$:

[TABLE]

Here $\lambda\rightarrow n$ is a non-negative integer. Remember that $b$ is only allowed to be positive, because otherwise either the wave functions or the probability density will have a singularity at $z=0$. This means that if $N$ is an natural number it is the number of bound states, if not then the integer part of $N+1$ is the number of bound states. The wave function is then given by:

[TABLE]

With proper normalisation this becomes:

[TABLE]

Here we have defined $z$ and $N_{n}$ as follows:

[TABLE]

Note that all these equations are still for a dimensionless Hamiltonian. When we reintroduce the parameters we rescaled the Hamiltonian with in the beginning we get:

[TABLE]

Where $\Lambda$ is defined as $\Lambda:=N+\frac{1}{2}$. Now for a given Morse Oscillator we want the part that is proportional to $n+\frac{1}{2}$ to be the harmonic part thus it should be $\hbar\omega(n+\frac{1}{2})$. So now we change the parameters of the Morse Oscillator from $\alpha$ and $D$ to $\omega$ and $\Lambda$:

[TABLE]

As mentioned earlier we obtain a set $\{\omega_{k},g_{k}\}$ by discretizing the spectral density. The $\omega$ in the above equations corresponds to the $\omega_{k}$ that we obtain by the discretization. For the implementation of the Morse environment we choose $\Lambda$ to be constant for all oscillators, thus they all have the same number of bound states.

Appendix B Matrix elements of the position operator

In our description, it is highly convenient to write the Morse Hamiltonian in its diagonal form. This in turn means that the interaction Hamiltonian should also be expressed in such basis. To this aim, it is necessary to write the matrix element of the position operator in such basis, a computation that as we will see in the following is not at all trivial. To proceed with it, we shall make use of the wave function of the bound states, as computed in the previous appendix. In the following, we consider the same strategy as in Ref. de Lima and Hornos (2005). That is we compute them as:

[TABLE]

Here $x$ is still supposed to be dimensionless, thus the rescaling of the last section is still in place. Actually, computing $\langle n|\exp(i\eta x)|m\rangle$ is pretty straight forward apart from the fact that one has to use the gaussian hypergeometric theorem and the Saalschütz theorem. Those will both be stated when we use them. Nevertheless, the most delicate part of the computation is taking the limit after the derivative. Since this important step was not detailed in de Lima and Hornos (2005), we will present it here in more detail.

Let us first consider that $\exp(i\eta x)=(2N+1)^{i\eta}z^{-i\eta}$, such that

[TABLE]

Then, considering the wave functions (85) we can now compute:

[TABLE]

where we have defined

[TABLE]

Using the definition of the $\Gamma$-function: $\Gamma(z)=\int_{0}^{\infty}dte^{-t}t^{z-1}$ one obtains:

[TABLE]

where we have defined a function that contains all terms that depend on ${m^{\prime}}$

[TABLE]

We now we consider that

[TABLE]

in order to re-express Eq. (99) in terms of Pochhammer symbols (69),

[TABLE]

This allows us to perform the sum over ${m^{\prime}}$ by using the gaussian hypergeometric theorem, which states (see de Lima and Hornos (2005) section 3.1)

[TABLE]

Note that $(-m)_{m^{\prime}}$ becomes 0, when ${m^{\prime}}>m$, which allows us to apply the theorem even though it is meant for summation up to infinity. When using this we obtain:

[TABLE]

Thus the matrix element becomes:

[TABLE]

Now we want to organise this expression into two terms, where one contains everything depending on $n^{\prime}$. Also we want to have the things that depend on $n^{\prime}$ to be formulated in terms of Pochhammer symbols again, so that we can ultimatly use Saalschütz theorem. First we have to note the following:

[TABLE]

With this we obtain:

[TABLE]

Where $A_{mn}(\eta)$ and $\Gamma_{mn}(\eta)$ are defined as follows:

[TABLE]

Where ${}_{3}F_{2}$ is a generalised hypergeometrical function as can also be seen in de Lima and Hornos (2005) section 3.1 eq. 18 and Bateman (1953) chapter 4. The derivative can now be computed as:

[TABLE]

Up to now we went pretty much along the lines of de Lima and Hornos (2005) apart from expressing things in a much more detailled manner of course. Taking this derivative and the limit afterwards is the most problematic part of this computation and since this step is skipped in de Lima and Hornos (2005) we are going to do this in all detail here.

Let us distinguish two cases. The offdiagonal elements namely $m>n$ and the diagonal elements. Note that the condition $m>n$ does not imply a loss of generality since x is hermitian. First we take a look at the offdiagonal elements. For $m>n$ one can see that in the denominator of $A_{mn}(\eta)$ there is a term $\Gamma(n-m+1+i\eta)$ which approaches a singularity in the limit $\eta\rightarrow 0$. Thus $A_{mn}(\eta)$ becomes 0 in the limit. If $\frac{d}{d\eta}\Gamma_{nm}(\eta)$ does not have a diverging term in the limit this means, that we only have to take into account the term containing the derivative of $A_{mn}(\eta)$. By looking at the form of $\Gamma_{mn}(\eta)$ we can already see that taking a derivative and sending $\eta$ to 0 will not give us something divergent. Thus we only care for the derivative of $A_{mn}(\eta)$.

[TABLE]

where $\psi(x):=\frac{d}{dx}\ln(\Gamma(x))=\frac{\frac{d}{dx}\Gamma(x)}{\Gamma(x)}$ is the digamma function. Here we skipped a few steps, but they are just taking the derivative and realising, that in the limit $\eta\rightarrow 0$ only this term survives, because in all the other terms there is no singularity in the enumerator. In order to obtain this limit one can do a nice trick. Namely one can use the following two identities(Bateman (1953) chapter 1.2 eq. 6 and chapter 1.7 eq. 11):

[TABLE]

With the use of these equations, the limit becomes:

[TABLE]

With this we obtain(for $m>n$):

[TABLE]

Thus the only thing that remains to be done in the offdiagonal case is to take the limit for $\Gamma_{mn}(\eta)$. If one takes this limit one obtains:

[TABLE]

A generalized hypergeometric function ${}_{3}F_{2}$ is Saalschützian if its parameters are in the following form and the Saalschütz theorem states that (see Bateman (1953) chapter 4.4):

[TABLE]

In the limit $\eta\rightarrow 0$ the parameters for our generalised hypergeometric function are of exactly that form, thus:

[TABLE]

Here, we have considered that

[TABLE]

With this we obtain:

[TABLE]

Now we can finally obtain the matrix elements for the offdiagonal term with $m>n$:

[TABLE]

The last two lines show our result for the offdiagonal elements. Now for the diagonal elements $A_{mn}(\eta)$ does not have a singularity in the denominator, thus both terms contribute. Let us first compute the contribution from $A_{nn}(\eta)$:

[TABLE]

Here we have not detailed the derivation, which basically consists on taking the derivative of each term of $A_{nn}(\eta)$ in Eq. (108) and then taking the limit. Please note that the derivation of $\lim_{\eta\rightarrow 0}\Gamma_{nm}(\eta)$ was not relying on $m>n$ so it also holds here, and thus $\lim_{\eta\rightarrow 0}\Gamma_{nn}(\eta)=1$. So now comes the tricky part. In order to obtain a nice form for the derivative of $\Gamma_{nn}(\eta)$ we would like to use the Saalschütz theorem of Eq. (LABEL:Saa_th). But if we are not in the limit where $\eta\rightarrow 0$ our ${}_{3}F_{2}$ is not Saalschützian. This is where we use a trick. We define another $\tilde{\Gamma}_{nn}(\eta)$, that is Saalschützian and that has the same derivative in the limit $\eta\rightarrow 0$. Define $\tilde{\Gamma}_{nn}(\eta)$ as follows:

[TABLE]

First we have to check that in the limit of $\eta\rightarrow 0$ this indeed has the same derivative as $\Gamma_{nn}(\eta)$. So lets take a look at how this is expressed in terms of Pochhammer symbols:

[TABLE]

Because of the term $(-i\eta)_{n^{\prime}}$ this is proportional to $\eta$, so the only term that contributes to the derivative in the limit of $\eta\rightarrow 0$ is the one where one takes the derivative of $(-i\eta)_{n^{\prime}}$ (in other words, the derivative of the other terms with $\eta$ will always be multiplied by $(-i\eta)_{n^{\prime}}$ and therefore vanish in the limit $\eta\rightarrow 0$). This leads to:

[TABLE]

Now let us compare this to $\Gamma_{nn}(\eta)$:

[TABLE]

With the same argument as for $\tilde{\Gamma}_{nn}(\eta)$ we obtain:

[TABLE]

Now we use the Saalschütz theorem to compute $\tilde{\Gamma}_{nn}(\eta)$:

[TABLE]

Thus:

[TABLE]

Also note that $\lim_{\eta\rightarrow 0}A_{nn}(\eta)=1$. With this we can compute the diagonal matrix elements:

[TABLE]

For the sake of completeness let us state the offdiagonal elements for $m>n$ here again:

[TABLE]

Here we still refer to the dimensionless $x$ thus in order to obtain the matrix elements of $b^{\dagger}+b$ one only has to rescale the ones from (133) and (B) by a factor $\sqrt{2/\Lambda}$, namely

[TABLE]