Zero-dimensional models for gravitational and scalar QED decoherence

Starting with the Lagrangian $L_{\mathrm{grav}}$ (11), the system and bath momentum coordinates are

	$\displaystyle p$	$\displaystyle=\frac{\partial L}{\partial\dot{x}}=M\dot{x}\left(1-\sum_{i}\lambda_{i}q_{i}\right),$		(13)
	$\displaystyle p_{i}$	$\displaystyle=\frac{\partial L}{\partial\dot{q}_{i}}=m\dot{q}_{i},$		(14)

where we omit the ‘grav’ subscript from now on. The model Hamiltonian is

	$\displaystyle H=$	$\displaystyle\frac{p^{2}}{2M}\left(1-\sum_{i}\lambda_{i}q_{i}\right)^{-1}+\frac{1}{2}M\Omega^{2}x^{2}\left(1+\sum_{i}\lambda_{i}q_{i}\right)$
		$\displaystyle+\sum_{i}\left(\frac{p_{i}^{2}}{2m}+\frac{1}{2}m\omega_{i}^{2}q_{i}^{2}\right).$		(15)

For weak system-environment (bath) coupling as in the case of weak field gravity, we can Taylor expand the kinetic energy coupled bath term to obtain approximately

	$\displaystyle H=$	$\displaystyle\left(\frac{p^{2}}{2M}+\frac{1}{2}M\Omega^{2}x^{2}\right)\left(1+\sum_{i}\lambda_{i}q_{i}\right)$
		$\displaystyle+\sum_{i}\left(\frac{p_{i}^{2}}{2m}+\frac{1}{2}m\omega_{i}^{2}q_{i}^{2}\right).$		(16)

Quantizing and expressing the Hamiltonian (16) in terms of the oscillator system and bath creation and annihilation operators, which are defined through the respective relations $x=\sqrt{\frac{\hbar}{2M\Omega}}(a+a^{{\dagger}})$, $p=i\sqrt{\frac{M\Omega\hbar}{2}}(a^{\dagger}-a)$, $q_{i}=\sqrt{\frac{\hbar}{2m\omega_{i}}}(a_{i}+a_{i}^{{\dagger}})$, $p_{i}=i\sqrt{\frac{m\omega_{i}\hbar}{2}}(a_{i}^{\dagger}-a_{i})$, the scalar gravity model Hamiltonian is

	$\displaystyle H=$	$\displaystyle\hbar\Omega\left(a^{\dagger}a+\frac{1}{2}\right)\left(1+\sum_{i}\lambda_{i}\left(a_{i}+a_{i}^{\dagger}\right)\right)$
	$\displaystyle+$	$\displaystyle\sum_{i}\hbar\omega_{i}\left(a_{i}^{\dagger}a_{i}+\frac{1}{2}\right),$		(17)

where the system-bath coupling is redefined as $\lambda_{i}=\sqrt{\frac{\hbar}{2m\omega_{i}}}\lambda_{i}$. We recognize in Eq. (17) the familiar form of the standard optomechanical Hamiltonian, but with a bath of mechanical oscillator modes (labelled by index $i$) in contrast to the usually considered situation of just one mechanical mode Aspelmeyer et al. (2014).

Solving for the quantum evolution, we will make the common assumption that the system and bath are in an initial product state $\rho_{s}\otimes\rho_{\mathrm{bath}}$. While the latter assumption facilitates solving for the quantum dynamics, it is not always justified experimentally, since it necessarily requires that the system quantum state can be sufficiently isolated and prepared quickly enough compared to the interaction time scale with the bath degrees of freedom. While such an approximation may be justified for an electromagnetic environment under certain conditions, a mass-energy system can never be isolated from its gravitational environment. Nevertheless, as we shall see, the ability to solve exactly for the scalar gravity model quantum dynamics will give insights into the consequences of assuming a product state.

It is convenient to work in a basis of system number states and bath coherent states $|n,\{\alpha_{i}\}\rangle$; the time evolution for such a state can be written as:

	$\displaystyle e^{-\frac{iHt}{\hbar}}|n,\{\alpha_{i}\}\rangle$	$\displaystyle=\exp\Bigg{(}-\frac{it}{\hbar}\Big{[}\hbar\Omega\left(n+\frac{1}{2}\right)\big{(}1+\sum_{i}\lambda_{i}(a_{i}$
	$\displaystyle+a_{i}^{\dagger})\big{)}+\sum_{i}\hbar$	$\displaystyle\omega_{i}\left(a_{i}^{\dagger}a_{i}+\frac{1}{2}\right)\Big{]}\Bigg{)}\left|n,\{\alpha_{i}\}\right\rangle.$		(18)

