Production of Entanglement Entropy by Decoherence

TL;DR

This paper rigorously analyzes how decoherence in a two-level system interacting with an oscillator environment leads to entanglement entropy production, identifying key oscillator frequencies involved.

Contribution

It establishes a precise link between decoherence and entanglement entropy production, highlighting the specific environment oscillators responsible.

Findings

Not all environment oscillators contribute equally to entanglement entropy.

Identifies oscillator frequency regions that significantly produce entanglement entropy.

Results are valid for arbitrary interaction strengths and are mathematically rigorous.

Abstract

We examine the dynamics of entanglement entropy of all parts in an open system consisting of a two-level dimer interacting with an environment of oscillators. The dimer-environment interaction is almost energy conserving. We find the precise link between decoherence and production of entanglement entropy. We show that not all environment oscillators carry significant entanglement entropy and we identify the oscillator frequency regions which contribute to the production of entanglement entropy. Our results hold for arbitrary strengths of the dimer-environment interaction, and they are mathematically rigorous.

Equations

$$|{\Psi }_0⟩=|{\psi }_{\mathrm{S}}⟩\otimes |{\alpha }_1⟩\otimes |{\alpha }_2⟩\cdots \otimes |{\alpha }_N⟩,$$

$$|{\alpha }_j⟩={\mathrm{e}}^{-|{\alpha }_j{|}^2/2}\sum _{n\ge 0}\frac{{\alpha }_j^n}{\sqrt{n!}}|n⟩$$

$$H_X=\frac{1}{2}\Omega {\sigma }_z+\sum _{j=1}^N{\omega }_j{a}_j^{†}{a}_j+\lambda {\sigma }_z\otimes \sum _{j=1}^N{g}_j(a_j^{†}+a_j),$$

$$H_D=\frac{1}{2}\Omega {\sigma }_z+\frac{1}{2}V{\sigma }_x+\sum _{j=1}^N{\omega }_j{a}_j^{†}{a}_j+\lambda {\sigma }_z\otimes \sum _{j=1}^N{g}_j{a}_j^{†}{a}_j+\mu {\sigma }_x\otimes \sum _{j=1}^N{f}_j{a}_j^{†}{a}_j.$$

$${\rho }_{\mathrm{S}}(t)={\mathrm{Tr}}_{\mathrm{oscillators}}{\mathrm{e}}^{-\mathrm{i}tH}|{\Psi }_0⟩⟨{\Psi }_0|{\mathrm{e}}^{\mathrm{i}tH},$$

$${\rho }_J(t)={\mathrm{Tr}}_{\mathrm{system}+\mathrm{oscillators}\mathrm{\ell }\mathrm{with}{\omega }_{\mathrm{\ell }}\notin J}{\mathrm{e}}^{-\mathrm{i}tH}|{\Psi }_0⟩⟨{\Psi }_0|{\mathrm{e}}^{\mathrm{i}tH}.$$

$$\big{|}\,[\rho_{\rm S}(t)]_{12}\,\big{|}=|s_{[0,\infty)}(t)|\ \big{|}\,[\rho_{\rm S}(0)]_{12}\,\big{|}\,,$$