Following the analysis of Ref. Bose et al. (1997), Eq. (18) can be evaluated as:

	$\displaystyle e^{-\frac{iHt}{\hbar}}|n,\{\alpha_{i}\}\rangle=$
	$\displaystyle\exp\Bigg{(}-it\left[\Omega\left(n+\frac{1}{2}\right)+\sum_{i}\left(\frac{\omega_{i}}{2}-\frac{\Omega^{2}\lambda_{i}^{2}\left(n+\frac{1}{2}\right)^{2}}{\omega_{i}}\right)\right]$
	$\displaystyle-\sum_{i}\frac{i\left(n+\frac{1}{2}\right)^{2}\lambda_{i}^{2}\Omega^{2}}{\omega_{i}^{2}}\sin\omega_{i}t+\frac{1}{2}\sum_{i}\frac{\lambda_{i}}{\omega_{i}}\left(n+\frac{1}{2}\right)\Omega$
	$\displaystyle\times\left[\alpha_{i}^{*}\left(1-e^{i\omega_{i}t}\right)-\alpha_{i}\left(1-e^{-i\omega_{i}t}\right)\right]\Bigg{)}$
	$\displaystyle\times\left|n,\left\{\alpha_{i}e^{-i\omega_{i}t}-\frac{\Omega\left(n+\frac{1}{2}\right)}{\omega_{i}}\lambda_{i}\left(1-e^{-i\omega_{i}t}\right)\right\}\right\rangle.$		(19)

Supposing the bath to initially be in a thermal state, we can express its initial density matrix in the coherent basis as follows Bose et al. (1999):

	$\displaystyle\rho_{\mathrm{bath}}=$	$\displaystyle\prod_{i}\frac{1}{\pi\left(e^{\beta\hbar\omega_{i}}-1\right)}\int d\alpha_{i}^{2}\exp\Big{(}-|\alpha_{i}|^{2}$
	$\displaystyle\times$	$\displaystyle\left(e^{\beta\hbar\omega_{i}}-1\right)\Big{)}|\alpha_{i}\rangle\langle\alpha_{i}|,$		(20)

where $\beta^{-1}=k_{B}T$, with $k_{B}$ Boltzmann’s constant and $T$ the bath temperature. Decomposing the initial system-environment state in the number state basis:

$\displaystyle\rho_{\mathrm{initial}}=\sum_{n,n^{\prime}}C_{nn^{\prime}}|n\rangle\langle n^{\prime}|\otimes\rho_{\mathrm{bath}},$

(21)

we have for the time evolution of the number state outer products after tracing out the bath degrees of freedom:

		$\displaystyle|n(t)\rangle\langle n^{\prime}(t)|=|n\rangle\langle n^{\prime}|\exp\Bigg{(}-it\Big{[}\Omega(n-n^{\prime})+\sum_{i}\frac{\Omega^{2}\lambda_{i}^{2}}{\omega_{i}}$
	$\displaystyle\times$	$\displaystyle(n^{\prime}+n+1)(n^{\prime}-n)\Big{]}+i\sum_{i}\frac{\lambda_{i}^{2}\Omega^{2}}{\omega_{i}^{2}}\sin(\omega_{i}t)(n^{\prime}+n+1)$
	$\displaystyle\times$	$\displaystyle(n^{\prime}-n)-2\sum_{i}\left(\frac{\Omega\lambda_{i}(n-n^{\prime})}{\omega_{i}}\right)^{2}\coth\left(\frac{\beta\hbar\omega_{i}}{2}\right)$
	$\displaystyle\times$	$\displaystyle\sin^{2}\left(\frac{\omega_{i}t}{2}\right)\Bigg{)}.$		(22)

As we shall show later below in Sec. III.2, the time evolution of an arbitrary initial reduced oscillator system state can be obtained by decomposing in terms of the number state outer product solutions (22).

In order to carry out the sum over bath degrees of freedom in Eq. (22), we will assume an “Ohmic” bath spectral density with exponential cut-off frequency $\omega_{c}$:

$\displaystyle\pi\sum_{i}\lambda_{i}^{2}\delta\left(\omega-\omega_{i}\right)=C\omega e^{-{\omega}/{\omega_{c}}}.$

(23)

We note that the bath spectral density is in general determined by the microscopic coupling strength $\lambda_{i}$ and the spectral structure of the bath modes; here we consider the linear in $\omega$ Ohmic dependence since a thermal graviton bath is Ohmic Blencowe (2013); a cut-off is necessary in order to tame infinities that would otherwise arise for $\omega\rightarrow\infty$ in the intermediate stages of the analysis; the exponential cut-off form is primarily motivated by calculational convenience, enabling exact solutions to be obtained for the oscillator system reduced dynamics in the number state basis. In the gravitational decoherence settings, the cut-off energy scale is related to the nature of the ‘Planckian’ physics, which is not known yet. From the perspective of the effective field theory approach to quantum gravity, this unknown high energy physics is assumed to decouple from the low energy dynamics, appearing only in renormalized parameters that are determined phenomenologically Donoghue (1994); as we shall see below, the cut-off dependence can be absorbed in part through a frequency and non-linear Kerr-type self-interaction renormalization. However, the cut-off does affect the initial decoherence dynamics, which is tied to the artificial initial product state assumption.

Using Eq. (23) to replace the sum in Eq. (22) with an integral over the continuous variable $\omega$, we obtain

	$\displaystyle|n(t)\rangle\langle n^{\prime}(t)|=|n\rangle\langle n^{\prime}|\exp\Bigg{(}-i\Omega(n-n^{\prime})t+i\frac{C\Omega^{2}}{\pi}$
	$\displaystyle\times(n-n^{\prime})(n+n^{\prime}+1)\left[\omega_{c}t-\tan^{-1}(\omega_{c}t)\right]-\frac{2C\Omega^{2}}{\pi}$
	$\displaystyle\times(n-n^{\prime})^{2}\int_{0}^{\infty}d\omega\omega\coth\left(\frac{\beta\hbar\omega}{2}\right)\frac{\sin^{2}(\frac{\omega t}{2})}{\omega^{2}}e^{-{\omega}/{\omega_{c}}}\Bigg{)}.$		(24)

Note that for an optomechanical system realization, the environment would have a finite extent resulting in a non-zero, lower frequency cut-off $\omega_{1}$: $0<\omega_{1}\ll\omega_{c}$; here we set the lower frequency cut-off $\omega_{1}=0$, reflecting the actual gravitational environment with effectively infinite spatial extent. In the case of a non-zero lower frequency cut-off, the system will not fully decohere in the long time limit Xu and Blencowe (2021), while for the $\omega_{1}=0$ gravitational setting, the dephasing dynamics displays exponential decay with time as we shall see in the following. The integral in the above expression can be evaluated analytically to give

		$\displaystyle\int_{0}^{\infty}d\omega\omega\coth\left(\frac{\beta\hbar\omega}{2}\right)\frac{\sin^{2}(\frac{\omega t}{2})}{\omega^{2}}e^{-{\omega}/{\omega_{c}}}$
	$\displaystyle=$	$\displaystyle\frac{1}{4}\ln\left(1+t^{2}\omega_{c}^{2}\right)$
	$\displaystyle+$	$\displaystyle\frac{1}{2}\ln\left[\frac{\Gamma^{2}\left(\frac{1}{\beta\hbar\omega_{c}}+1\right)}{\Gamma\left(\frac{1-it\omega_{c}}{\beta\hbar\omega_{c}}+1\right)\Gamma\left(\frac{1+it\omega_{c}}{\beta\hbar\omega_{c}}+1\right)}\right].$		(25)

Taking the limit ${\beta\hbar\omega_{c}}\rightarrow\infty$ (i.e, upper cut-off frequency large compared to the bath temperature), we have

$\displaystyle\frac{\Gamma^{2}\left(\frac{1}{\beta\hbar\omega_{c}}+1\right)}{\Gamma\left(\frac{1-it\omega_{c}}{\beta\hbar\omega_{c}}+1\right)\Gamma\left(\frac{1+it\omega_{c}}{\beta\hbar\omega_{c}}+1\right)}\rightarrow\frac{\beta\hbar}{\pi t}\sinh\left(\frac{\pi t}{\beta\hbar}\right).$

(26)

With approximation (26), Eq. (24) becomes

	$\displaystyle|n(t)\rangle\langle n^{\prime}(t)|=|n\rangle\langle n^{\prime}|\exp\bigg{[}-i\Omega(n-n^{\prime})t+i\frac{C\Omega^{2}}{\pi}$
	$\displaystyle\times(n-n^{\prime})(n+n^{\prime}+1)\left[\omega_{c}t-\tan^{-1}(\omega_{c}t)\right]-\frac{C\Omega^{2}}{\pi}$
	$\displaystyle\times(n-n^{\prime})^{2}\left(\frac{1}{2}\ln\left(1+t^{2}\omega_{c}^{2}\right)+\ln\left[\frac{\beta\hbar}{\pi t}\sinh\left(\frac{\pi t}{\beta\hbar}\right)\right]\right)\bigg{]}.$		(27)

We now discuss the various terms appearing in Eq. (27). First, note that the outer product is time-independent for $n=n^{\prime}$, a consequence of the fact that the system oscillator Hamiltonian commutes with the system-bath interaction Hamiltonian. The first, pure imaginary term $-i\Omega(n-n^{\prime})t$ in the argument of the exponential is just the free oscillator system evolution. The second pure imaginary term

$\displaystyle i\frac{C\Omega^{2}}{\pi}(n-n^{\prime})(n+n^{\prime}+1)\left[\omega_{c}t-\tan^{-1}(\omega_{c}t)\right]$

(28)

is upper cut-off dependent and comprises both a linear term in system number, which renormalizes the system oscillator frequency $\Omega$, and a quadratic term in system number that is in fact of the same form as the free evolution of a Kerr nonlinear oscillator expressed in the number state basis:

$\displaystyle H=\hbar\Omega a^{\dagger}a+\hbar\Lambda_{\mathrm{kerr}}(a^{\dagger}a)^{2}.$

(29)