$$s_{J}(t)=\left\{\begin{array}[]{ll}{\rm e}^{-4\lambda^{2}\int_{J}h_{X}(\omega)\ \,\frac{1-\cos(\omega t)}{\omega^{2}}\,d\omega},&\mbox{$X$-interaction model}\\ {\rm e}^{-\int_{J}h_{D}(\omega)\,[1-{\rm e}^{2{\rm i}t\lambda g(\omega)}]\,d\omega},&\mbox{$D$-interaction model}\end{array}\right.$$

$$h_X(\omega )=\frac{dn}{d\omega }|g(\omega ){|}^2,\text{mbox}and{h}_D(\omega )=\frac{dn}{d\omega }|\alpha (\omega ){|}^2$$

$$|s_{[0,\infty)}(t)|=\left\{\begin{array}[]{ll}{\rm e}^{-2\pi t\lambda^{2}h_{X}(0)\,[1+o_{t}]},&\mbox{$X$-interaction model}\\ {\rm e}^{-\int_{0}^{\infty}h_{D}(\omega)d\omega}\ \big{[}\,1+O\big{(}(t_{\rm pd}/t)^{2}\big{)}\,\big{]},&\mbox{$D$-interaction model}\end{array}\right.$$

$$\displaystyle{\rm spec}\big{(}\rho_{J}(t)\big{)}$$

$$\displaystyle{\rm spec}\big{(}\rho_{\rm S}(t)\big{)}$$

$$r_J(t)=\sqrt{1-4p(1-p)(1-|s_J(t){|}^2)},$$

$$S_{J}(t)\equiv S\big{(}\rho_{J}(t)\big{)}=-\big{(}\tfrac{1}{2}+\tfrac{1}{2}r_{J}(t)\big{)}\ln\big{(}\tfrac{1}{2}+\tfrac{1}{2}r_{J}(t)\big{)}-\big{(}\tfrac{1}{2}-\tfrac{1}{2}r_{J}(t)\big{)}\ln\big{(}\tfrac{1}{2}-\tfrac{1}{2}r_{J}(t)\big{)}.$$

$$\sum _{j=1}^{N}f({\omega }_j)=\sum _{j=1}^{N^{\prime }}n({\omega }_j)f({\omega }_j)=\sum _{j=1}^{N^{\prime }}\frac{n({\omega }_j)}{\Delta {\omega }_j}f({\omega }_j)\Delta {\omega }_j⟶{\int }_0^{\infty }\frac{dn(\omega )}{d\omega }f(\omega )d\omega ,$$

$$|\psi (t)⟩={\mathrm{e}}^{-\mathrm{i}tH}|\psi ⟩=\sum _{k=1}^d{c}_k|{\varphi }_k⟩\otimes {\mathrm{e}}^{-\mathrm{i}t{H}_k}|{\psi }_{\mathrm{R}}⟩,$$

$${\rho }_{\mathrm{R}}(t)={\mathrm{Tr}}_{\mathrm{S}}|\psi (t)⟩⟨\psi (t)|=\sum _{k=1}^d|c_k{|}^2{\mathrm{e}}^{-\mathrm{i}t{H}_k}|{\psi }_{\mathrm{R}}⟩⟨{\psi }_{\mathrm{R}}|{\mathrm{e}}^{\mathrm{i}t{H}_k}.$$

$${\rho }_{\mathrm{R}}(t)=\sum _{k=1}^d|c_k{|}^2{⨂}_{j=1}^N{\mathrm{e}}^{-\mathrm{i}t{H}_{k,j}}|{\psi }_{\mathrm{R},j}⟩⟨{\psi }_{\mathrm{R},j}|{\mathrm{e}}^{\mathrm{i}t{H}_{k,j}}.$$

$${\rho }_J(t)={\mathrm{Tr}}_{\{1,\dots ,N\}\backslash J}{\rho }_{\mathrm{R}}(t)=\sum _{k=1}^d|c_k{|}^2{⨂}_{j\in J}{\mathrm{e}}^{-\mathrm{i}t{H}_{k,j}}|{\psi }_{\mathrm{R},j}⟩⟨{\psi }_{\mathrm{R},j}|{\mathrm{e}}^{\mathrm{i}t{H}_{k,j}}.$$

$$s_{J}(t)=\left\{\begin{array}[]{ll}\exp\Big{[}-4\lambda^{2}\sum_{\{j\,:\,\omega_{j}\in J\}}|g_{j}|^{2}\,\frac{1-\cos(\omega_{j}t)}{\omega_{j}^{2}}\Big{]},&\mbox{$X$-interaction model}\\ &\\ \exp\Big{[}-\sum_{\{j\,:\,\omega_{j}\in J\}}|\alpha_{j}|^{2}(1-{\rm e}^{2{\rm i}t\lambda g_{j}})\Big{]},&\mbox{$D$-interaction model}\end{array}\right.$$

$${\rho }_J(t)=\left(\begin{array}{cc}p+(1-p)|s_J(t){|}^2& (1-p)s_J\sqrt{1-|s_J(t){|}^2}\\ (1-p){\bar{s}}_J(t)\sqrt{1-|s_J(t){|}^2}& (1-p)(1-|s_J(t){|}^2)\end{array}\right),$$

$$[\rho_{\rm S}(t)]_{12}={\rm e}^{-{\rm i}\Omega t}s_{[0,\infty)}(t)\,\exp\big{\{}4{\rm i}\lambda\,{\rm Im}\,\textstyle\sum_{j=1}^{N}\bar{\alpha}_{j}g_{j}\tfrac{1-{\rm e}^{{\rm i}\omega_{j}t}}{\omega_{j}}\big{\}}\ [\rho_{\rm S}(0)]_{12}\,,$$

$$[{\rho }_{\mathrm{S}}^0(t){]}_{12}={\mathrm{e}}^{-\mathrm{i}\Omega t}{\bar{s}}_{[0,\infty )}(t)[{\rho }_{\mathrm{S}}(0){]}_{12},$$

$$|s_J(t)|={\mathrm{e}}^{-{\int }_J{h}_D(\omega )d\omega }{\mathrm{e}}^{\mathrm{Re}{\int }_J{h}_D(\omega ){\mathrm{e}}^{2\mathrm{i}t\lambda g(\omega )}d\omega }⟶{\mathrm{e}}^{-{\int }_J{h}_D(\omega )d\omega },t\to \infty ,$$

$$h_D(\omega )=\frac{dn(\omega )}{d\omega }|\alpha (\omega ){|}^2.$$

$${\int }_0^{\infty }h(\omega ){\mathrm{e}}^{-2\mathrm{i}t\lambda g(\omega )}d\omega$$

$$t_{\mathrm{pd}}=\frac{\sqrt{|\xi |}}{\lambda },$$

$$\xi=\Big{(}\frac{h_{D}(g^{-1}(\omega))}{g^{\prime}(g^{-1}(\omega))}\Big{)}^{\prime}\,\Big{|}_{\omega=g(0)}.$$

$$\big{|}[\rho^{0}_{\rm S}(t)]_{12}\big{|}=\big{|}[\rho_{\rm S}(0)]_{12}\big{|}\ {\rm e}^{-\int_{0}^{\infty}h_{D}(\omega)d\omega}\ \big{(}1+O(t_{\rm pd}^{2}/t^{2})\big{)}$$

$$\displaystyle\big{\|}\rho_{\rm S}(t)-\rho_{\rm S}^{0}(t)-\tfrac{V}{\Omega}\,{\rm e}^{-\sum_{j=1}^{N}|\alpha_{j}|^{2}}\rho^{1}_{\rm S}(t)\big{\|}_{1}=o(V/\Omega),$$

$$\rho_{\rm S}^{1}(t)=\begin{pmatrix}{\rm Re}\ \big{(}1-{\rm e}^{-{\rm i}\Omega t}\big{)}\ [\rho_{\rm S}(0)]_{12}&(p-\tfrac{1}{2})(1-{\rm e}^{-{\rm i}\Omega t})\\ (p-\tfrac{1}{2})(1-{\rm e}^{{\rm i}\Omega t})&-{\rm Re}\ \big{(}1-{\rm e}^{-{\rm i}\Omega t}\big{)}\ [\rho_{\rm S}(0)]_{12}\end{pmatrix}.$$

Full text

LA-UR-17-20823

1 Introduction

Quantum systems in contact with an environment (reservoir) typically undergo decoherence. An initial pure system state becomes mixed and hence ‘more classical’ due to the system-reservoir interaction [16]. This loss of purity (or, ‘quantumness’) is mediated by entanglement between the system and its environment and can be quantified by the notion of entanglement entropy (EE). First introduced in [4], the examination of EE is receiving enormous attention from researches in different fields of science, ranging from solid state physics and the theory of quantum computation and information to quantum field- and black hole theory [7, 8, 3, 14, 29, 9, 17, 26, 33, 1, 15] (see also references therein).

The goal of the present work is to study the time evolution of EE in a simple open system, a two-level quantum system interacting with a field of quantum oscillators. We adopt two premises:

• (1)

The initial state of the combined total system plus environment is a pure state.

• (2)

The interaction between the system and the environment is (almost) energy conserving, meaning that the system Hamiltonian is (almost) a constant of motion.

Assumption (1) allows us to study quantum correlations between different parts of the system-environment complex using entanglement entropy (see the paragraph below). In the general theory of open quantum systems, reservoirs are often postulated to be initially in thermal equilibrium states, which are mixed states and which do not fall into the category of assumption (1). However, in our view, it is not necessarily clear why one should start off with initial states in thermal equilibrium, as indeed, equilibrium can be shown to emerge by reducing a pure state to a subsystem. In view of this, premise (1) is rather natural. We comment on this phenomenon in Section 4.

It is well known that decoherence happens more quickly, sometimes much more quickly, than relaxation (dynamics of populations) in many open systems [6, 16, 2]. Since relaxation is produced by energy-exchange processes, but decoherence can be brought about by energy-conserving ones, our assumption (2) is reasonable for time scales not exceeding the system decoherence time. It follows that our approach allows a good description of the relation between entanglement production and decoherence, which is the main topic of the present work.

Before giving an outline of our main results in Section 1.1, we explain the notion of EE.

**Entanglement entropy (EE). ** A quantum state (density matrix) $\rho$ is called a pure state if it has rank one, or equivalently, if $\rho=|\psi\rangle\langle\psi|$ for some vector $\psi$. A state which is not pure is called a mixed state. The von Neumann entropy of a state $\rho$ is defined by $S(\rho)=-{\rm Tr}(\rho\ln\rho)$, which is a non-negative number. The equivalence

$\rho$ is pure $\Longleftrightarrow$ $S(\rho)=0$

makes the von Neumann entropy a good measure for the purity of a state. Given a bipartite quantum system $A+B$, a pure state determined by a vector $\psi$ is called entangled if it is not possible to find vectors $\psi_{A}$ and $\psi_{B}$ of the individual subsystems, such that $\psi=\psi_{A}\otimes\psi_{B}$. Equivalently, $\psi$ is called disentangled (or not entangled), or a product state, if $\psi=\psi_{A}\otimes\psi_{B}$ for some subsystem vectors $\psi_{A}$, $\psi_{B}$.

The following is a basic relation between purity and entanglement. Let $\psi$ be a pure state of a bipartite quantum system $A+B$. The reduced states (density matrices) of subsystems $A$ and $B$ are defined by $\rho_{A}={\rm Tr}_{B}|\psi\rangle\langle\psi|$ and $\rho_{B}={\rm Tr}_{A}|\psi\rangle\langle\psi|$ (partial traces). Then the following equivalence holds,

$\rho_{A}$ is pure $\Longleftrightarrow$ $\psi$ is not entangled.

This is also equivalent to $\rho_{B}$ being a pure state. It is important that the state of the system $A+B$ be pure for the above equivalence to hold; indeed if $\rho=\sigma_{A}\otimes\sigma_{B}$ is the product of two mixed states, then the reduced state of $A$, $\rho_{A}=\sigma_{A}$, is not pure. Combining the above two equivalences leads to the following result.

Let $\psi$ be a pure state of a bipartite quantum system $A+B$ and let $\rho_{A}={\rm Tr}_{B}|\psi\rangle\langle\psi|$. Then

$S(\rho_{A})=0$* $\Longleftrightarrow$ $\psi$ is not entangled.*

This motivates the definition of entanglement entropy [4, 28] to be $S(\rho_{A})$, the von Neumann entropy of the reduction of the pure state of a bipartite system to one of its parts. As it turns out, when $\rho_{A}$ and $\rho_{B}$ are the reductions to the subsystems $A$ and $B$, obtained from a pure state of $A+B$, then $S(\rho_{A})=S(\rho_{B})$. So it does not matter of which subsystem state we take the von Neumann entropy, in the definition of entanglement entropy.

1.1 Outline of main results

In this section, we present the initial state and the Hamiltonians of our models (1.1)-(1.3) and we define the reduced density matrices for the dimer (1.4) and (all parts of) the reservoir (1.5). We then present concrete expressions for the dimer decoherence in subsection 1.1.1 and for the dynamics of the entanglement entropy in subsection 1.1.2. In the latter, we show that an increase in the dimer EE is equivalent to an increase in the dimer decoherence. In subsection 1.1.3, we mention applications to a noisy qubit and to the photosynthetic light harvesting complex.

We consider a dimer (a two-level system) in contact with $N$ oscillators. The initial state of the entire system is

[TABLE]

where $|\psi_{\rm S}\rangle=a|\!\uparrow\rangle+b|\!\downarrow\rangle$ (with $|a|^{2}=p$, $|b|^{2}=1-p$) is an arbitrary pure state of the dimer and

[TABLE]

is a single-oscillator coherent state, $\alpha_{j}\in\mathbb{C}$. While we discuss here an environment of oscillators initially in coherent states, many of our formulas are readily expressed for general, pure or mixed initial states for the oscillators (see Section 6.1). We analyze two different models. In the first one, called the $X$-interaction model, the total Hamiltonian is given by

[TABLE]

where $\Omega,\omega_{j}>0$, $\lambda\in\mathbb{R}$ and $\sigma_{z}$ is the Pauli $z$-matrix. The operators $a^{\dagger}_{j}$ and $a_{j}$ are the creation and annihilation operators associated with the oscillator labeled by $j$, satisfying the canonical commutation relations $a_{j}a^{\dagger}_{k}-a^{\dagger}_{k}a_{j}=\delta_{kj}$ (Kronecker delta). The collection of complex numbers $\{g_{j}\}_{j=1}^{N}$ is called the form factor. The size of $|g_{j}|$ measures how strongly the oscillator labeled by $j$ is coupled to the dimer.666To our knowledge, a mathematically controllable theory for the full system-reservoir dynamics in the presence of an additional term $V\sigma_{x}$ in the Hamiltonian (1.2) is not available yet. Nevertheless, the exact non-Markovian master equation for the dimer alone, evolving according to (1.2) with $V\sigma_{x}$ added, and where initially the dimer and reservoir are disentangled and the latter is in thermal equilibrium, has been derived recently in [13]. The second model we discuss is described by a Hamiltonian with a density-, or $D$-interaction, given by

[TABLE]

Here, $\sigma_{x}$ is the Pauli $x$-matrix and we assume in this model that $\Omega,\omega_{j},\lambda,\mu,g_{j},f_{j}>0$.

Note that the dimer energy, $\tfrac{1}{2}\Omega\sigma_{z}$, is a conserved quantity for the dynamics generated by $H_{X}$ as well as for that generated by $H_{D}$ with $V=\mu=0$. Models in which the dimer energy is conserved are called energy conserving. The oscillator energy, $\sum_{j=1}^{N}\omega_{j}a^{\dagger}_{j}a_{j}$, is also conserved for $H_{D}$, but not for $H_{X}$.

We denote by $\rho_{\rm S}(t)$ be the reduced dimer density matrix, obtained by tracing out all oscillator degrees of freedom. It is given by

[TABLE]

where $\Psi_{0}$ is the initial pure dimer-reservoir state (1.1). The matrix elements are denoted by $[\rho_{\rm S}(t)]_{ij}$ ($i=1$ being associated with the eigenvector $|\!\uparrow\rangle$ of $\sigma_{z}$ with eigenvalue $+1$). The reduced oscillator density matrix is defined as follows. Let $J\subseteq\mathbb{R}$ be a window of frequencies of interest. For a given distribution of the oscillator frequencies $\{\omega_{1},\ldots,\omega_{N}\}$, we have, say, $K$ discrete frequencies lying inside $J$. Denote by $\rho_{J}(t)$ the reduced oscillator density matrix of all oscillators with frequencies inside $J$. Similarly to (1.4), $\rho_{J}(t)$ is obtained from the full density matrix by tracing out the degrees of freedom of the dimer and of all oscillators with (the $N-K$ discrete) frequencies lying outside $J$,

[TABLE]

The sets of non-vanishing eigenvalues of $\rho_{\rm S}(t)$ and of $\rho_{J}(t)$ are called the entanglement spectra of the respective reduce density matrices [19].

1.1.1 Dimer decoherence

For finitely many oscillators (frequencies), the dynamics of the dimer and of the oscillators is quasi-periodic in time. As the oscillator frequency spectrum becomes denser, recurrence times are longer and in the (idealized) limit of a continuum of oscillator frequencies (see Section 1.2), one observes irreversible dynamics. Since the system energy is conserved for the Hamiltonians (1.2) and (1.3) with $V=\mu=0$ (the Hamiltonian $\tfrac{1}{2}\Omega\sigma_{z}$ commutes with the total Hamiltonian), the populations of the system – the diagonal matrix elements in the energy basis – are constant in time. In Theorems 2.1 and 2.2 we analyze the dimer energy conserving dynamics of $\rho_{\rm S}(t)$ and $\rho_{J}(t)$, for any frequency window $J$. We show that

• •

The oscillator density matrix $\rho_{J}(t)$ has rank two ($t>0$). Its entanglement spectrum consists of two eigenvalues which we find explicitly for any $t>0$ and any $N$.

• •

The off-diagonal dimer density matrix element evolves according to the law

[TABLE]

where, for any $J\subseteq[0,\infty)$,

[TABLE]

Here,

[TABLE]

and $\tfrac{dn}{d\omega}(\omega)$ is the frequency density of the reservoir initial state at the frequency $\omega$. In the continuous mode limit, the decoherence factor $0\leq|s_{[0,\infty)}(t)|\leq 1$ satisfies (see Propositions 3.1 and 3.2)

[TABLE]

and $o_{t}$ is a quantity converging to zero as $t\rightarrow\infty$. Relations (1.6) and (1.9) show that decoherence is full for the $X$-interaction (provided $h_{X}(0)\neq 0$), meaning that $[\rho_{\rm S}(t)]_{12}\rightarrow 0$ as $t\rightarrow\infty$. For the $D$-interaction, decoherence is only partial and happens on a partial decoherence time scale $t_{\rm pd}$, explicitly given in (2.11) below.

• •

While (1.6) and (1.9) hold for dimer energy conserving models, we examine in Theorems 2.3 and 7.1 the evolution generated by (1.3) when additionally, $V/\Omega$ takes small non-vanishing values. We isolate the main term in the dynamics as well as the first order correction $\propto V/\Omega$, and we control the error on time scales up to $t_{\rm pd}$. We show that the energy-exchange interaction produces a correction to the (strictly positive) entanglement entropy caused by the energy conserving interaction, whose sign fluctuates in time and depends on the initial dimer state (see Section 2.3)

1.1.2 Dynamics of the entanglement entropy

We show in Theorems 2.1 and 2.2 that the (nonzero) spectrum of the reduced density matrices $\rho_{\rm S}(t)$ and $\rho_{J}(t)$, in the energy-conserving situations, are given by

[TABLE]

where

[TABLE]

with $p=|\langle\uparrow\!|\psi_{\rm S}\rangle|^{2}$. Note that (1.11) follows immediately from (1.10). Indeed, by partitioning the whole system into parts $A$=dimer and $B$=oscillators, the entanglement entropy of $A$ automatically equals that of $B$ (see the introduction). The entanglement entropy is thus

[TABLE]

This gives us the following results.

The entanglement entropy $S_{J}$ is a monotonically decreasing function of the decoherence factor $|s_{J}|$. The parameter $|s_{J}|$ (see (1.7)) takes values in the closed interval $[0,1]$. The entanglement entropy $S_{J}$ is monotonically decreasing in the difference of the eigenvalues, $|r_{J}|$, hence by (1.12), $S_{J}$ is a monotonically decreasing function of $|s_{J}|$ (meaning that $\frac{d}{d|s_{J}|}S_{J}(|s_{J}|)<0$). The value of $S_{J}$ is maximal ($S_{J}=\ln 2$) when both eigenvalues are equal, for $r_{J}=0$. It is minimal ($S_{J}=0$) if $r_{J}=1$. In particular, if $p=0$ or $p=1$, then $r_{J}(t)=1$ for all $J$ and all $t$ and no entanglement is ever created. We have $s_{J}(0)=1=r_{J}(0)$ and $S_{J}(0)=0$, reflecting the fact that the initial state is pure and entirely disentangled.

• 2.

Dimer decoherence and dimer EE production. According to (1.6) and (1.9), we say the dimer decoherence is large provided $|s_{[0,\infty)}(t)|$ (off-diagonal) is small. This shows that an increase in the dimer entanglement entropy is necessarily accompanied by an increase in the dimer decoherence. We say the dimer undergoes full decoherence if $\lim_{t\rightarrow\infty}|s_{[0,\infty)}(t)|=0$ and partial decoherence if $\lim_{t\rightarrow\infty}|s_{[0,\infty)}(t)|>0$. If the dimer undergoes full decoherence, then its final EE, $S_{[0,\infty)}(\infty)$, is given by (1.13) with $r_{J}=r(p)=\sqrt{1-4p(1-p)}$. As a function of $p\in[0,1]$, $S_{[0,\infty)}(\infty)$ is symmetric around the point $p=1/2$ and increases monotonically from the value [math] for $p=0$ to the value $\ln 2$ for $p=1/2$. Full decoherence happens in the energy conserving $X$-interaction system (Proposition 3.1), while for the energy-conserving $D$-interaction, we have partial decoherence (Proposition 3.2).

• 3.

**Localization of EE in the oscillators. ** Consider first the $X$-interaction model. We show in Proposition 3.1 (B) below that the value of $|s_{J}(t)|$ for $J=[0,\omega_{1}]$ with arbitrary $\omega_{1}>0$, and hence the entanglement entropy $S(\rho_{J}(t))$, are determined for large times entirely by the coupling function at the frequency zero, $h_{X}(0)$ (see (1.8)). In Proposition 3.1 (C) we show that as the interval $J$ moves away from the low frequency oscillator region, the entanglement entropy of $\rho_{J}$ decreases. More precisely, let $J=[\omega_{0},\omega_{1}]$ with $0<\omega_{0}<\omega_{1}$. When $\omega_{0}$ in increases, the quantity $|s_{J}(t)|$ increases as well and hence the entanglement entropy of $\rho_{J}(t)$ decreases. This shows that for the $X$-interaction, the entanglement entropy in the oscillators is concentrated in the low-frequency region. For the $D$-interaction, a different picture emerges from our Proposition 3.2 below. By point 1. above, $S_{J}$ decreases monotonically in $s_{J}$, which by (1.7) itself is a decreasing function in the size of the considered frequency window. Assuming for the sake of the argument that $h_{D}(\omega)\approx h_{D,0}$ is roughly constant over a region $\omega\in[0,\omega_{c}]$, we see from (1.9) that for large times, we have $s_{J}\approx{\rm e}^{-|J|h_{D,0}}$, where $|J|$ is the length of the interval $J\subset[0,\omega_{c}]$. It follows that the oscillator entanglement entropy $S(\rho_{J}^{0})$ is translation invariant for large times. Namely, for the $D$-interaction, the entanglement entropy in the oscillators with frequencies in $J$ only depends on the size of $J$, not on the location of $J$. This suggests that in the $D$-interaction model, all oscillators contribute to the dynamical process.

1.1.3 Applications

(A) The decohering qubit

In quantum computation and information theory, the dimer represents a qubit and our energy conserving interaction is the coupling to a purely dephasing noise. There is a huge literature on the dynamics on qubits in quantum computation, we mention only the paper [25] and the books [6, 16, 31]. Characteristic times scales for the dephasing of superconducting qubits in laboratories have been pushed to $\tau_{D}\sim 20\mu$s, while the relaxation time scale is \tau_{R}\sim 60$$\mus, [32]. The dynamics can be viewed as purely dephasing, driven by an energy conserving Hamiltonian, for a long time $t\lesssim\tau_{D}$.

(B) The photosynthetic light harvesting complex

In light harvesting antenna complexes of photosynthetic organisms, dimers based on two chlorophyll molecules, Chl$b$ (donor) and Chl$a$ (acceptor), have been proposed to play a critical role in regulating energy (exciton) transfer and dissipation [30, 24, 11, 23, 12, 18]. The chlorophyll molecules are embedded in a protein-solvent environment and form an open quantum system. The chlorophyll dimer is described by a two-level system representing the two excited energy levels of the chlorophyll molecules. The ground states of the chlorophylls are not included in the description, as it is assumed that the time of the excitation transfer between the two excited levels is much shorter than either the time of fluorescence or the time of decay of the exciton to the environment [23]. When a chlorophyll dimer is populated by light absorption, it undergoes the processes of decoherence and relaxation. Decoherence happens much faster than relaxation (evolution of the populations, i.e., diagonal density matrix elements), and so we are in the energy conserving regime for the dynamics and our theory applies. Typical decoherence time scales are $\tau_{D}\sim~{}300$fs while relaxation times are $\tau_{R}\sim 10$ps.

Current technologies allow detailed observations and measurements of decoherence effects in photosynthetic light harvesting complexes [23, 18]. Therefore, a measurement of entanglement entropy in the dimer is experimentally accessible. Since the dimer entanglement entropy equals the entanglement entropy of the environment (see the introduction above), one has experimental access to characteristics of the environment as well.

1.2 Continuous mode limit

The dynamics generated by either of $H_{X}$ or $H_{D}$ is quasi-periodic, and stays so even when $N\rightarrow\infty$, if the oscillator frequencies $\omega_{j}$ form a discrete set. Irreversible dynamics (of the dimer and of the oscillators) emerges only in the presence of a continuum of oscillator frequencies. We will refer to the procedure of taking the oscillator frequencies closer and closer together as the continuous mode limit.777This does not mean that we create ‘modes’ in the sense of delocalized standing waves in the reservoir. Rather, it means that we increase the number of atoms at varying frequencies making up the environment, thus getting a very dense frequency spectrum. The continuous mode limit is taken as follows. Let $f_{j}=f(\omega_{j})$ be a function of the frequency $\omega_{j}$. Then

[TABLE]

where $N^{\prime}$ is the number of distinct, discrete frequencies and $n(\omega_{j})$ is the number of oscillators having the same frequency $\omega_{j}$, $\Delta\omega_{j}=\omega_{j+1}-\omega_{j}$ and $\tfrac{dn(\omega)}{d\omega}$ is the density of oscillators at the given continuous frequency $\omega$. In what follows, when we refer to functions of the continuous frequency, $\alpha(\omega)$, $f(\omega)$, $g(\omega)$, we assume that the continuous mode limit has been performed.

The procedure (1.14) tells us how to obtain the continuous mode limit of suitable sums of a frequency dependent function. However, not all quantities of interest are of this form! For example, the reduced density matrix $\rho_{J}(t)$ of all oscillators with frequencies lying inside $J$ is an operator acting on the Hilbert space $\otimes_{j=1}^{K}{\mathcal{H}}_{\rm osc}$, where ${\mathcal{H}}_{\rm osc}$ is the pure state space of a single oscillator and $K$ counts the number of elements in $J$. We derive in Lemma 6.2 an explicit form of $\rho_{J}(t)$. It cannot be expressed merely in terms of quantities of the form (1.14). What then is the continuous mode limit of the state $\rho_{J}(t)$? This is not an easy question. It is not even clear on what Hilbert space this continuous mode state should act. Indeed, as we make the oscillator frequency spectrum denser and denser, the window $J$ will contain more and more of the discrete frequencies. In the wanted limit, the cardinality $K$ of the number of modes inside $J$ will increase to ${\mathfrak{c}}=|\mathbb{R}|$, the cardinality of $\mathbb{R}$. Then we do not even have a candidate for a Hilbert space, since we do not know how to make sense of $\otimes_{j=1}^{K}{\mathcal{H}}_{\rm osc}$ in this limit! (How should we take a tensor product of uncountably many Hilbert spaces?) Nevertheless, there is a method of constructing the state and its Hilbert space in the continuous mode limit. It is based on the Gelfand-Naimark-Segal construction for linear functionals on $C^{*}$-algebras (i.e., states). This construction has been carried out in [17] for a coherent state reservoir. However, in the present context, as we will show below, the entanglement spectrum of $\rho_{J}(t)$ can be expressed purely in terms of quantities to which we can apply the procedure (1.14). Therefore, the mathematically more involved formalism of [17] is not needed for the considerations of the present work.

2 Dynamics of subsystems for finite $N$

2.1 Energy conserving interactions

**Reduction of the rank of the oscillator density matrix. ** The structure of the reduced density matrices is easy to understand in general for the energy conserving situation. Suppose a $d$-dimensional system is coupled to a reservoir, described by an interacting Hamiltonian $H=H_{\rm S}+H_{\rm R}+H_{{\rm S}{\rm R}}$, where $H_{{\rm S}{\rm R}}$ commutes with $H_{\rm S}$. Let $|\psi\rangle=|\psi_{\rm S}\rangle\otimes|\psi_{\rm R}\rangle$ be a pure system-reservoir initial state. Expanding $|\psi_{\rm S}\rangle=\sum_{k=1}^{d}c_{k}|\phi_{k}\rangle$, where $\{|\phi_{k}\rangle\}_{k=1}^{d}$ is an orthonormal basis of eigenvectors of $H_{\rm S}$, we get

[TABLE]

with $H_{k}$ a Hamiltonian acting on the reservoir alone. Indeed, because $H_{\rm S}$ and $H_{{\rm S}{\rm R}}$ commute, the eigenvectors $|\phi_{k}\rangle$ of $H_{\rm S}$ also diagonalize the system part of the operator $H_{{\rm S}{\rm R}}$, meaning that $H_{{\rm S}{\rm R}}(|\phi_{k}\rangle\otimes|\psi_{\rm R}\rangle)=|\phi_{k}\rangle\otimes(H^{\prime}_{{\rm S}{\rm R},k}|\psi_{\rm R}\rangle)$ for any reservoir state $|\psi_{\rm R}\rangle$, and where $H^{\prime}_{{\rm S}{\rm R},k}$ is an operator acting purely on the reservoir Hilbert space.888 If $H_{{\rm S}{\rm R}}=G\otimes\Phi$ with $G$ and $\Phi$ acting on the system and on the reservoir, respectively, then $H_{{\rm S}{\rm R},k}^{\prime}=g_{k}\Phi$, where $G|\phi_{k}\rangle=g_{k}|\phi_{k}\rangle$. Then $H_{k}=E_{k}+H_{\rm R}+H_{{\rm S}{\rm R},k}^{\prime}$, where $E_{k}$ is the energy of $|\phi_{k}\rangle$ (i.e., $H_{\rm S}|\phi_{k}\rangle=E_{k}|\phi_{k}\rangle$).

It follows from (2.1) that the reduced reservoir density matrix has the form

[TABLE]

Even though $\rho_{\rm R}(t)$ is a density matrix acting on an (in general) infinite-dimensional reservoir Hilbert space, (2.2) shows that its rank is at most $d$, the dimension of the small system. Since the initial state is pure, this fact is also immediately implied by the Schmidt decomposition theorem, which shows that the ranks of $\rho_{\rm S}(t)$ and $\rho_{\rm R}(t)$ are equal, where $\rho_{\rm S}(t)$ is the reduced system density matrix. Assume now that the reservoir consists of $N$ uncoupled subsystems (e.g., oscillators), so that $H_{\rm R}=\sum_{j=1}^{N}H_{{\rm R},j}$, and that the coupling to the system is of the form $\sum_{j=1}^{N}G_{j}\otimes\Phi_{j}$, where $G_{j}$ and $\Phi_{j}$ are coupling operators acting on the system and on the $j$th reservoir subsystem (oscillator), respectively. Then the operator $H_{k}$ in (2.1) has the form $H_{k}=\sum_{j=1}^{N}H_{k,j}$, where $H_{k,j}$ acts on the $j$th reservoir subsystem only (see also footnote 8). For unentangled initial states $|\psi_{\rm R}\rangle=|\psi_{{\rm R},1}\rangle\otimes\cdots\otimes|\psi_{{\rm R},N}\rangle$ expression (2.2) becomes

[TABLE]

Therefore, reducing $\rho_{\rm R}(t)$ to a collection $J$ of reservoir subsystems (oscillators) results in the density matrix

[TABLE]

This shows that the reduced density matrix of any part of the reservoir has rank at most $d$ as well. For general systems with energy exchange, starting with a pure initial state, it is still true that whole reservoir density matrix $\rho_{\rm R}(t)$ has rank $\leq d$ (again, just by the Schmidt decomposition theorem). However, the rank of $\rho_{J}(t)$ might exceed $d$. Indeed, $\rho_{J}(t)$ is obtained from $\rho_{R}(t)$ by tracing out parts of the reservoir degrees of freedom, an operation which may increase the rank (for instance, the partial trace of a pure state is generally a mixed state, hence taking the partial trace generally increases the rank). This can also be understood from a different point of view: switching on an energy exchange term (even a small one) shifts the spectrum of $\rho_{J}(t)$ and may cause some of its zero eigenvalues to become nonzero, hence increasing the rank of $\rho_{J}(t)$.

**Energy conserving dynamics. ** We examine the reduced dynamics of the dimer and (parts of) the reservoir under the dynamics generated by (1.2) and (1.3) for $V=\mu=0$. Given $J\subseteq[0,\infty)$, we define the quantity (the discrete frequency analogue of (1.7))

[TABLE]

We prove in Theorems 2.1, 2.2 below that the reduced reservoir density matrix of all oscillators with frequencies inside $J$ is a rank-two operator, described in a suitable basis by the matrix

[TABLE]

with $p$ and $s_{J}(t)$ given after (1.1) and in (2.5), respectively. Note that the reservoir entanglement spectrum is entirely determined by (2.6).

Theorem 2.1** (Energy conserving $X$-interaction)**

Consider the system with the $X$-interaction, (1.2) for a fixed, finite $N$.

• (A)

Let $\rho_{\rm S}(t)$ be the reduced dimer density matrix. Its diagonal is time-independent and the off-diagonal is given by

[TABLE]

where $s_{[0,\infty)}(t)$ is given in (2.5) ($X$-interaction) with $J=[0,\infty)$.

• (B)

Let $\rho_{J}(t)$ be the reduced density matrix of all oscillators having frequencies inside $J$. The rank of $\rho_{J}(t)$ is two ($t>0$) and in a suitable orthonormal basis, it has the matrix representation (2.6).

The basis in which $\rho_{J}(t)$ takes the form (2.6) depends on time $t$, the number $N$ of oscillators and the frequency window $J$. For the energy conserving $D$-interaction ($V=\mu=0$) we denote the reduced dimer and oscillator density matrices by $\rho^{0}_{\rm S}(t)$ and $\rho_{J}^{0}(t)$. The superscript [math] indicates that these are the unperturbed states, as later on, we will switch on the energy exchange interactions.

Theorem 2.2** (Energy conserving $D$-interaction)**

Let $N$ be fixed and finite. Denote by $\rho^{0}_{\rm S}(t)$ and $\rho_{J}^{0}(t)$ the reduced system and oscillator density matrices evolving according to the $D$-interaction, (1.3), for $V=\mu=0$.

• (A)

The diagonal of $\rho^{0}_{\rm S}(t)$ is time-independent and the off-diagonal is

[TABLE]

where $\bar{s}_{[0,\infty)}(t)$ is the complex conjugate of $s_{[0,\infty)}(t)$, given in (2.5) ($D$-interaction) with $J=[0,\infty)$.

• (B)

The rank of $\rho^{0}_{J}(t)$ is two ($t>0$) and in a suitable orthonormal basis, it has the matrix representation (2.6).

2.2 $D$-interaction: partial decoherence

As mentioned in Section 1.2, irreversible dynamics emerges in the continuous mode limit only. An example is the evolution of coherences, encoded in the time dependence of $s_{J}(t)$, (2.5). For the $D$-interaction model, we obtain in the continuous mode limit (1.14),

[TABLE]

where

[TABLE]

The convergence is due to fast oscillations, causing $\int_{J}h_{D}(\omega){\rm e}^{2{\rm i}t\lambda g(\omega)}d\omega\rightarrow 0$ as $t\rightarrow\infty$, a process we can control by a ‘Riemann-Lebesgue Lemma’ argument as follows. For definiteness, consider $h$ to be a function of $\omega\in[0,\omega_{c}]$, where $\omega_{c}$ is some cutoff (at and beyond which $h$ and $h^{\prime}$ vanish) and suppose that the coupling function $g$ is invertible on $\omega\in[0,\omega_{c}]$. Using that ${\rm e}^{-2{\rm i}t\lambda\omega}=\frac{1}{-2{\rm i}t\lambda}\partial_{\omega}{\rm e}^{-2{\rm i}t\lambda\omega}$ and integrating by parts, we expand

[TABLE]

It follows from (2.10) and (2.8) that $|[\rho^{0}_{\rm S}(t)]_{12}|$ decays partially (but not monotonically in time) over the partial decoherence time

[TABLE]

with $\xi$ given by

[TABLE]

We point out that the decoherence is only partial, since

[TABLE]

does not vanish as $t\rightarrow\infty$. In contrast, we will show that for the $X$-interaction model, we have full decoherence (see Proposition 3.1).

2.3 $D$-interaction: almost energy conserving regime

Even for $V$ and $\mu$ small (not vanishing), the energy exchange interaction terms in $H_{D}$, (1.3), have a big influence over long time scales. For instance, they are responsible for the process of relaxation (thermalization, if the reservoir is in a thermal state). We consider here time scales of the order $t_{\rm pd}$, (2.11), over which partial decoherence happens and we show that the effects of the non-energy conserving terms are merely a perturbation of this dynamics for small $V,\mu$. We call this the almost energy conserving regime.

In the result below, we use the notation $o(x)$ for a real number s.t. $\lim_{x\rightarrow 0}o(x)/x=0$. Also, we denote by $\|\cdot\|_{1}$ the trace norm, namely, $\|A\|_{1}={\rm Tr}|A|$, where $|A|=\sqrt{A^{*}A}$ is the absolute value of an operator $A$.

Theorem 2.3** (Energy exchange $D$-interaction)**

Let $N$ be finite and consider the $D$-interaction with $\mu=0$ and $V\geq 0$. Suppose that $V/\Omega<\!\!<1$, $V/\lambda<\!\!<1$ and that $\Omega/\lambda<\!\!<\min_{1\leq j\leq N}g_{j}$. Then, for times $0\leq t\leq t_{\rm pd}$, we have the following.

• (A)

The dimer density matrix satisfies

[TABLE]

where $\rho_{\rm S}^{0}(t)$ is given in (2.8) and

[TABLE]

The remainder is uniform in $t$ for $0\leq t\leq t_{\rm pd}$.

• (B)

The oscillator density matrix satisfies

[TABLE]

where $\rho_{J}^{00}(t)$ is the $3\times 3$ matrix obtained by adding a (third) row and a (third) column consisting of zeroes to the matrix in (2.6), and, expressed in the same basis, the hermitian $\rho_{J}^{1}(t)$ has matrix elements

[TABLE]

The parameters are given by

[TABLE]

We refer to Theorem 7.1 below for results when in addition, $\mu>0$.

3 Dynamics in the continuous mode limit

3.1 Dynamics of the dimer and of $s_{J}$

In the continuous mode limit (1.14), the discrete mode expressions for $s_{J}(t)$, (2.5), become (1.7).

Proposition 3.1** ($X$-interaction)**

The dynamics of $\rho_{\rm S}$ and $s_{J}$ for the $X$-interaction model, in the continuous mode limit, satisfies the following. In what follows, $o_{t}$ denotes a real number satisfying $\lim_{t\rightarrow\infty}o_{t}=0$.

• (A)

Suppose that $h_{X}(\omega)$, (1.8), is continuous as $\omega\rightarrow 0_{+}$. Then

[TABLE]

• (B)

Let $J=[0,\omega_{1}]$, where $0<\omega_{1}\leq\infty$. Then for all $t\geq 0$,

[TABLE]

• (C)

Let $J=[\omega_{0},\omega_{1}]$, where $0<\omega_{0}<\omega_{1}<\infty$ and set $m_{J}=\inf_{\omega\in J}h_{X}(\omega)$, $M_{J}=\sup_{\omega\in J}h_{X}(\omega)$. We have for all $t\geq 0$,

[TABLE]

According to (A), the dimer undergoes full decoherence for couplings $h_{X}(0)\neq 0$, exponentially quickly in time.

Proposition 3.2** ($D$-interaction, $V=\mu=0$)**

The dynamics of $\rho_{\rm S}$ and $s_{J}$ for the energy conserving $D$-interaction model, in the continuous mode limit, satisfies the following.

• (A)

We have

[TABLE]

where $t_{\rm pd}\propto 1/\lambda$ is the partial decoherence time (c.f. (2.11)).

• (B)

For any $J\subseteq[0,\infty)$, the quantity $s_{J}(t)$, given in (1.7), satisfies

[TABLE]

The expansions are correct provided $t_{\rm pd}>0$, which translates into a condition of effective coupling, see (7.7), (7.8). The dimer undergoes partial decoherence.

We close the discussion by explaining how the energy-exchange term $\propto V$ of the $D$-interaction model, (1.3), influences the entanglement entropy. For some results including additionally the term $\propto\mu$, see Section 7.

Proposition 3.3** (Entanglement entropies, $D$-interaction, $V>0$)**

Consider the dynamics associated with the $D$-interaction Hamiltonian, for $V>0$ (and $\mu=0$), in the continuous mode limit. Define

[TABLE]

Suppose that $V/\Omega<\!\!<1$, $V/\lambda<\!\!<1$ and that $\Omega/\lambda<\!\!<\inf_{\omega\geq 0}g(\omega)$. Then, for times $0\leq t\leq t_{\rm pd}$, we have:

• (A)

The eigenvalues of $\rho_{\rm S}(t)$ are

[TABLE]

• (B)

The eigenvalues of $\rho_{J}(t)$ are, in the continuous mode limit,

[TABLE]

The small-o notation means that $o(x)$ is a real number s.t. $\lim_{x\rightarrow 0}o(x)/x=0$. Also, $r_{J}(t)$ is given in (1.12).

**Discussion. ** 1. Theorem 2.3 and Proposition 3.3 show that the energy-exchange perturbation $\propto V$ in (1.3) leaves the rank of the oscillator density matrices $\rho_{J}(t)$ unchanged, up to order $V/\Omega$ (see (3.8)).

2. For $p=1/2$, the term on the right side of (3.7), which is linear in $V/\Omega$, vanishes. The entanglement entropy of the dimer decreases if the gap between the two eigenvalues of $\rho_{\rm S}(t)$ increases, which happens exactly if the term $\propto V/\Omega$ in the parentheses in (3.7) has a positive sign. This shows that the energy-exchange interaction produces a correction to the (strictly positive) entanglement entropy caused by the energy conserving interaction, whose sign fluctuates in time and depends on the initial dimer state.

4 Energy exchange, dissipation, equilibration

4.1 Systems close to thermal equilibrium

The theory of open quantum systems analyzes the dynamics of ‘small systems’ coupled to ‘environments’ (also called reservoirs or quantum noises). The environment itself is a quantum system but it shows irreversible dynamical effects, both in its own evolution and in that of a small system coupled to it. A reservoir is obtained by a limiting procedure from a ‘normal’ (finite) quantum system, for instance by taking a continuous mode limit. The limits turn purely discrete energy spectrum of finite systems into continuous spectrum, which causes irreversibility in the dynamics. It is plausible to expect that a system with continuous energy (frequency) spectrum shows dissipative dynamics, suggesting that reservoirs have a unique stationary state. The same thing happens when a transition of discrete to continuous frequencies is caused by a thermodynamic limit, where the volume of a reservoir increases, see [17, 22, 5].

A non-interacting system environment complex, described by a Hamiltonian $H_{0}=H_{\rm S}+H_{\rm R}$, has a manifold of stationary states, spanned by $|\psi_{E}\rangle\langle\psi_{E}|\otimes\rho_{\rm R}$, where $\psi_{E}$ are the system eigenfunctions and $\rho_{\rm R}$ is the reservoir invariant state. What happens as the system and reservoir start to interact, due to a term $H_{{\rm S}{\rm R}}$ in the Hamiltonian? In the situation where the reservoir is in (or close to) thermal equilibrium (the thermodynamic limit of finite-volume Gibbs states), one can show that, for very generic interactions $H_{{\rm S}{\rm R}}$, the interacting system has a unique stationary state, which is the coupled (interacting) equilibrium state. The fact that such a coupled equilibrium exists is guaranteed by the stability theory of KMS states.999For continuous modes (frequencies) systems, the Hamiltonian $H_{\rm R}$ has continuous spectrum and one does not know how to interpret the candidate ${\rm e}^{-\beta H_{\rm R}}$ for an equilibrium density matrix, since its trace is infinite. A better characterization of equilibrium is the KMS characterization (Kubo-Martin-Schwinger) [5], which reduces to the usual Gibbs state formulation for equilibrium in finite systems. This theory says, quite generally, that if an unperturbed Hamiltonian $H_{0}$ generates a dynamics which has an equilibrium (KMS) state $\rho_{0}$ then a perturbation $H_{0}+H_{{\rm S}{\rm R}}$ generates a dynamics with associated equilibrium state $\rho_{{\rm S}{\rm R}}$, which varies continuously in $H_{{\rm S}{\rm R}}$. The fact that $\rho_{{\rm S}{\rm R}}$ is the only stationary state can be proven for many open systems, provided an efficient coupling condition is satisfied. The latter requires that energy exchange takes place between the system and the environment, and it will not be satisfied in general for models where the system energy is conserved.101010Suppose $H=H_{\rm S}+H_{\rm R}+H_{{\rm S}{\rm R}}$ is energy conserving, meaning that $[H_{\rm S},H_{{\rm S}{\rm R}}]=0$. Then all initial states $|\psi_{E}\rangle\langle\psi_{E}|\otimes\rho_{\rm R}$, where $H_{\rm S}\psi_{E}=E\psi_{E}$ and $\rho_{\rm R}$ is an arbitrary reservoir density matrix produce stationary system states, ${\rm Tr}_{\rm R}\{{\rm e}^{-{\rm i}tH}(|\psi_{E}\rangle\langle\psi_{E}|\otimes\rho_{\rm R}){\rm e}^{{\rm i}tH}\}=|\psi_{E}\rangle\langle\psi_{E}|$. This observation hints at the need of energy exchange to guarantee uniqueness of a stationary state for the whole, coupled system-reservoir complex. We sum up the situation close to equilibrium as follows:

The coupled system-reservoir complex has an equilibrium state.

• 2.

Due to efficient coupling, requiring in particular energy-exchange, the coupled equilibrium is the only stationary state.

• 3.

The coupled dynamics is dissipative and drives any initial state towards the coupled equilibrium.

For small $H_{{\rm S}{\rm R}}$ the reduction of the coupled equilibrium to the system alone is the usual Gibbs state $\propto{\rm e}^{-\beta H_{\rm S}}$ (modulo $O(H_{{\rm S}{\rm R}})$ terms) at the temperature inherited from the reservoir. For strong system-reservoir coupling, the reduced state is a Gibbs state (again at inverse temperature $\beta$) with respect to a Hamiltonian having renormalized energies [20]. The renormalizations are $\propto(H_{{\rm S}{\rm R}})^{2}$ and can be large. In either case, the small system becomes equilibrated due to the contact with a thermal reservoir, in the course of time.

We have developed the dynamical resonance theory [21], allowing a rigorous proof of 1.-3., for initial states of the form $\rho_{\rm S}\otimes\rho_{\rm R}$ (or perturbations thereof), where $\rho_{\rm R}$ is the reservoir equilibrium, hence a mixed state. However, for the analysis of entanglement entropy, one has in mind pure system-reservoir initial states. A modification at the very core of the resonance theory is necessary to make it applicable to the study of entanglement entropy. While we are working to expand the theory in this direction, the present paper investigates the almost energy conserving situation only.

4.2 Emergence of Gibbs states from reducing pure states

In the “standard setup of open systems” described in Section 4.1, the starting point is a state where the reservoir is in (or close to) equilibrium. A different question is how a system equilibrium (Gibbs) state emerges due to the interaction with a reservoir when initially the system-reservoir state is pure.

Suppose $H=H_{\rm S}+H_{\rm R}+H_{{\rm S}{\rm R}}$ where the reservoir has discrete, but close lying energy levels, spaced by maximally $\Delta B$, while the system energy spacing is $\Delta\epsilon$. Denote by $\rho(B)$ the reservoir density of states and define the effective temperature by $\beta\equiv\beta(B)=\tfrac{d}{dB}\ln\rho(B)$. In the regime $\Delta B<\!\!<\|H_{{\rm S}{\rm R}}\|<\!\!<\Delta\epsilon$ (and under additional assumptions), it is shown in [33, 34] that for an initial pure state $\Phi(0)$ of the system plus reservoir, having a (coupled) energy distribution peaked at a value $E$, one has

[TABLE]

for all system observables $A$, where $\Phi(t)={\rm e}^{-{\rm i}tH}\Phi(0)$ and $\langle A\rangle_{\beta}^{\rm can}={\rm Tr}_{\rm S}({\rm e}^{-\beta H_{\rm S}}A)/{\rm Tr}_{\rm S}{\rm e}^{-\beta H_{\rm S}}$ is the system Gibbs canonical equilibrium state. Relation (4.1) is valid for sufficiently large, typical times $t$ (it does not hold for all times, since the reservoir is not taken to be in the continuous mode limit and hence the overall dynamics is quasi-periodic). It would be interesting to give a proof of (4.1) for our $X$-interaction model with a coherent state initial reservoir, especially for strong system-reservoir couplings (relevant, for instance, in quantum biological processes). In view of the techniques developed in [20], we might expect such a result to hold, with $H_{\rm S}$ replaced with a renormalized Hamiltonian.

The following two arguments are often used to heuristically derive the emergence of the system Gibbs state, obtained by reducing a full (pure state) density matrix.

(Decoherence in the energy basis.) After the system and reservoir have interacted for some time, the system will be in a state which is diagonal in its (uncoupled) energy basis, and which does not depend on the initial system state. It follows from this (assumption, or principle) that the system state is entirely described by the populations (occupation probabilities) $p_{i}$ of the energy levels, and that those populations are functions of the energy only, $p_{i}=f(E_{i})$. It is implicitly assumed here that the system and reservoir exchange energy since if they were not, then the system populations would be time-independent and they would not settle to values independent of the initial system state.

• 2.

When coupling two identical, independent systems to the same reservoir (no direct interaction), the population probabilities of the joint two-system complex should split into products of single system probabilities, due to their independence. Here, it is assumed that the systems-reservoir interaction is weak, so that independence of the two systems is approximately preserved even while they are in contact with the same reservoir. This leads to the constraint $f(x+y)=f(x)f(y)$, whose solution has the Gibbs distribution form $f(x)\propto{\rm e}^{-\beta x}$, where $\beta$ is a reservoir-dependent constant.

In the arguments of [33, 34], decoherence is derived by observing that for the models considered there, the stationary Schrödinger equation of the interacting system can be interpreted as that of a single quantum particle in a potential. A weak system-energy interaction makes this potential vary slowly in space and as a consequence, eigenstates of the interacting system associated with energies differing by more than the interaction energy have roughly disjoint support (quasiclassical analysis), which leads to decoherence (c.f. equations (6) in [33]). It is not clear to us if this reasoning is applicable to situations we are interested in here, namely a dimer strongly coupled to an environment.

While the above considerations are based on the dynamics, there is also a ‘static’ link between reduced states of a large system and Gibbs states. In particular, for a collection of solvable many-particle free lattice particle systems and spin chain models in their ground state, the reduced density matrices of parts of the system and their entanglement spectra have been analyzed in [26, 27, 28, 10]. The key finding is that the reduced density matrix $\rho_{\alpha}$ (where $\alpha$ labels a lattice subregion) is of thermal form, $\rho_{\alpha}\propto{\rm e}^{-H_{\alpha}}$, where $H_{\alpha}$ is an effective, local, free particle Hamiltonian, but $H_{\alpha}$ is not the original Hamiltonian restricted to the region $\alpha$. In these models, a thermal structure thus emerges generically by reduction of the ground state to a subregion of space.

5 Numerical simulations

In our numerical simulations we use Eq. (1.13) for entanglement entropy of the environment, written as

[TABLE]

where $r(t)=\sqrt{1-4p(1-p)(1-|s(t)|^{2})}$ and

[TABLE]

We adopt units in which $\hbar=1$.

5.1 $\mathbf{X}$-interaction

Let us consider the function $h_{X}(\omega)\equiv h(\omega)$, given in (1.8), of the form

[TABLE]

with an exponential high frequency cutoff $\omega_{c}>0$, where $A_{q}>0$ and where $p\in\mathbb{R}$ determines the low frequency behavior. It is convenient to introduce a new variable, $z=\omega/\omega_{c}$, and a dimensionless time, $\tau=\omega_{c}t$. Then, Eq. (5.2) can be written as,

[TABLE]

where $\varepsilon=4\lambda^{2}A_{q}\omega_{c}^{2q+1}$ and

[TABLE]

Below, we consider two cases: $q>-1/2$ and $q=-1$. Performing the integration in (5.5) gives the following.

• •

Case: $q>-1/2$

[TABLE]

where $\Gamma(u)=\int_{0}^{\infty}{\rm e}^{-z}z^{u-1}{\rm d}z$ is the gamma function, defined for $u\in\mathbb{C}$ with $\Re u>0$.

• •

Case: $q=-1$

[TABLE]

For both these cases, the EE, $S_{E}(\tau)$, is presented in Fig. 1. For $q=1$, we have: $h(0)=0$, and we see from Fig. 1a that only partial decoherence occurs in the X-interaction model. The function $|s(t)|$ does not decay to zero for large times. Correspondingly, the EE $S_{E}(t)$ does not approach its maximum for a given value of $p$. In particular, for $p=1/2$, $\lim_{\tau\rightarrow\infty}S_{E}(\tau)\approx 0.3$. At the same time, for $q=-1$ and $p=1/2$, the EE reaches its maximum: $\lim_{\tau\rightarrow\infty}S_{E}(\tau)=\ln 2\approx 0.7$, as in this case full decoherence takes place ($\lim_{\tau\rightarrow\infty}s(\tau)=0$).

5.1.1 Contribution of low and high frequencies to the EE

To estimate the contribution of low and high frequencies, we consider lower, $\omega_{0}$, and upper, $\omega_{c}$, cutoff frequencies in the function $\mathcal{Q}_{2}(\tau)$ by introducing,

[TABLE]

where $z_{0}=\omega_{0}/\omega_{c}$. We have ${\mathcal{Q}_{2}}(\tau)={\mathcal{Q}_{2}}(\tau,z_{0})+\tilde{\mathcal{Q}_{2}}(\tau,z_{0})$, but the EE, $S_{E}(\tau)$, is of course not a linear function of ${\mathcal{Q}_{2}}$. In Fig. 2 we evaluate $S_{E}(\tau)$, given in (5.1), in which the parameter $|s(\tau)|$ is chosen to be either $|s(\tau)|={\rm e}^{-\varepsilon\mathcal{Q}_{2}(\tau,z_{0})}$ or $|\tilde{s}(\tau)|={\rm e}^{-\varepsilon\tilde{\mathcal{Q}_{2}}(\tau,z_{0})}$ (the low frequency or high frequency contributions to the EE, respectively).

5.1.2 Contribution of low and high frequencies for $h(\omega)=h_{0}$ inside of the frequency interval $(\omega_{0},\omega_{1})$

Let us consider the function $h(\omega)$ given by,

[TABLE]

where $\Theta(x)$ is the Heaviside step function. We define $\omega_{h}=2\lambda^{2}\pi h_{0}$, and introduce the following dimensionless variables: $\nu_{0}=\omega_{0}/\omega_{h}$, $\nu_{1}=\omega_{1}/\omega_{h}$ and $\tau=\omega_{h}t$. In the new variables we have $|s(\tau)|={\rm e}^{-(1/\nu_{0})\mathcal{Q}_{2}(\tau)}$, where

[TABLE]

We use the asymptotic formula ${\rm Si}(z)\sim\frac{\pi}{2}-\frac{\cos z}{z}-\frac{\sin z}{z^{2}},\quad z\gg 1$, to obtain that for $\nu_{0}\tau,\nu_{1}\tau\gg 1$,

[TABLE]

For fixed $\delta=\nu_{1}-\nu_{0}$, we study numerically the behavior of the functions $|s(\tau,\nu_{0})|={\rm e}^{-(1/\nu_{0}){\mathcal{Q}_{2}}(\tau,\nu_{0})}$ and ${\mathcal{Q}_{2}}(\tau,\nu_{0})$, where

[TABLE]

Fig. 3 shows the result for $\delta=0.1$. The entanglement entropy as a function of $\tau$ and $p$ is depicted in Fig. 4, for a fixed value of the frequency domain, $\delta=\nu_{1}-\nu_{0}$, and for increasing value of $\nu_{0}$. These results show hat the small frequencies, near $\omega=0$ give the main contribution to the EE in the X-interaction model.

5.2 $\mathbf{D}$-interaction

For the $D$-interaction model, the function, $s(t)$, given by Eq. (1.7), is written here as,

[TABLE]

For illustrative purposes, in our numerical simulations we choose $h_{D}(\omega)=h_{0}(\Theta(0)-\Theta(\omega_{c}))$ and $g(\omega)=g_{0}\omega^{k}$, where $k>0$. Performing the integration in (5.12), we obtain,

[TABLE]

where $\mu=2\lambda g_{0}h_{0}\omega^{k}_{c}$, $\varepsilon=h_{0}2\lambda g_{0}$, $\tau=h_{0}t$ is the dimensionless time and $\gamma(x,\alpha)=\int_{0}^{\alpha}e^{-z}z^{x-1}{\rm d}z$ denotes the incomplete gamma function. Writing $\gamma(x,\alpha)=\Gamma(x)-\Gamma(x,\alpha)$ and using $\Gamma(a,z)\sim z^{a-1}e^{-z}$, we obtain

[TABLE]

The entanglement entropy vs $\tau$ and $p$ is depicted in Fig. 5. As one can see, in this case, only partial decoherence occurs, and the function, $s(t)$, does not decay to zero. Correspondingly, the EE does not reach its potential maximum $\ln 2\approx 0.7$ for any values of parameters.

6 $X$-interaction, proofs

6.1 Finite mode reservoirs

The Hamiltonian $H_{X}$, (1.2), leaves the dimer populations stationary. Namely, setting

[TABLE]

where $P_{1}$ is the projection onto the ‘donor level’ (state $\varphi_{1}$ with eigenvalue $+1$ of $\sigma_{z}$) and $P_{2}$ is the projection onto the ‘acceptor level’ $\varphi_{2}$, the commutators $[H,P_{1}]=[H,P_{2}]=0$ vanish. We take an initial state of the combined system which is pure and completely disentangled, of the form

[TABLE]

Here,

[TABLE]

is an arbitrary pure initial state of the dimer, with $|a|^{2}+|b|^{2}=1$ and $\psi_{j}$ is a pure initial state of the $j$th oscillator. Below, we specialize to the coherent state reservoir initial state (1.1). We set

[TABLE]

The diagonal entries of the density matrix $|\psi_{\rm S}\rangle\langle\psi_{\rm S}|$ are $p$ and $1-p$.

6.1.1 Reduced states

The state of the total system at time $t$ is

[TABLE]

where

[TABLE]

By a polaron transformation, we can diagonalize $U_{j}^{\pm}(t)$ as follows,

[TABLE]

where $D_{j}(z)$, $z\in\mathbb{C}$, are the (coherent state) displacement operators

[TABLE]

(Note that the second line in (6.7) is obtained from the first one by switching $\lambda\rightarrow-\lambda$.)

Lemma 6.1** (Reduced dimer state)**

The reduced dimer density matrix has constant diagonals and (with $[\rho_{\rm S}(0)]_{12}=a\bar{b}$),

[TABLE]

*Proof. * The reduced density matrix of the dimer is obtained immediately from (6.5),

[TABLE]

Using (6.7) and simply writing $D_{j}$ for $D_{j}(\lambda g_{j}/\omega_{j})$ and $S_{j}$ for $S_{j}(t)$,

[TABLE]

The displacement operators satisfy

[TABLE]

so that $(D_{j}^{*})^{2}=D_{j}(-2\lambda g_{j}/\omega_{j})$. Next, ${\rm e}^{{\rm i}\omega_{j}ta^{\dagger}_{j}a_{j}}D_{j}(z){\rm e}^{-{\rm i}\omega_{j}ta^{\dagger}_{j}a_{j}}=D_{j}({\rm e}^{{\rm i}\omega_{j}t}z)$ and hence $S_{j}^{*}(D^{*}_{j})^{2}S_{j}=D_{j}(-2{\rm e}^{{\rm i}\omega_{j}t}\lambda g_{j}/\omega_{j})$. Then using again (6.12) twice we find

[TABLE]

Using expression (6.13) in (6.11) gives (6.9). This shows Lemma 6.1 $\blacksquare$

Next we obtain the reduced density matrix of the collection of operators with given indices.

Lemma 6.2** (Reduced oscillators state)**

Let $J\subseteq\mathbb{R}_{+}$ and denote by $\rho_{J}(t)$ the reduced density matrix of the oscillators having frequencies inside $J$, obtained from the total density matrix by tracing out all other oscillator degrees of freedom as well as those of the dimer. Then

[TABLE]

where $p=[\rho_{\rm S}(0)]_{11}$ and

[TABLE]

Depending on time $t$, the matrix $\rho_{J}(t)$ has rank either one or two. It has rank one exactly when

[TABLE]

Apart from $t=0$, condition (6.15) is difficult to satisfy, so generically, the rank of $\rho_{J}(t)$ is two.

*Proof. * The reduced density matrix in question is obtained from (6.5),

[TABLE]

This shows the form of $\rho_{J}(t)$ as given in Lemma 6.2. Since the vectors $\Psi^{\pm}_{J}(t)$ are normalized, their span is one-dimensional exactly if $|\left\langle{\Psi_{J}^{-}(t)},{\Psi_{J}^{+}(t)}\right\rangle|=1$. This equation is equivalent to (c.f. (6.11), (6.13)) $|\left\langle{\psi_{j}},{D_{j}(2\lambda g_{j}(1-{\rm e}^{{\rm i}t\omega_{j}t})/\omega_{j})\psi_{j}}\right\rangle|=1$, $\forall j\in\widetilde{J}$. $\blacksquare$

6.1.2 Entanglement entropy

Lemma 6.3** (Entanglement entropy)**

The entanglement entropy of $\rho_{\rm S}(t)$ and $\rho_{J}(t)$, given in Lemmas 6.1 and 6.2, respectively, are as follows.

[TABLE]

where $p$ is the (time-independent) population probability of the first level of the dimer and

[TABLE]

**Remarks. **

1. If the dimer is populated entirely on one level, i.e., $p=0$ or $p=1$, then $r_{J}(t)=1$ and $S(\rho_{J}(t))=0$ regardless of $J$ and $t$. That is, no entanglement entropy is created during the dynamics in any subsystem. This can be understood in the following way: for $\Psi_{\rm S}=\varphi_{1}$ the evolution of $\Psi(0)$, (6.2), under the dynamics generated by $H_{X}$, (1.2), is

[TABLE]

Therefore, the state is of product form for all times, so the reduction to any subsystem is again a pure state and has zero entropy.

2. The initially pure dimer state has maximal entropy for $p=1/2$. Then $r_{J}(t)=|s_{J}(t)|$.

*Proof of Lemma 6.3. * The entanglement entropy $S(\rho_{J}(t))=-{\rm Tr}\rho_{J}(t)\log\rho_{J}(t)$ of the state $\rho_{J}(t)$, (6.14), is expressed purely through the two non-zero eigenvalues of $\rho_{J}(t)$. (The operator $\rho_{J}(t)$ acts on an infinite-dimensional Hilbert space of $|\widetilde{J}|$ independent harmonic oscillators, however, all but two of its eigenvalues are zero.)

To find the two eigenvalues, we write $\rho_{J}=p|\Psi_{J}^{+}\rangle\langle\Psi_{J}^{+}|+(1-p)|\Psi_{J}^{-}\rangle\langle\Psi_{J}^{-}|$ (not exhibiting the $t$ dependence for simplicity of notation) and decompose

[TABLE]

Then we express $\rho_{J}$ in the ordered orthonormal basis $\{\Psi_{J}^{+},\hat{\eta}_{J}\}$, where $\hat{\eta}_{J}=\eta_{J}/\|\eta_{J}\|$. Namely, it has the form (2.6). The nonzero eigenvalues of $\rho_{J}$ are then immediately found,

[TABLE]

While $p$ is the time-independent population of the first level of the dimer, the overlap $s_{J}$ depends on time,

[TABLE]

Note that $|s_{J}(t)|\leq 1$ and $|s_{j}(t)|=1\Leftrightarrow$ $\rho_{J}(t)$ has rank one (c.f. Lemma 6.2). This completes the proof of Lemma 6.3. $\blacksquare$

6.2 Coherent state reservoir

Consider the initial reservoir state (1.1). Since $D(\xi)|\alpha\rangle={\rm e}^{{\rm i}{\rm Im}\xi\bar{\alpha}}|\alpha+\xi\rangle$ we obtain $\left\langle{\alpha_{j}},{D(\xi_{j})\alpha_{j}}\right\rangle={\rm e}^{{\rm i}{\rm Im}\xi_{j}\bar{\alpha}_{j}}\left\langle{\alpha_{j}},{\alpha_{j}+\xi_{j}}\right\rangle={\rm e}^{-2{\rm i}{\rm Im}\alpha_{j}\bar{\xi}_{j}}\ {\rm e}^{-\frac{1}{2}|\xi_{j}|^{2}}$, where $\xi_{j}=2(1-{\rm e}^{{\rm i}\omega_{j}t})\lambda g_{j}/\omega_{j}$. It follows that (see Lemma 6.3 and (6.19))

[TABLE]

We note that $|s_{J}(t)|$ is independent of the $\alpha_{1},\ldots,\alpha_{N}$ determining the initial reservoir state! Moreover, since the entanglement entropy only depends on this quantity, it too will not depend on the $\alpha_{j}$.

6.2.1 Proofs of Propositions 3.3 and 3.1

The results of Proposition 3.3 follow immediately from the expressions for the discrete oscillators by using (1.14). We now turn to the proof of Proposition 3.1.

In (2.7), we write $1-\cos(\omega t)=2\sin^{2}(\omega t/2)$ and use the fact that $\lim_{t\rightarrow\infty}\frac{1}{t}\int_{0}^{\infty}h_{X}(\omega)\frac{\sin^{2}(\omega t/2)}{\omega^{2}}d\omega=\pi h_{X}(0)/4$ to arrive at (3.1). This proves (A). The proof of (B) is obtained in the same way. To prove (C) we use the estimate

[TABLE]

Next

[TABLE]

Combining this with (6.21) yields the upper bound in (3.3). To get the lower bound in (3.3), we use $\int_{\omega_{0}}^{\omega_{1}}\frac{1-\cos(\omega t)}{\omega^{2}}d\omega\leq\frac{1}{\omega_{0}^{2}}(\omega_{1}-\omega_{0}+2/t)$ in (6.21). This completes the proof of Proposition 3.1. $\blacksquare$

7 $D$-interaction, proofs

7.1 Dimer dynamics for $V,\mu>0$, small

Theorem 7.1 below describes the change to this process under the influence of (small) energy-exchange interaction, i.e., for $V$ and $\mu$ in (1.3) small. The smallness of the energy-exchange terms is measured by the parameter

[TABLE]

We consider here the coupling function to satisfy

[TABLE]

If $\mu=V=0$ then $\eta=0$.

Theorem 7.1** (Dynamics of the dimer)**

Let $C_{\alpha}=\int_{0}^{\infty}h_{D}(\omega)f(\omega)d\omega$, where $h_{D}$ is given in (1.8) and $f$ is the coupling function in (1.3). The reduced dimer density matrix $\rho_{\rm S}(t)$ has the expansion

[TABLE]

Here, $\rho_{\rm S}^{0}(t)$ is given in (2.8) (with constant diagonal) and

[TABLE]

In (7.3), $C$ is a constant not depending on any of the parameters $V,\Omega,\eta,t,\mu$, nor on the functions $h$, $\alpha$.

The remainder in (7.3) is not uniform in time. However, up to times $t\sim t_{\rm pd}$, when the partial decoherence is completed (c.f. (2.11)), the remainder is at most $O(B_{1})$, where

[TABLE]

(use $t=\sqrt{|\xi|}/\lambda$ in the right side of (7.3)). The expansion (7.3) is meaningful provided the remainder is smaller than the main term. Therefore, in order to resolve the effect of the term $\propto\frac{V}{\Omega}$ on the dimer dynamics (the left side of (7.3)) on the time-scale $t_{\rm pd}$, we need $B_{1}<\!\!<\frac{V}{\Omega}$. Using the definition of $\eta$, (7.1) (parameter regime of the second line), we see that

[TABLE]

provided that

[TABLE]

This shows the following result.

Corollary 7.2** (Validity of the perturbation expansion)**

In the parameter regime (7.6), the right side of (7.3) is $<\!\!<V/\Omega$. This means that the perturbation expansion (7.3) resolves the partial decoherence process and the $O(V/\Omega)$ corrections to the dynamics, and in particular, the associated change in the dimer entanglement entropy.

**Remarks. ** 1. If $g^{\prime}(0)\neq 0$ then

[TABLE]

2. If $\varepsilon>0$, $\delta\geq 0$ and $g(\omega)=\varepsilon+\omega^{\delta}$ as $\omega\sim 0$, then we obtain

[TABLE]

Note that strictly speaking, we need $\varepsilon>0$ to have $\inf_{\omega}g(\omega)>0$, however, the value of $\xi$ does not depend on $\varepsilon$. This independence can be understood as follows: shifting the coupling function $g(\omega)$ in (1.3) by a constant value $\varepsilon$ amounts to renormalizing the frequency $\Omega$ as $\Omega+2\lambda\varepsilon\sum_{j=1}^{N}a^{\dagger}_{j}a_{j}$. However, the energy-conserving partial decoherence process does not depend on that frequency.

3. If one is only interested in the closeness of $\rho_{\rm S}(t)$ to $\rho^{0}_{\rm S}(t)$ then the term $O(V/\Omega)$ can be included in the remainder, and the remainder stays small for a wider parameter regime than (7.6).

4. We derive a more precise estimate than (7.3). Namely, we can omit the remainder term $tV^{3}/\Omega^{2}$ on the right side of (7.3) by replacing $\rho^{1}_{\rm S}(t)$ in the left side of (7.3) by

[TABLE]

Due to the bound ${\rm e}^{-it\sqrt{\Omega^{2}+V^{2}}}={\rm e}^{-i\Omega t}(1+O(tV^{2}/\Omega))$, we have $\widetilde{\rho}_{\rm S}^{1}(t)=\rho_{\rm S}^{1}(t)+O(tV^{2}/\Omega)$.

7.2 Dynamics of the oscillators for $V>0$ small, $\mu=0$

For simplicity of the presentation, we consider in this section the Hamiltonian (1.3) with $\mu=0$ only. Let $J\subseteq[0,\infty)$ be a window of continuous oscillator frequencies of interest.

Theorem 7.3** (Dynamics of the oscillators)**

Consider the system for finite $N$, with Hamiltonian (1.3) and $\mu=0$. Denote by $\rho_{J}(t)$ the reduced density matrix of the oscillators having (discrete) frequencies inside $J\subseteq\mathbb{R}$. We have

[TABLE]

where $C$ is a constant independent of $N,J,t,V,\Omega,\lambda$. Here, $m$ is given in (7.2) and $\rho^{00}_{J}$ and $\rho^{1}_{J}$ are as in Theorem 2.3, (B). On the left side of (7.10), $\|\cdot\|_{1}$ denotes the trace norm of operators acting on the Hilbert space of $K$ harmonic oscillators, where $K$ is the number of discrete frequencies $\omega_{j}$ lying inside $J$.

7.3 Proof of Theorem 7.1

The reduced dimer density matrix is

[TABLE]

where

[TABLE]

and

[TABLE]

The spectral representation of ${\rm e}^{-{\rm i}t\bar{H}}$ is

[TABLE]

where

[TABLE]

and, with ${\cal N}=2(\bar{\omega}^{2}+\bar{V}^{2}+\bar{\omega}\sqrt{\bar{V}^{2}+\bar{\omega}^{2}})$,

[TABLE]

In the sum in (7.11) we split off the term with $n_{1}=\ldots=n_{N}=0$,

[TABLE]

where $F(0)=\prod_{j=1}^{N}{\rm e}^{-|\alpha_{j}|^{2}}$,

[TABLE]

and where the $*$ indicates that we omit the term in which all $n_{1},\ldots,n_{N}$ vanish in the summation. The point of this splitting is that in the remaining sum at least one of the $n_{1},\ldots,n_{N}$ is $\geq 1$, so that we can estimate the parameter

[TABLE]

by (see also (7.1))

[TABLE]

Then we use the decomposition (7.14) and the expansions

[TABLE]

and

[TABLE]

to estimate

[TABLE]

It follows that

[TABLE]

where the remainder has the estimate (c.f. (7.12))

[TABLE]

and where $C$ is a constant independent of $t,\eta,V,\mu,\alpha_{j},N$. Finally we ‘remove the *’ in the sum on the right side of (7.3) by adding (and subtracting) the summand for $n_{1}=\ldots=n_{N}=0$, and we obtain

[TABLE]

Here we have set

[TABLE]

and $\psi_{\rm S}=a|\!\uparrow\rangle+b|\!\downarrow\rangle$, $p=|a|^{2}$. Combining (7.17) and (7.3) yields

[TABLE]

Note that the main term has trace one and, for $t=0$, it reduces to the correct value $\rho_{\rm S}(0)$. We may now diagonalize $H^{0}$ given in (7.18) exactly analogously to (7.14)-(7.16). However, we resolve the $O(V/\Omega)$-term in the expansion, namely,

[TABLE]

Using these expressions in the spectral representation ${\rm e}^{-{\rm i}tH^{0}}={\rm e}^{-{\rm i}tr^{0}}|\chi_{+}\rangle\langle\chi_{+}|+{\rm e}^{{\rm i}tr^{0}}|\chi_{-}\rangle\langle\chi_{-}|$ with $r_{0}=\frac{1}{2}\sqrt{\Omega^{2}+V^{2}}$, an easy calculation leads to

[TABLE]

Combining (7.3) with (7.32) and taking the continuous mode limit shows (7.3). This finishes the proof of Theorem 7.1. $\blacksquare$

7.4 Proofs of Theorem 7.3 and Proposition 3.3

We examine the reduced oscillator density matrix $\rho_{r}$ for oscillators with indices $j=1,\ldots,r$, where $1\leq r\leq N$.

Let $H$ be the Hamiltonian (1.3) with $\mu=0$ and let $|M_{1},\ldots,M_{r}\rangle$ and $|\bar{M}_{1},\ldots\bar{M}_{r}\rangle$ be number states of the first $r$ oscillators ($M_{j},\bar{M}_{j}\in\mathbb{N}$). The matrix elements of $\rho_{r}(t)$ are

[TABLE]

where we recall that the initial state is (1.1) and where we have set

[TABLE]

and

[TABLE]

Of course, $H_{\bar{M}}$ is (7.35) with $M_{j}$ replace with $\bar{M}_{j}$. It is readily seen that for $V=0$, the scalar product factor in the sum of (7.33) is simply $p{\rm e}^{-{\rm i}\lambda t\sum_{j=1}^{r}g_{j}(M_{j}-\bar{M}_{j})}+(1-p){\rm e}^{{\rm i}\lambda t\sum_{j=1}^{r}g_{j}(M_{j}-\bar{M}_{j})}$ (where $p=|\left\langle{\psi_{\rm S}},{\uparrow}\right\rangle|^{2}$). A little algebra on (7.33) then gives the following evolution of $\rho_{r}$, for $V=0$:

[TABLE]

with $\rho_{r}(0)=|\alpha_{1},\ldots,\alpha_{r}\rangle\langle\alpha_{1},\ldots,\alpha_{r}|$ and

[TABLE]

In particular, $\rho_{r}^{(V=0)}(t)$ acts non-trivially on ${\rm span}\{U_{t}^{+}|\alpha_{1},\ldots,\alpha_{r}\rangle,U_{t}^{-}|\alpha_{1},\ldots,\alpha_{r}\rangle\}$ and has rank two ($t\neq 0$). By the Gram-Schmidt procedure we construct an orthonormal basis of this span and express $\rho_{r}^{(V=0)}(t)$ as a $2\times 2$ matrix (see (7.54) below). To deal with the situation $V\neq 0$, we diagonalize the matrices $H_{M}$ and $H_{\bar{M}}$ in (7.33),

[TABLE]

We may use the formulas (7.14)-(7.16) with suitable substitutions. In particular,

[TABLE]

with a remainder uniform in the values of $M_{1},\ldots,M_{N}$, and

[TABLE]

Our strategy is now similar to what we did in the derivation of the dimer density matrix. We split off the term in (7.33) in which $n_{r+1}=\cdots=n_{N}=0$, which we can calculate explicitly (to its main contribution). The remaining terms in the sum have at least one $n_{j}\geq 1$ for some $j\geq r+1$, and this will make those terms small, provided that $\min_{j\geq r+1}g_{j}$ is relatively large.

Using formula (7.38) for each exponential and expanding the eigenvectors as in (7.39), it is not hard (albeit a bit tedious) to obtain the following bound,

[TABLE]

Here, we have set $W_{M}=2\sum_{j=1}^{r}g_{j}M_{j}$. Noticing that

[TABLE]

we see that the term with $n_{r+1}=\cdots=n_{N}=0$ on the right side of (7.33) is the $\langle M_{1},\ldots,M_{r}|\cdot|\bar{M}_{1},\ldots,\bar{M}_{r}\rangle$-matrix element of the selfadjoint operator

[TABLE]

So far, the remainder in (7.43) is estimated in the maximum norm for matrices. More precisely, let $\rho_{r,0}(t)$ be the operator obtained by taking only the term $n_{r+1},\ldots,n_{N}=0$ in (7.33), and let $\mu(t)$ be the operator (7.43), without the error term $O$. What we have shown so far is that

[TABLE]

where $C$ is a constant independent of time $t$ and of $r$. We show now that we actually have the stronger estimate $\|\rho_{r,0}-\mu(t)\|_{1}\leq C(1+t)V^{2}/\Omega$ for the same constant $C$, where $\|A\|_{1}={\rm Tr}|A|$ is the trace norm, with the absolute value of an operator $A$ defined by $|A|=\sqrt{A^{*}A}$. To do so, we recall that the error contains terms brought about by replacing ${\rm e}^{{\rm i}tH_{\bar{M}}}$ and ${\rm e}^{-{\rm i}tH_{M}}$ in the term $\left\langle{\psi_{\rm S}},{{\rm e}^{{\rm i}tH_{\bar{M}}}{\rm e}^{-{\rm i}tH_{M}}\psi_{\rm S}}\right\rangle$ in (7.33) by using (7.39) and (7.40). It is then not difficult to see that this remainder is a sum of at most 16 terms, each of the form $A|\alpha\rangle\langle\alpha|B$, for some operators $A$ and $B$ satisfying $\|A\|\|B\|\leq C(1+t)V^{2}/\Omega$, for a $C$ independent of $t,r$. To estimate this remainder in trace norm, we note that

[TABLE]

where $|U\rangle$ is the normalization of $B^{*}|\alpha\rangle$. Hence $|(A|\alpha\rangle\langle\alpha|B)|$ is rank-one and thus its trace, being the sum of its eigenvalues, equals $\|A|\alpha\rangle\|\,\|B^{*}|\alpha\rangle\|\leq\|A\|\,\|B^{*}\|\leq C(1+t)V^{2}/\Omega$. We conclude that

[TABLE]

In other words, the remainder term in (7.43) is understood in the trace-norm topology.

Let

[TABLE]

be the orthogonal projection onto the state $|n_{1}=0,\ldots,n_{r}=0\rangle$. We have $U_{t}^{\pm}P_{0}=P_{0}$, $\hat{W}P_{0}=0$ and $P^{\perp}_{0}\hat{W}\geq mP^{\perp}_{0}$, where $m=2\min_{j\geq 1}g_{j}$. (Here, $P^{\perp}_{0}=\mathchoice{\rm 1\mskip-4.0mul}{\rm 1\mskip-4.0mul}{\rm 1\mskip-4.5mul}{\rm 1\mskip-5.0mul}-P_{0}$.) Therefore,

[TABLE]

Using this bound in (7.43), we arrive at

[TABLE]

The error $\tfrac{V}{\Omega+\lambda m}$ will be small compared to $V/\Omega$ for sizable $m$. Again, this error is estimated in the trace-norm topology (as in (7.45)).

Now we treat the terms in the sum (7.33) in which at least one $n_{j}$ is $\geq 1$. We expand $\left\langle{\psi_{\rm S}},{{\rm e}^{{\rm i}tH_{\bar{M}}}{\rm e}^{-{\rm i}tH_{M}}\psi_{\rm S}}\right\rangle$ using (7.38)-(7.40), as we did to arrive at (7.41). We allow for errors of the order $V/(\Omega+\lambda T_{M})\leq V/(\Omega+\lambda m)$, where $m=2\min_{1\leq j\leq N}g_{j}$. For instance, we use (c.f. (7.39)) $\chi_{+}^{M}=\textstyle\begin{pmatrix}1\\ 0\end{pmatrix}+O(V/(\Omega+\lambda m))$. We then obtain, whenever at least one $n_{j}\geq 1$ for some $j=r+1,\ldots,N$,

[TABLE]

The right side of (7.48) is independent of $n_{j}$ with $j\geq r+1$. Moreover,

[TABLE]

where the star ∗ indicates that we sum over all terms except the one with all $n_{j}$ equal to zero. Combining (7.48), (7.49) and the definition of $F$, (7.34), we see that the contribution to the sum (7.33) coming from all terms where at least one $n_{j}\geq 1$, is the $\langle M_{1},\ldots,M_{N}|\cdot|\bar{M}_{1},\ldots,\bar{M}_{N}\rangle$-matrix element of the operator

[TABLE]

As above, the error is estimated in the trace-norm topology. The whole sum in (7.33) is thus the $\langle M_{1},\ldots,M_{N}|\cdot|\bar{M}_{1},\ldots,\bar{M}_{N}\rangle$-matrix element of the sum of the two operators (7.47) plus (7.50). We conclude that

[TABLE]

(error estimated in trace-norm topology) where we have defined

[TABLE]

This shows that, modulo the remainder in (7.51),

[TABLE]

where $|f_{+}\rangle=(\tfrac{1}{2}pU_{t}^{+}+\bar{v}P_{0})|\alpha_{1},\ldots,\alpha_{r}\rangle$, $|f_{-}\rangle=(\tfrac{1}{2}(1-p)U_{t}^{-}-vP_{0})|\alpha_{1},\ldots,\alpha_{r}\rangle$ and

[TABLE]

Since $P_{0}|\alpha_{1},\ldots,\alpha_{r}\rangle\propto|0\rangle\equiv|0,\ldots,0\rangle$ this shows that (the main part of) $\rho_{r}(t)$ acts nontrivially on the three-dimensional space spanned by the vectors $|e_{\pm}\rangle$ and $|0\rangle\equiv|0,\ldots,0\rangle$. By the Gram-Schmidt procedure we construct an orthonormal basis of ${\rm span}\{|e_{+}\rangle,|e_{-}\rangle,|0\rangle\}$ and then express $\rho_{r}(t)$ as a $3\times 3$ matrix. The orthonormal basis is given by $\{|e_{+}\rangle,|\hat{\eta}\rangle,|\hat{\chi}\rangle\}$, where

[TABLE]

Here we have set

[TABLE]

It is now a straightforward calculation to express $\rho_{t}(t)$, (7.53), as a matrix which in the continuous mode limit is precisely (2.17). Note that for simplicity of notation, we have derived the reduced density matrix of the oscillators with indices $j=1,\ldots,r$ in this proof; the general case resulting in the reduced density matrix of oscillators with frequencies in a window $J\subseteq\mathbb{R}$ is directly obtained from the expressions derived here, just as we did it in Lemma 6.2. This completes the proof of Theorem 7.3. $\blacksquare$

**Proof of Proposition 3.3. ** We write the $3\times 3$ matrix associated to the operator $\rho_{J}^{0}(t)+\frac{V}{\Omega}\rho_{J}^{1}(t)$ (c.f. Theorem 7.3) as $M_{0}+\frac{V}{\Omega}M_{1}$ and solve the equation $\det(M_{0}+\frac{V}{\Omega}M_{1}-z)=0$ for $z$. For $V=0$ we obtain the three values $z_{0}=0,\frac{1}{2}\pm\frac{1}{2}r_{J}$, where $r_{J}$ is given in (1.12). Then, for each $z_{0}$, we make the Ansatz $z=z_{0}+\frac{V}{\Omega}t+\frac{V^{2}}{\Omega^{2}}t^{\prime}+\cdots$ and solve the equation $\det(M_{0}+\frac{V}{\Omega}M_{1}-z)=0$ for $t$. We obtain one value for each $z_{0}$. $\blacksquare$

Acknowledgment

This work was carried out under the auspices of the National Nuclear Security Administration of the U.S. Department of Energy at Los Alamos National Laboratory under Contract No. DE-AC52-06NA25396. M.M. has been supported by NSERC through a Discovery Grant and a Discovery Accelerator Supplement. M.M. is grateful for the hospitality and financial support of LANL, where part of this work was carried out. A.I.N. acknowledges support from the CONACyT, Grant No. 15349 and partial support during his visit from the Biology Division, B-11, at LANL. G.P.B. and R.T.S. acknowledge support from the LDRD program at LANL.