Thus, we should properly include a Kerr-type nonlinearity in our starting Hamiltonian (17), with the environmentally induced term $i\frac{C\Omega^{2}}{\pi}(n-n^{\prime})(n+n^{\prime})\omega_{c}t$ renormalizing the nonlinear interaction strength $\Lambda_{\mathrm{kerr}}$. The latter term may be thought of as somewhat analogous to the Newtonian gravitational self-interaction arising from the interaction of a matter system with its gravitating environment. Since we are primarily concerned with decoherence in the present work, we will neglect the quadratic in number term, supposing that it renormalizes an existing Kerr nonlinearity with resulting negligible renormalized coupling strength $\Lambda_{\mathrm{kerr}}$. For $t\gg\omega_{c}^{-1}$, the $\tan^{-1}(\omega_{c}t)$ term in (28) tends to $\pi/2$; this term can be absorbed through a shift in the time coordinate: $t\rightarrow\tilde{t}=t-\pi/(2\omega_{c})$.

Taking into account the system frequency and Kerr nonlinearity renormalizations as just described, Eq. (27) simplifies to

	$\displaystyle|n(t)\rangle\langle n^{\prime}(t)|=|n\rangle\langle n^{\prime}|\exp\bigg{(}-i\Omega(n-n^{\prime})t-\frac{C}{\pi}$
	$\displaystyle\times(n-n^{\prime})^{2}\left(\frac{1}{2}\ln\left(1+t^{2}\omega_{c}^{2}\right)+\ln\left[\frac{\beta\hbar}{\pi t}\sinh\left(\frac{\pi t}{\beta\hbar}\right)\right]\right)\bigg{)},$		(30)

where we have redefined the coupling constant as $C\rightarrow\tilde{C}=C\Omega^{2}$ and dropped the tilde. The real term on the second line of the argument of the exponential results in decoherence, i.e., exponential decay of the outer product for $n\neq n^{\prime}$. In the high temperature (equivalently long time) limit corresponding to $t\gg\beta\hbar\gg\omega_{c}^{-1}$, Eq. (30) can be approximated asymptotically as

	$\displaystyle|n(t)\rangle\langle n^{\prime}(t)|=|n\rangle\langle n^{\prime}|\exp\bigg{(}-i\Omega(n-n^{\prime})t$
	$\displaystyle-(n-n^{\prime})^{2}C\left[\frac{1}{\pi}\ln\left(\frac{\beta\hbar\omega_{c}}{2\pi}\right)+(\beta\hbar)^{-1}t\right]\bigg{)}.$		(31)

From Eq. (31), it is clear that the outer product terms for $n\neq n^{\prime}$ decay exponentially with rate given by $(n-n^{\prime})^{2}Ck_{B}T/\hbar$. Note however, that for early, ‘Planckian’ (by analogy with gravity) times $t\lesssim\omega_{c}^{-1}$, the rate of decoherence is governed by the upper cut-off $\omega_{c}$, resulting in the logarithm term appearing in Eq. (31); depending on the magnitude of the ratio $\hbar\omega_{c}/{k_{B}T}\gg 1$, there may already be a significant ‘burst’ of decoherence during the ‘Planckian’ regime before the later, high temperature exponential decoherence regime. The fact that the decoherence rate depends on the upper cut-off frequency $\omega_{c}$ is a consequence of assuming an initial system-environment product state Anglin et al. (1997). The latter assumption is tantamount to supposing that the system initial state can be prepared on time scales shorter than $\omega_{c}^{-1}$ (or equivalently, the system-environment interaction is switched on over a time scale shorter than $\omega_{c}^{-1}$). While this may be possible for low energy, solid state system environments (i.e., phonons), for an actual gravitational environment with corresponding characteristic Planck time scale, the system state cannot be similarly isolated from the gravitational environment; an analysis which accounts for the system remaining correlated with the environment while its state is being prepared on timescales that are long compared with $\omega_{c}^{-1}$, is expected to result in a subsequent decoherence rate that does not depend on the upper cut-off frequency of the environment.

In the following, we will use Eq. (30) to determine the decoherence dynamics of the oscillator system for initial coherent state superpositions of the form

$\displaystyle\rho_{\mathrm{init}}=N\left(|\alpha_{1}\rangle+|\alpha_{2}\rangle\right)\left(\langle\alpha_{1}|+\langle\alpha_{2}|\right)\otimes\rho_{\mathrm{bath}},$

(32)

where $|\alpha_{1}\rangle$ and $|\alpha_{2}\rangle$ denote coherent states and $N$ is the normalization constant. We consider coherent states since they describe most closely a cooled down, macroscopic oscillator center of mass system. We emphasise that we do not rely on any of the approximations that are often invoked in the study of open quantum system dynamics (beyond assuming an initial product state). In particular, the following analysis is valid for both short/long time scales and high/low temperatures.

Note from the form of (the exact) Eq. (30), that the system will evolve into a classical mixture of number states with probability coefficients that are identical to the coefficients of the initial system state. In contrast to other types of system-bath interaction where the final steady state of the system is usually temperature dependent, for the present model the temperature only determines how fast the system decoheres–not its long time limit steady state.

From Eq. (30), we also see that the decoherence rate is proportional to $(n-n^{\prime})^{2}$ for a superposition of two number states $|n\rangle$ and $|n^{\prime}\rangle$. Thus, for a superposition of coherent states, we expect that the larger the average energy difference, the more rapid the decoherence. This trend is apparent in the oscillator system Wigner function Case (2008) snapshots shown in Fig. 2. For the initial, example superposition state with $\alpha_{1}=3$, $\alpha_{2}=-7$, the negative Wigner function regions disappear in the long time limit (signifying loss of quantum coherence). On the other hand, for the initial example superposition states with nearby coherent state parameter magnitudes: $\alpha_{1}=-\alpha_{2}=3$, $\alpha_{1}=3$ and $\alpha_{2}=-5$, negative Wigner function regions remain in the long time limit (signifying remaining quantum coherence), as is seen more clearly for the zoomed-in Fig. 3. Such trends are consistent with decoherence only resulting for initial spatial superpositions where the states making up the superposition have sufficiently distinct average energies; initial spatial superpositions with the same (or nearby) average energies for the states making up the superposition do not completely decohere.

Note also from the Wigner function snapshots in Fig. 2 and Fig. 3 that the initial coherent superpositions phase-diffuse first into crescent-like regions and then eventually into rings. This is consistent with the fact that, as mentioned above, the final state is always a mixture of number states.

The above findings are in accord with first investigations on the gravitational decoherence of massive scalar quantum field initial superposition states Blencowe (2013); Anastopoulos and Hu (2013), where it is found that superpositions comprising distinct energy states decohere.

Following from the discussion in Sec. I, an operational way (i.e., in principle measurement procedure) to quantify the coherence is through the system oscillator position detection probability density $P(x,t)=\langle x|\rho(t)|x\rangle$ [c.f. the full field-theoretic counterpart Eq. (8)] when the two (initially coherent) wavefunctions making up the superposition pass through each other at $x=0$; these time instants are $\tau_{n}=\Omega t_{n}={\pi}(n+1/2),\,n=0,1,2,\dots$ for the initial coherent state superposition examples considered above, as can be seen for the early time snapshots in Fig 2. The presence of coherence is manifested in $P(x,t)$ having an oscillatory dependence about $x=0$. The latter operational approach corresponds to a two-slit inteference measurement, where the harmonic potential plays the role of the slits by (periodically) bringing the wavefunction components in the initial superposition together. Figure 4 shows the position probability distribution function in the long time limit, steady state for the various example initial coherent state superpositions; we can see that the probability density indicates interference fringes in the vicinity of $x=0$ consistent with the presence of negative-valued Wigner function regions shown in Fig. 2; the snapshots can be interpreted as the marginal probability distributions obtained by integrating over the momentum coordinate Wigner function distributions. In particular, the interference remains for $\alpha_{1}\simeq-\alpha_{2}$, where the average energies of each coherent state making up the initial superposition are not too dissimilar. Note that the other, larger scale scale probability variations in Fig. 4 are due to the final, steady state being a mixture of different number states, as mentioned earlier above.

We adopt the commonly used ‘visibility’ as a measure of the size of the interference fringes, defined as

$\displaystyle\nu=\frac{P_{{\mathrm{max}}}-P_{{\mathrm{min}}}}{P_{{\mathrm{max}}}+P_{{\mathrm{min}}}},$

(33)

where $P_{{\mathrm{max}}}$ is the central maximum of the probability density $P(x)$ at $x=0$ , and $P_{{\mathrm{min}}}$ is the first local minimum of the probability to the right (or left) of the central maximum. The decrease in visibility over time starting from the initial superposition state (32), provides an operational, quantitative measure of decoherence; Fig. 5 gives the visibility as a function of time for various example, initial coherent state, bath temperature, and system-bath coupling parameters. As to be expected, the visibility decreases more rapidly the higher the temperature and the stronger the coupling. Also, the more dissimilar in magnitude $\alpha_{2}<0$ is from $\alpha_{1}>0$ (and hence the larger the average energy difference) in the initial coherent state superposition, the more rapid is the decrease in visibility.

We start with the Lagrangian $L_{\mathrm{qed}}$ (12) and will first derive the Langevin equation that describes the classical oscillator system dynamics interacting with the oscillator bath following the approach of Ref. Cortés et al. (1985). Expanding out the kinetic energy term of Eq. (12), we have

	$\displaystyle L$	$\displaystyle=\frac{1}{2}M\dot{x}^{2}-\frac{1}{2}M\Omega^{2}x^{2}+\sum_{i}\left(\frac{1}{2}m\dot{q}_{i}^{2}-\frac{1}{2}m\omega_{i}^{2}q_{i}^{2}\right)$
		$\displaystyle+Mx\dot{x}\sum_{i}\lambda_{i}q_{i}+\frac{1}{2}Mx^{2}\sum_{i,j}\lambda_{i}\lambda_{j}q_{i}q_{j},$		(34)

where we omit the ‘qed’ subscript from now on. The system and bath momentum coordinates are

	$\displaystyle p=\frac{\partial L}{\partial\dot{x}}=M\dot{x}+Mx\sum_{i}\lambda_{i}q_{i},$		(35)
	$\displaystyle p_{i}=\frac{\partial L}{\partial\dot{q}_{i}}=m\dot{q}_{i},$		(36)

and the model Hamiltonian is

	$\displaystyle H=$	$\displaystyle\frac{p^{2}}{2M}+\frac{1}{2}M\Omega^{2}x^{2}+\sum_{i}\left(\frac{p_{i}^{2}}{2m}+\frac{1}{2}m\omega_{i}^{2}q_{i}^{2}\right)$
	$\displaystyle-$	$\displaystyle xp\sum_{i}\lambda_{i}q_{i}.$		(37)

Hamiltonian (37) is to be compared with the gravity toy model Hamiltonian (16), which differs solely in the form of the system coordinate part of the interaction term; both models have in common a quadratic system coordinate coupling, to be contrasted with the usually studied oscillator system-oscillator bath model with interaction term that is linear in the coupled system and bath coordinates.

Hamilton’s equations for the system and bath coordinates are

	$\displaystyle\dot{p}_{i}$	$\displaystyle=-m\omega_{i}^{2}q_{i}+\lambda_{i}xp,$		(38)
	$\displaystyle\dot{q}_{i}$	$\displaystyle=\frac{p_{i}}{m},$		(39)
	$\displaystyle\dot{p}$	$\displaystyle=-M\Omega^{2}x+p\sum_{i}\lambda_{i}q_{i},$		(40)
	$\displaystyle\dot{x}$	$\displaystyle=\frac{p}{M}-x\sum_{i}\lambda_{i}q_{i}.$		(41)

Formally integrating the equations of motion (38), (39) for the bath coordinates and expressing in terms of the system coordinates:

	$\displaystyle q_{i}(t)$	$\displaystyle-\frac{\lambda_{i}x(t)p(t)}{m\omega_{i}^{2}}=$
		$\displaystyle\left[q_{i}(0)-\frac{\lambda_{i}}{m\omega_{i}^{2}}x(0)p(0)\right]\cos\omega_{i}t+\frac{p_{i}(0)}{m\omega_{i}}\sin\omega_{i}t$
		$\displaystyle-\frac{\lambda_{i}}{m\omega_{i}^{2}}\int_{0}^{t}d\tau\cos\omega_{i}(t-\tau)\frac{d}{d\tau}\left(x(\tau)p(\tau)\right),$		(42)

where we have performed an integration by parts that allows to identify system renormalization and damping terms as we shall see below. Substituting the solution (42) for $q_{i}(t)$ into the equations of motion (40), (41) for the system coordinates leads to the following non-linear Langevin equations:

	$\displaystyle\dot{x}$	$\displaystyle=\frac{\partial H^{m}}{\partial p}+x\int_{0}^{\tau}d\tau K(t-\tau)\frac{d}{d\tau}\left(x(\tau)p(\tau)\right)-xF(t),$		(43)
	$\displaystyle\dot{p}$	$\displaystyle=-\frac{\partial H^{m}}{\partial x}-p\int_{0}^{t}d\tau K(t-\tau)\frac{d}{d\tau}\left(x(\tau)p(\tau)\right)+pF(t),$		(44)

where the renormalized system Hamiltonian is

$\displaystyle H^{m}=\frac{p^{2}}{2M}+\frac{1}{2}M\Omega^{2}x^{2}-\sum_{i}\frac{\lambda_{i}^{2}}{2m\omega_{i}^{2}}x^{2}p^{2},$

(45)

the bath memory kernel is

$\displaystyle K(t-\tau)=\sum_{i}\frac{\lambda_{i}^{2}}{m\omega_{i}^{2}}\cos\omega_{i}(t-\tau),$

(46)

and the bath random force function is

	$\displaystyle F(t)$	$\displaystyle=\sum_{i}\lambda_{i}\left(\left[q_{i}(0)-\frac{\lambda_{i}}{m\omega_{i}^{2}}x(0)p(0)\right]\cos\omega_{i}t\right.$
		$\displaystyle\left.+\frac{p_{i}(0)}{m\omega_{i}}\sin\omega_{i}t\right).$		(47)

In particular, the first term on the right hand side of the equals sign in the Langevin equations (43),(44) describes the Hamiltonian evolution, the second term describes nonlinear damping, and the third term the random force. Note that the bath induces a quartic anharmonic potential in the system Hamiltonian (45). Such a term is analogous to a Coulomb self-interaction potential in the scalar QED field system. After making the rotating wave approximation (RWA), the interaction term reduces to a Kerr-type nonlinearity [c.f. Eq. (29)]. Together with ‘two-photon’ damping [see Eq. (74) below], the resulting open system quantum dynamics can generate quantum states with negative-valued Wigner function regions in the long-time limit, starting from initial Gaussian states Leghtas et al. (2015). In the following, with decoherence dynamics our main subject of interest, we will neglect this induced potential energy term, supposing that it renormalizes an existing anharmonic potential with resulting negligible renormalized coupling strength.

Assuming a thermal equilibrium canonical ensemble distribution for the initial bath coordinates $q_{i}(0),p_{i}(0)$ in Eq. (47), it can be shown that the fluctuation dissipation relation (FDR) between the memory kernel and the random force follows:

$\displaystyle\langle F(t)F(\tau)\rangle=k_{B}TK(t-\tau),$

(48)

where $k_{B}$ is Boltzmann’s constant and $T$ is the bath temperature. We shall assume that the bath responds rapidly on the time-scale of the system oscillator dynamics, so that memory kernel is approximated as $K(t-\tau)=\lambda^{2}k_{0}\delta(t-\tau)$, where $k_{0}$ is a constant. The Langevin equations (43), (44) then become

	$\displaystyle\dot{x}=\frac{p}{M}+\frac{1}{2}\lambda^{2}k_{0}x\frac{d}{dt}(xp)-xF,$		(49)
	$\displaystyle\dot{p}=-M\Omega^{2}x-\frac{1}{2}\lambda^{2}k_{0}p\frac{d}{dt}(xp)+pF,$		(50)

with the FDR (48) taking the form

$\displaystyle\langle F(t)F(\tau)\rangle=k_{B}Tk_{0}\lambda^{2}\delta(t-\tau).$

(51)

The above delta function-approximated memory kernel can be obtained from a bath spectral density $n(\omega)$ with upper cut-off frequency $\omega_{c}$ in the limit $\omega_{c}\rightarrow\infty$. In particular, for a dense bath spectrum, we can approximate the sum over bath degrees of freedom with a bath spectral frequency integral:

$$\sum_{i}\lambda_{i}^{2}\left(\cdots\right)\rightarrow\int_{0}^{\infty}d\omega\lambda^{2}n(\omega)\left(\cdots\right).$$

(52)

Assuming a Lorentzian spectral density

$$n(\omega)=\frac{mk_{0}}{\pi}\frac{\omega^{2}\omega_{c}^{2}}{\omega^{2}+\omega_{c}^{2}},$$

(53)

the memory kernel (46) then becomes

$\displaystyle K(t-\tau)=\frac{\lambda^{2}k_{0}\omega_{c}}{2}e^{-\omega_{c}|t-\tau|}.$

(54)

Taking the infinite limit $\omega_{c}\rightarrow+\infty$, we obtain the above delta function-approximated memory kernel:

$$\lim_{\omega_{c}\rightarrow+\infty}K(t-\tau)=\lambda^{2}k_{0}\delta(t-\tau).$$

(55)

Note that we could equally well have assumed a spectral density with exponential cut-off function instead, as for the gravity toy model [c.f. Eq. (23)]; while the calculations are somewhat more straightforward for the Lorentzian spectral density, we do not expect any qualitative differences in the resulting system quantum dynamics. The motivation to use the Lorentzian spectral density here is purely calculational convenience. The classical, non-linear Langevin equations (49), (50) can be numerically solved as stochastic differential equations as we show in the following sections when comparing with the corresponding quantum dynamics.

The quantum description is obtained through the correspondence principle where $x$, $p$ and $p_{i}$, $q_{i}$ become operators satisfying the canonical commutation relations:

$$[x,p]=i\hbar,~{}[x_{i},p_{j}]=i\hbar\delta_{ij},$$

(56)

with all other commutators vanishing. From Eq. (37), the quantum Hamiltonian operator is

	$\displaystyle H=$	$\displaystyle\frac{p^{2}}{2M}+\frac{1}{2}M\Omega^{2}x^{2}+\sum_{i}\left(\frac{p_{i}^{2}}{2m}+\frac{1}{2}m\omega_{i}^{2}q_{i}^{2}\right)$
		$\displaystyle-\frac{1}{2}\sum_{i}\lambda_{i}q_{i}(xp+px),$		(57)

where the interaction term on the second line is symmetrized in $x$ and $p$ in order that $H$ is Hermitian. Formally integrating Heisenberg’s equations of motion for the bath operators, we obtain the following quantum Langevin equations for the system position and momentum operators:

	$\displaystyle\dot{x}$	$\displaystyle=\frac{p}{M}-\sum_{i}\frac{\lambda_{i}^{2}}{2m\omega_{i}^{2}}x\left(xp+px\right)$
		$\displaystyle+\frac{1}{2}x\int_{0}^{\tau}d\tau K(t-\tau)\frac{d}{d\tau}\left(xp+px\right)-F(t)x,$		(58)

	$\displaystyle\dot{p}$	$\displaystyle=-M\Omega^{2}x+\sum_{i}\frac{\lambda_{i}^{2}}{2m\omega_{i}^{2}}p\left(xp+px\right)$
		$\displaystyle-\frac{1}{2}p\int_{0}^{\tau}d\tau K(t-\tau)\frac{d}{d\tau}\left(xp+px\right)+F(t)p,$		(59)

where the force noise operator is given by Eq. (47) with the system/bath coordinates and momenta replaced by their corresponding operators, and $K(t-\tau)$ is the memory kernel given by Eq. (46). It is convenient to express the quantum Langevin equations in terms of the system creation and annihilation operators which are defined through the usual relations $x=\sqrt{\frac{\hbar}{2M\Omega}}(a+a^{{\dagger}})$, $p=i\sqrt{\frac{M\Omega\hbar}{2}}(a^{\dagger}-a)$:

	$\displaystyle\dot{a}$	$\displaystyle=-i\Omega a-\frac{i\hbar}{2}\sum_{i}\frac{\lambda_{i}^{2}}{m\omega_{i}^{2}}a^{\dagger}\left(a^{{\dagger}2}-a^{2}\right)-F(t)a^{\dagger}$
		$\displaystyle+\frac{i\hbar}{2}a^{\dagger}\int_{0}^{t}d\tau K(t-\tau)\frac{d}{d\tau}\left(a^{{\dagger}2}-a^{2}\right).$		(60)

Under conditions of weak system-environment coupling, Eq. (60) can be simplified by applying the RWA as we now show. Making the substitution $a(t)=A(t)e^{-i\Omega t}$ in Eq. (60), we obtain

	$\displaystyle\dot{A}$	$\displaystyle=-\frac{i\hbar}{2}\sum_{i}\frac{\lambda_{i}^{2}}{m\omega_{i}^{2}}A^{{\dagger}}\left(e^{4i\Omega t}A^{{\dagger}2}-A^{2}\right)$
		$\displaystyle+\frac{i\hbar}{2}A^{{\dagger}}e^{2i\Omega t}\int_{0}^{t}d\tau K(t-\tau)\frac{d}{d\tau}\left(A^{{\dagger}2}e^{2i\Omega\tau}\right.$
		$\displaystyle\left.-A^{2}e^{-2i\Omega\tau}\right)-e^{2i\Omega t}F(t)A^{{\dagger}}.$		(61)

Dropping fast rotating terms, neglecting time derivatives of $A(\tau)$ (since $A$ evolves at much slower rates than $\Omega$), and setting $A(\tau)=A(t)$ (Markov approximation), Eq. (61) becomes approximately

	$\displaystyle\dot{A}$	$\displaystyle=\frac{i\hbar}{2}\sum_{i}\frac{\lambda_{i}^{2}}{m\omega_{i}^{2}}A^{{\dagger}}A^{2}-e^{2i\Omega t}F(t)A^{{\dagger}}$
		$\displaystyle-{\hbar\Omega}\int_{0}^{t}d\tau K(t-\tau)e^{2i\Omega(t-\tau)}A^{{\dagger}}(t)A(t)^{2}.$		(62)

Utilizing the Lorentzian spectral density (53) [with Eq. (49)], Eq. (62) becomes

$$\dot{A}=\frac{i\gamma\omega_{c}}{2\Omega}A^{\dagger}A^{2}-\frac{\gamma\omega_{c}}{\omega_{c}-2i\Omega}A^{{\dagger}}A^{2}-e^{2i\Omega t}F(t)A^{\dagger},$$

(63)

where we have dropped fast oscillating terms and where $\gamma=\hbar\Omega\lambda^{2}k_{0}/2$. For $\omega_{c}\gg\Omega$ and neglecting the anharmonic interaction term, Eq. (63) simplifies to

$$\dot{A}=-\gamma A^{\dagger}A^{2}-e^{2i\Omega t}F(t)A^{\dagger}.$$

(64)

Defining the noise operator as

$$b(t)=\frac{-e^{2i\Omega t}F(t)}{2\sqrt{\gamma}}$$

(65)

and utilizing Eqs. (47), (53) and the RWA, the usual noise operator (anti)commutation rules follow:

$$\begin{split}[b(t),b^{{\dagger}}(t^{\prime})]&=\delta(t-t^{\prime}),\\ \{b(t),b^{{\dagger}}(t^{\prime})\}&=\delta(t-t^{\prime})\big{[}2n(2\Omega)+1\big{]},\end{split}$$

(66)

where the Bose-Einstein thermal average occupation number is evaluated at twice the system oscillator frequency: $n(2\Omega)=\left(e^{2\hbar\Omega/k_{B}T}-1\right)^{-1}$. Finally, transforming back to the non-rotating frame, $A(t)=a(t)e^{i\Omega t}$, we obtain our desired, RWA quantum Langevin equation:

$$\dot{a}=-i\Omega a-\gamma a^{\dagger}a^{2}+2\sqrt{\gamma}e^{-2i\Omega t}ba^{\dagger}.$$

(67)

From Eq. (67), we see that the parameter $\gamma$ has the dimensions of inverse time and characterizes the strength of a nonlinear damping term, while the third term is the nonlinear force noise operator. Equation (67) can be solved numerically as a quantum stochastic differential equation or approximately by first deriving the equations for the various moments in $a$, $b$, and their Hermitian conjugates and truncating at some order.