Relaxation to Gaussian Generalized Gibbs Ensembles in Quadratic Bosonic Systems in the Thermodynamic Limit

TL;DR

This paper demonstrates that quadratic bosonic systems in the thermodynamic limit relax to Gaussian generalized Gibbs ensembles, revealing a general mechanism for equilibration based on diagonal singularity.

Contribution

It introduces a universal mechanism explaining relaxation to GGEs in quadratic bosonic systems, including coupled systems, under certain physical conditions.

Findings

Quadratic bosonic systems relax to Gaussian GGEs.

Relaxation occurs from general initial states, not necessarily Gaussian.

Explicit analytical demonstration for coupled harmonic oscillator systems.

Abstract

Integrable quantum many-body systems are considered to equilibrate to generalized Gibbs ensembles (GGEs) characterized by the expectation values of integrals of motion. We study the dynamics of exactly solvable quadratic bosonic systems in the thermodynamic limit, and show a general mechanism for the relaxation to GGEs, in terms of the diagonal singularity. We show analytically and explicitly that a free bosonic system relaxes from a general (not necessarily Gaussian) initial state under certain physical conditions to a Gaussian GGE. We also show the relaxation to a Gaussian GGE in an exactly solvable coupled system, a harmonic oscillator linearly interacting with bosonic reservoirs.

Equations

$${\hat{\rho }}_{\text{GGE}}\propto \exp \left(-\int {\text{d}}^D\mathbf{k}\ln (1+f_{\mathbf{k}}^{-1}){\hat{b}}_{\mathbf{k}}^{†}{\hat{b}}_{\mathbf{k}}\right)$$

$$\hat{H}=\int {\text{d}}^D\mathbf{k}{\omega }_k{\hat{b}}_{\mathbf{k}}^{†}{\hat{b}}_{\mathbf{k}},$$

$${\hat{\rho }}_0\propto \exp \left[-\frac{1}{2}\int {\text{d}}^D\mathbf{k}\int {\text{d}}^D{\mathbf{k}}^{\prime }\left(\begin{array}{cc}{\hat{b}}_{\mathbf{k}}^{†}-⟨{\hat{b}}_{\mathbf{k}}^{†}⟩& {\hat{b}}_{\mathbf{k}}-⟨{\hat{b}}_{\mathbf{k}}⟩\end{array}\right){\Theta }_{\mathbf{k}{\mathbf{k}}^{\prime }}\left(\begin{array}{c}\text{vskip}6.0ptplus2.0ptminus2.0pt{\hat{b}}_{{\mathbf{k}}^{\prime }}-⟨{\hat{b}}_{{\mathbf{k}}^{\prime }}⟩\\ {\hat{b}}_{{\mathbf{k}}^{\prime }}^{†}-⟨{\hat{b}}_{{\mathbf{k}}^{\prime }}^{†}⟩\end{array}\right)\right],$$

$${\chi }_0[J,J^{\ast }]$$

$$=\exp \left[-\frac{1}{2}\int {\text{d}}^D\mathbf{k}\int {\text{d}}^D{\mathbf{k}}^{\prime }\left(\begin{array}{cc}{J}_{\mathbf{k}}^{\ast }& J_{\mathbf{k}}\end{array}\right)V_{\mathbf{k}{\mathbf{k}}^{\prime }}\left(\begin{array}{c}\text{vskip}6.0ptplus2.0ptminus2.0pt{J}_{{\mathbf{k}}^{\prime }}\\ J_{{\mathbf{k}}^{\prime }}^{\ast }\end{array}\right)+\int {\text{d}}^D\mathbf{k}(J_{\mathbf{k}}⟨{\hat{b}}_{\mathbf{k}}^{†}⟩-J_{\mathbf{k}}^{\ast }⟨{\hat{b}}_{\mathbf{k}}⟩)\right],$$

$$V_{\mathbf{k}{\mathbf{k}}^{\prime }}=\left(\begin{array}{cc}\text{vskip}6.0ptplus2.0ptminus2.0pt\frac{1}{2}⟨\{{\hat{b}}_{\mathbf{k}}-⟨{\hat{b}}_{\mathbf{k}}⟩,{\hat{b}}_{{\mathbf{k}}^{\prime }}^{†}-⟨{\hat{b}}_{{\mathbf{k}}^{\prime }}^{†}⟩\}⟩& -⟨({\hat{b}}_{\mathbf{k}}-⟨{\hat{b}}_{\mathbf{k}}⟩)({\hat{b}}_{{\mathbf{k}}^{\prime }}-⟨{\hat{b}}_{{\mathbf{k}}^{\prime }}⟩)⟩\\ -⟨({\hat{b}}_{\mathbf{k}}^{†}-⟨{\hat{b}}_{\mathbf{k}}^{†}⟩)({\hat{b}}_{{\mathbf{k}}^{\prime }}^{†}-⟨{\hat{b}}_{{\mathbf{k}}^{\prime }}^{†}⟩)⟩& \frac{1}{2}⟨\{{\hat{b}}_{\mathbf{k}}^{†}-⟨{\hat{b}}_{\mathbf{k}}^{†}⟩,{\hat{b}}_{{\mathbf{k}}^{\prime }}-⟨{\hat{b}}_{{\mathbf{k}}^{\prime }}⟩\}⟩\end{array}\right)$$

$$\Theta =2Z{\coth }^{-1}(2ZV)=Z\ln \frac{1+(2ZV{)}^{-1}}{1-(2ZV{)}^{-1}}$$

$$⟨({\hat{b}}_{\mathbf{k}}^{†}-⟨{\hat{b}}_{\mathbf{k}}^{†}⟩)({\hat{b}}_{{\mathbf{k}}^{\prime }}-⟨{\hat{b}}_{{\mathbf{k}}^{\prime }}⟩)⟩$$

$$⟨({\hat{b}}_{\mathbf{k}}-⟨{\hat{b}}_{\mathbf{k}}⟩)({\hat{b}}_{{\mathbf{k}}^{\prime }}-⟨{\hat{b}}_{{\mathbf{k}}^{\prime }}⟩)⟩$$

$$V_{\mathbf{k}{\mathbf{k}}^{\prime }}=\left(\begin{array}{cc}\text{vskip}6.0ptplus2.0ptminus2.0pt(\frac{1}{2}+f_{\mathbf{k}}){\delta }^D(\mathbf{k}-{\mathbf{k}}^{\prime })+F_{\mathbf{k}{\mathbf{k}}^{\prime }}^{\ast }& -g_{\mathbf{k}}{\delta }^D(\mathbf{k}+{\mathbf{k}}^{\prime })-G_{\mathbf{k}{\mathbf{k}}^{\prime }}\\ -g_{\mathbf{k}}^{\ast }{\delta }^D(\mathbf{k}+{\mathbf{k}}^{\prime })-G_{\mathbf{k}{\mathbf{k}}^{\prime }}^{\ast }& (\frac{1}{2}+f_{\mathbf{k}}){\delta }^D(\mathbf{k}-{\mathbf{k}}^{\prime })+F_{\mathbf{k}{\mathbf{k}}^{\prime }}\end{array}\right).$$

$${\chi }_t[J,J^{\ast }]$$

$$=⟨{\text{e}}^{\int {\text{d}}^D\mathbf{k}(J_{\mathbf{k}}{\hat{b}}_{\mathbf{k}}^{†}{\text{e}}^{\text{i}{\omega }_{k}t}-J_{\mathbf{k}}^{\ast }{\hat{b}}_{\mathbf{k}}{\text{e}}^{-\text{i}{\omega }_{k}t})}⟩$$

$$={\chi }_0[J{\text{e}}^{\text{i}\omega t},J^{\ast }{\text{e}}^{-\text{i}\omega t}],$$

$$\displaystyle\chi_{t}[J,J^{*}]=\exp\biggl{(}$$

$$\displaystyle{}-\int\text{d}^{D}\bm{k}\int\text{d}^{D}\bm{k}^{\prime}\,[J_{\bm{k}}F_{\bm{k}\bm{k}^{\prime}}J_{\bm{k}^{\prime}}^{*}\text{e}^{\text{i}(\omega_{k}-\omega_{k^{\prime}})t}-\mathop{\text{Re}}\nolimits(J_{\bm{k}}G_{\bm{k}\bm{k}^{\prime}}^{*}J_{\bm{k}^{\prime}}\text{e}^{\text{i}(\omega_{k}+\omega_{k^{\prime}})t})]+2\text{i}\mathop{\text{Im}}\nolimits\int\text{d}^{D}\bm{k}\,J_{\bm{k}}\text{e}^{\text{i}\omega_{k}t}\langle\hat{b}_{\bm{k}}^{\dagger}\rangle\biggr{)}.$$

$${\chi }_t[J,J^{\ast }]$$

$$={\chi }_{\text{GGE}}[J,J^{\ast }],$$

$$\hat{\rho }(t)\stackrel{t\to \infty }{\to }{\hat{\rho }}_{\text{GGE}}\propto \exp \left(-\int {\text{d}}^D\mathbf{k}\ln (1+f_{\mathbf{k}}^{-1}){\hat{b}}_{\mathbf{k}}^{†}{\hat{b}}_{\mathbf{k}}\right),$$

$$\ln {\chi }_0[J,J^{\ast }]$$

$$=\cdots +\int {\text{d}}^D{\mathbf{k}}_1\int {\text{d}}^D{\mathbf{k}}_2\int {\text{d}}^D{\mathbf{k}}_3{K}_{{\mathbf{k}}_1{\mathbf{k}}_2{\mathbf{k}}_3}{J}_{{\mathbf{k}}_1}{J}_{{\mathbf{k}}_2}{J}_{{\mathbf{k}}_3}^{\ast }$$

$$+\cdots$$

$$K_{{\mathbf{k}}_1{\mathbf{k}}_2{\mathbf{k}}_3}={\bar{K}}_{{\mathbf{k}}_1{\mathbf{k}}_2}{\delta }^D({\mathbf{k}}_1+{\mathbf{k}}_2-{\mathbf{k}}_3)+{\tilde{K}}_{{\mathbf{k}}_1{\mathbf{k}}_2{\mathbf{k}}_3},$$

$$\ln {\chi }_t[J,J^{\ast }]$$

$$=\cdots +\int {\text{d}}^D{\mathbf{k}}_1\int {\text{d}}^D{\mathbf{k}}_2\int {\text{d}}^D{\mathbf{k}}_3{K}_{{\mathbf{k}}_1{\mathbf{k}}_2{\mathbf{k}}_3}{J}_{{\mathbf{k}}_1}{J}_{{\mathbf{k}}_2}{J}_{{\mathbf{k}}_3}^{\ast }$$

$$\times {\text{e}}^{\text{i}({\omega }_{k_1}+{\omega }_{k_2}-{\omega }_{k_3})t}$$

$$+\cdots .$$

$$⟨\hat{A}\hat{B}(t)\hat{C}{⟩}_{\text{GGE}}\stackrel{t\to \infty }{\to }⟨\hat{A}\hat{C}{⟩}_{\text{GGE}}⟨\hat{B}{⟩}_{\text{GGE}}$$

$${\hat{\rho }}_0$$

$$\mapsto \hat{\rho }(t)={\text{e}}^{-\text{i}\hat{H}t}{\hat{\rho }}_0{\text{e}}^{\text{i}\hat{H}t}\stackrel{t\to \infty }{\to }{\hat{\rho }}_{\text{GGE}},$$

$$⟨\hat{A}{⟩}_t$$

$$=\sum _j⟨{\hat{L}}_j^{†}\hat{A}(t){\hat{L}}_j{⟩}_{\text{GGE}}$$

$$\stackrel{t\to \infty }{\to }\sum _j⟨{\hat{L}}_j^{†}{\hat{L}}_j{⟩}_{\text{GGE}}⟨\hat{A}{⟩}_{\text{GGE}}$$

$$=⟨\hat{A}{⟩}_{\text{GGE}}$$

$$\hat{H}={\hat{H}}_S+\sum _{\ell }{\hat{H}}_{\ell }+\lambda \hat{V}$$

$${\hat{H}}_S=\Omega {\hat{a}}^{†}\hat{a},{\hat{H}}_{\ell }=\int {\text{d}}^D\mathbf{k}{\omega }_{k\ell }{\hat{b}}_{\mathbf{k}\ell }^{†}{\hat{b}}_{\mathbf{k}\ell },$$

$$\hat{V}=\sum _{\ell }\int {\text{d}}^D\mathbf{k}(u_{\mathbf{k}\ell }^{\ast }{\hat{a}}^{†}{\hat{b}}_{\mathbf{k}\ell }+u_{\mathbf{k}\ell }{\hat{b}}_{\mathbf{k}\ell }^{†}\hat{a}).$$

$$[\hat{a},\hat{a}]=0,[\hat{a},{\hat{a}}^{†}]=1,$$

$$[{\hat{b}}_{\mathbf{k}\ell },{\hat{b}}_{{\mathbf{k}}^{\prime }{\ell }^{\prime }}]=0,[{\hat{b}}_{\mathbf{k}\ell },{\hat{b}}_{{\mathbf{k}}^{\prime }{\ell }^{\prime }}^{†}]={\delta }_{\ell {\ell }^{\prime }}{\delta }^D(\mathbf{k}-{\mathbf{k}}^{\prime }),$$

$$[\hat{a},{\hat{b}}_{\mathbf{k}\ell }]=[\hat{a},b_{\mathbf{k}\ell }^{†}]=0.$$

Full text

I Introduction

Recent advances in experimental studies of relaxation processes with quantum atomic gases Kinoshita1 ; Kinoshita2 ; Gring1 ; Schmiedmayer1 ; Langen1 ; Kaufman1 have been stimulating theoretical studies on equilibration and prethermalization of isolated quantum many-body systems ref:PolkovnikovRMP ; Yukalov1 ; ref:EisertFriesdorfGogolin-ReviewNatPhys ; ref:GooldReviewJPA ; Goglin1 ; Vidmar1 ; ref:CalabresePhysicaA ; ref:MoriIkedaKaminishiUeda-ReviewJPB . Among various issues, it is recognized that an integrable large quantum system would relax to a nonthermal steady state described by a generalized Gibbs ensemble (GGE) Rigol2 ; ref:CramersDawsonEisertOsbornePRL2008 ; Rigol1 ; Caux1 ; Wouters1 , which is characterized by the expectation values of a set of integrals of motion.

In Ref. Rigol2 , the GGE was conjectured on the basis of the principle of maximum entropy Jaynes1 and verified numerically, and then, in Refs. Caux1 ; Wouters1 , it was studied with integrable systems under initial quench scenarios. There exist many studies on this issue: see reviews ref:PolkovnikovRMP ; ref:EisertFriesdorfGogolin-ReviewNatPhys ; ref:GooldReviewJPA ; Goglin1 ; Vidmar1 ; ref:CalabresePhysicaA ; ref:MoriIkedaKaminishiUeda-ReviewJPB and references therein. Nonetheless, it remains a challenging open problem to clarify a concrete scenario of the microscopic mechanism of the relaxation to GGE, e.g. in the context of the Liouville integrability in Hamiltonian systems, where we essentially perform canonical transformations to an assembly of noninteracting oscillators, i.e. action-angle variables ref:KuboTextbookI-Chap5 . For the unitary evolution of nonintegrable systems, on the other hand, there are intensive studies, e.g. to give foundations to the quantum ergodic theorem vonNeumann1 ; vonNeumann2 ; Typicality1 ; Goldstein2 ; Typicality3 ; Typicality2 and to understand the relaxation times relaxation1 ; Goldstein3 ; Monnai2 ; relaxation2 ; Goldstein1b ; relaxation3 .

In Ref. ref:CramerEisert-CentralLimitTheorem (see also Ref. ref:CramersDawsonEisertOsbornePRL2008 ), it was shown that a one-dimensional bosonic lattice system with a quadratic Hamiltonian with nearest-neighbor hopping terms locally equilibrates from an arbitrary initial state fulfilling conditions such as clustering and the absence of anomalous correlations to a Gaussian GGE for a finite time duration. Similarly, in Ref. ref:GogolinEisert-GaussificationFermi , the local equilibration to a Gaussian GGE during a finite time interval was shown for a fermionic $d_{\mathcal{L}}$-dimensional cubic lattice system with a quadratic Hamiltonian with finite-range interactions, under a clustering condition on the initial state. The key to the proofs is the Lieb-Robinson bound Goglin1 . In addition, power-law relaxation towards Gaussian GGE was studied for similar bosonic/fermionic lattice models with quadratic Hamiltonians, by making use of the Kusmin-Landau bound Farrelly1 and the stationary-phase approximation Murthy1 under the clustering condition on the initial state. See also Ref. Sotiriadis1 for further arguments on the clustering decomposition in the context of Gaussian GGE. These results are interesting because it is shown that the systems equilibrate to the GGEs during time intervals (not simply on average in time Goglin1 ) and that only subsets of the sets of conserved quantities of the integrable systems are relevant to the GGEs.

In this paper, we provide another contribution to this issue, studying exactly solvable quadratic bosonic models, namely, a class of models which can be mapped to free bosonic fields in the thermodynamic limit. The existence of such a mapping corresponds to the Liouville integrability for classical systems note:QuantumIntegrability . We solve the evolution of the state of the whole system in the thermodynamic limit exactly and explicitly, and observe that it relaxes to a simple Gaussian GGE in the long-time limit, under certain conditions on the initial state.

More specifically, we are going to show the following. We consider solvable bosonic systems, with canonical normal modes say $\hat{b}_{\bm{k}}$, in the thermodynamic limit in $D$-dimensional space. We assume the following physical conditions on the initial state.

1. (i)

The initial state of the system, which is non-Gaussian in general, is prepared irrespective of the Hamiltonian of the system (like in a quench scenario). 2. (ii)

We allow the correlations in the initial state to possess translationally invariant components (particles can be distributed all over the space and correlations can exist everywhere in space in the initial state). 3. (iii)

But the correlations in the initial state are assumed to be of finite range. 4. (iv)

The first moment $\langle\hat{b}_{\bm{k}}\rangle$ in the initial state is assumed to be nonvanishing only locally, i.e. free from a translationally invariant component.

In addition:

1. (v)

We exclude observables extending all over the space, since the expectation values of such quantities diverge in the thermodynamic limit and are not actually measurable.

Under these conditions, we show that the system relaxes to a Gaussian GGE of the form

[TABLE]

in the long-time limit, where $f_{\bm{k}}$ is the translationally invariant component of the single-particle correlation $\langle(\hat{b}_{\bm{k}}^{\dagger}-\langle\hat{b}_{\bm{k}}^{\dagger}\rangle)(\hat{b}_{\bm{k}^{\prime}}-\langle\hat{b}_{\bm{k}^{\prime}}\rangle)\rangle$ in the initial state. Observe that the GGE in (1) admits only the occupation numbers $\hat{I}_{\bm{k}}=\hat{b}_{\bm{k}}^{\dagger}\hat{b}_{\bm{k}}$ as integrals of motion (cf. Sels1 ; Polkovnikov1 ). The state $\hat{\rho}_{\text{GGE}}$ in (1) is a Gaussian state, since its Wigner function is Gaussian ref:ContVarQI ; ref:GaussianQI . Therefore, we call the state $\hat{\rho}_{\text{GGE}}$ in (1) Gaussian GGE. The correlations among the occupation numbers such as $I_{\bm{k}}I_{\bm{k}^{\prime}}$ do not survive in the equilibrium state $\hat{\rho}_{\text{GGE}}$.

In contrast to the previous works ref:CramerEisert-CentralLimitTheorem and ref:GogolinEisert-GaussificationFermi , in which the results are rigorously proved for large but finite lattice systems, we directly go to the thermodynamic limit and deal with field-theoretical Hamiltonians with continuous spectrum. This greatly simplifies the analysis and the picture. Our results are valid for any spatial dimension, anomalous correlations are allowed in the initial states, and we do not focus on a part of the system (we do not take partial trace) Note1 . In this approach, the relaxation to a GGE is understood in terms of the diagonal singularity VANHOVE1955901 ; ANTONIOU1997737 , or equivalently, as a consequence of the Riemann-Lebesgue lemma ref:FourierAnalysisKoerner , smearing the spectrum due to the rapid oscillations in the long-time limit. We will see that some “regularity” of the correlation functions and the observables is important, for the Riemann-Lebesgue lemma to work. In particular, we stress that observables are relevant to the relaxation process (e.g. for the relaxation time).

We first study the dynamics of a free bosonic field in Sec. II. We show the relaxation from Gaussian initial states to the Gaussian GGE in (1), and generalize the result to non-Gaussian initial states. We also study a solvable coupled system, i.e. a harmonic oscillator interacting with bosonic reservoirs, in Sec. III. We show the relaxation from factorized initial states to a Gaussian GGE like (1), and generalize the result to correlated initial states. Concluding remarks are given in Sec. IV, and some involved calculations demonstrating the decay of the cumulants leading to the Gaussification and the proofs of some mixing properties of the Gaussian GGE, which are key to the relaxation, are shown in Appendices A–C.

The equilibration of an exactly solvable quadratic model similar to the one we study in Sec. III was discussed in Ref. ref:EquilibrationExact-PagelAlvermannFehskePRE on the basis of its exact solution. But only the reduced dynamics of the central harmonic oscillator, with the reservoirs’ degrees of freedom traced out, was analyzed there with a factorized initial state with the reservoirs being in a Gaussian state. In Sec. III, in contrast, we will analyze the evolution of the state of the whole system, including the reservoirs, from a correlated initial state, which is not Gaussian in general. The total system relaxes to a GGE like (1).

II Free Bosonic Field

Let us start with a free bosonic field in $D$-dimensional space. The annihilation and creation operators $\hat{b}_{\bm{k}}$ and $\hat{b}_{\bm{k}}^{\dagger}$ for bosons with momentum $\bm{k}$ obey the canonical commutation relations $[\hat{b}_{\bm{k}},\hat{b}_{\bm{k}^{\prime}}]=0$, $[\hat{b}_{\bm{k}},\hat{b}_{\bm{k}^{\prime}}^{\dagger}]=\delta^{D}(\bm{k}-\bm{k}^{\prime})$, and the Hamiltonian of the system is given by

[TABLE]

with a dispersion relation $\omega_{k}\geq 0$. We set $\hbar=1$.

II.1 Gaussian Initial State

We first focus on Gaussian initial states ref:ContVarQI ; ref:GaussianQI . The density operator for a Gaussian state is formally given by ref:GaussFidelity

[TABLE]

where $\langle{}\cdots{}\rangle$ denotes the expectation value in the state $\hat{\rho}_{0}$, and $\Theta_{\bm{k}\bm{k}^{\prime}}$ is a $2\times 2$ matrix satisfying $(\Theta_{\bm{k}\bm{k}^{\prime}})^{\dagger}=\Theta_{\bm{k}^{\prime}\bm{k}}$. Note however that, for an infinitely extended system like the present bosonic field in the thermodynamic limit, such a density operator is not normalizable and is mathematically ill-defined. Instead, in the $C^{*}$-algebraic approach to quantum field theory and quantum statistical mechanics Bratteli1 ; ref:NESS-AschbacherPillet-JSP ; ref:NESS-AschbacherJaksicPautratPillet ; ref:NESS-Tasaki , states are rigorously characterized by characteristic functionals, i.e. the generating functionals of correlation functions. The characteristic functional of a Gaussian state is Gaussian, and for the present bosonic field it reads as

[TABLE]

where $V_{\bm{k}\bm{k}^{\prime}}$ is the covariance matrix of the Gaussian state, defined by

[TABLE]

and satisfying $(V_{\bm{k}\bm{k}^{\prime}})^{\dagger}=V_{\bm{k}^{\prime}\bm{k}}$. The covariance matrix $V_{\bm{k}\bm{k}^{\prime}}$ is formally related to the matrix $\Theta_{\bm{k}\bm{k}^{\prime}}$ in the Gaussian density operator (3) through

[TABLE]

with $Z=\begin{pmatrix}1&0\\ 0&-1\end{pmatrix}$ ref:GaussFidelity , where $\Theta$ and $V$ are understood as infinite-dimensional matrices with their rows and columns labeled by the continuous indices $\bm{k}$ and $\bm{k}^{\prime}$ as well note:Matrix .

We impose physical conditions on the initial state $\hat{\rho}_{0}$. We assume that the correlation functions involved in the covariance matrix $V_{\bm{k}\bm{k}^{\prime}}$ are endowed with the following structures:

[TABLE]

where $f_{\bm{k}}$, $F_{\bm{k}\bm{k}^{\prime}}$, $g_{\bm{k}}$, and $G_{\bm{k}\bm{k}^{\prime}}$ are “regular” functions of $\bm{k}$ and $\bm{k}^{\prime}$ (so that the Riemann-Lebesgue lemma works later note:RiemannLebesgueLemma ). We also assume that the first moment $\langle\hat{b}_{\bm{k}}\rangle$ is regular in $\bm{k}$. For instance, in the case of the canonical ensemble $\hat{\rho}_{0}\propto\text{e}^{-\hat{H}/k_{B}T}$, we have $f_{\bm{k}}=1/(\text{e}^{\omega_{k}/k_{B}T}-1)$, which is the Bose distribution function, with the other components $F_{\bm{k}\bm{k}^{\prime}}$, $g_{\bm{k}}$, $G_{\bm{k}\bm{k}^{\prime}}$, and $\langle\hat{b}_{\bm{k}}\rangle$ vanishing. We are however interested in more general Gaussian states than the canonical ensemble.

The first contribution $f_{\bm{k}}\delta^{D}(\bm{k}-\bm{k}^{\prime})$ to the correlation (7) represents the translationally invariant component in space, in the corresponding correlation function $\langle[\hat{\psi}^{\dagger}(\bm{r})-\langle\hat{\psi}^{\dagger}(\bm{r})\rangle][\hat{\psi}(\bm{r}^{\prime})-\langle\hat{\psi}(\bm{r}^{\prime})\rangle]\rangle$ in the configuration space, where $\hat{\psi}(\bm{r})$ is the field operator defined by $\hat{\psi}(\bm{r})=\int\text{d}^{D}\bm{k}\,\hat{b}_{\bm{k}}\text{e}^{-\text{i}\bm{k}\cdot\bm{r}}/\sqrt{(2\pi)^{D}}$. In particular, it yields a particle distribution $\langle\hat{\psi}^{\dagger}(\bm{r})\hat{\psi}(\bm{r})\rangle$ uniform over the space. Note that in general a sharp and narrow spectrum in the momentum space like the delta function in (7) corresponds to a widely spread distribution in the configuration space. The second contribution $F_{\bm{k}\bm{k}^{\prime}}$, on the other hand, adds a nonuniformity to the particle distribution, and rules the single-particle coherence. The regularity of $F_{\bm{k}\bm{k}^{\prime}}$ implies that the coherence length of the single-particle correlation is finite, i.e. $\langle[\hat{\psi}^{\dagger}(\bm{r})-\langle\hat{\psi}^{\dagger}(\bm{r})\rangle][\hat{\psi}(\bm{r}^{\prime})-\langle\hat{\psi}(\bm{r}^{\prime})\rangle]\rangle$ decays as $|\bm{r}-\bm{r}^{\prime}|$ increases. In other words, particles are distributed all over the space (the total number of particles is infinite), while the single-particle coherence length is finite.

Similarly, the first contribution $g_{\bm{k}}\delta^{D}(\bm{k}+\bm{k}^{\prime})$ to the other correlation (8) represents the translationally invariant component in the pair correlation $\langle[\hat{\psi}(\bm{r})-\langle\hat{\psi}(\bm{r})\rangle][\hat{\psi}(\bm{r}^{\prime})-\langle\hat{\psi}(\bm{r}^{\prime})\rangle]\rangle$ in the configuration space. This pair correlation decays as $|\bm{r}-\bm{r}^{\prime}|$ increases, due to the regularity of $G_{\bm{k}\bm{k}^{\prime}}$. That is, the pairing is allowed everywhere in space (the total number of pairs can be infinite), while the size of each pair is finite.

Under these conditions, the covariance matrix (5) of the initial Gaussian state is reduced to

[TABLE]

Finally, the regularity of $\langle\hat{b}_{\bm{k}}\rangle$ means that a nonvanishing first moment $\langle\hat{\psi}(\bm{r})\rangle$ is allowed only locally in space. These are our physical conditions on the initial state $\hat{\rho}_{0}$. The particles are distributed all over the space, while the correlations are assumed to be of finite range, excluding peculiar long-range correlations.

We note some conditions that have to be satisfied by $f_{\bm{k}}$, $F_{\bm{k}\bm{k}^{\prime}}$, $g_{\bm{k}}$, and $G_{\bm{k}\bm{k}^{\prime}}$. First, it is clear from the definitions of these functions in (7) and (8) that they possess the symmetries $f_{\bm{k}}=f_{\bm{k}}^{*}$, $F_{\bm{k}\bm{k}^{\prime}}=F_{\bm{k}^{\prime}\bm{k}}^{*}$, $G_{\bm{k}\bm{k}^{\prime}}=G_{\bm{k}^{\prime}\bm{k}}$. Second, the covariance matrix $V_{\bm{k}\bm{k}^{\prime}}$ should satisfy $V+Z/2\geq 0$ as an infinite-dimensional matrix with its rows and columns labeled by the continuous indices $\bm{k}$ and $\bm{k}^{\prime}$ as well, and with $Z$ here representing $Z\delta^{D}(\bm{k}-\bm{k}^{\prime})$. It is an uncertainty relation expressed in terms of the covariance matrix ref:ContVarQI ; ref:GaussianQI .

II.2 Relaxation to GGE

We are now ready to study the dynamics of the system evolving according to the Hamiltonian (2) from the Gaussian initial state (4) with the covariance matrix (9). In the Heisenberg picture, the characteristic functional of the state at time $t$ is calculated as

[TABLE]

where $\langle{}\cdots{}\rangle_{t}$ denotes the expectation value in the state $\hat{\rho}(t)=\text{e}^{-\text{i}\hat{H}t}\hat{\rho}_{0}\text{e}^{\text{i}\hat{H}t}$ at time $t$. For the Gaussian initial state (4) with the covariance matrix (9), it reads

[TABLE]

In the long-time limit $t\to\infty$, the system relaxes to

[TABLE]

by the Riemann-Lebesgue lemma note:RiemannLebesgueLemma . In terms of the density operator, it formally means

[TABLE]

recalling the conversion formula (6) applied to the diagonal covariance matrix $V_{\bm{k}\bm{k}^{\prime}}^{(\text{GGE})}=(\frac{1}{2}+f_{\bm{k}})\delta^{D}(\bm{k}-\bm{k}^{\prime})\openone_{2}$ in (12). If $f_{\bm{k}}$ is the Bose distribution function $f_{\bm{k}}=1/(\text{e}^{\omega_{k}/k_{B}T}-1)$, this equilibrium state coincides with the canonical ensemble $\hat{\rho}_{\text{can}}\propto\text{e}^{-\hat{H}/k_{B}T}$ at temperature $T$. If $f_{\bm{k}}$ is different from the Bose distribution function, the equilibrium state $\hat{\rho}_{\text{GGE}}$ is a GGE, with a set of integrals of motion $\hat{I}_{\bm{k}}=\hat{b}_{\bm{k}}^{\dagger}\hat{b}_{\bm{k}}$. That is why we have named the equilibrium density operator $\hat{\rho}_{\text{GGE}}$ in (13).

The mechanism for the relaxation to the GGE in this simple example is clear. It is due to the Riemann-Lebesgue lemma with the diagonal singularity $\delta^{D}(\bm{k}-\bm{k}^{\prime})$ VANHOVE1955901 ; ANTONIOU1997737 in the correlation (7) in the covariance matrix $V_{\bm{k}\bm{k}^{\prime}}$ of the initial Gaussian state. The translationally invariant component with the delta function $\delta^{D}(\bm{k}-\bm{k}^{\prime})$ in the normal correlation (7) survives in the long-time limit, while the other components decay away. We notice that the equilibration time depends not only on $F_{\bm{k}\bm{k}^{\prime}}$, $g_{\bm{k}}$, $G_{\bm{k}\bm{k}^{\prime}}$, and $\langle\hat{b}_{\bm{k}}\rangle$, characterizing the initial state, but also on $J_{\bm{k}}$, related to observables. In taking the limit in (12), we have assumed the regularity of $J_{\bm{k}}$ so that the Riemann-Lebesgue lemma works. It physically means that our observables are assumed to be spatially localized (of finite size) note:Locality . If, for instance, one considers a nonlocal observable spreading over a very large region in space, it corresponds to taking a very narrow function $J_{\bm{k}}$ in the momentum space, and the relaxation due to the Riemann-Lebesgue lemma becomes very slow. In this way, the time scale for the relaxation to the GGE is ruled by the locality of the observables of interest ($J_{\bm{k}}$), as well as the locality of the correlations in the initial state ($F_{\bm{k}\bm{k}^{\prime}}$, $g_{\bm{k}}$, $G_{\bm{k}\bm{k}^{\prime}}$, and $\langle\hat{b}_{\bm{k}}\rangle$).

In Appendix A, we provide explicit formulas for the time evolutions of the correlations for Gaussian spectra with a quadratic dispersion relation $\omega_{k}$. The correlations actually decay except for the translationally invariant component of the normal correlation.

II.3 Non-Gaussian Initial State

We have so far focused on Gaussian initial states. Let us generalize the analysis to non-Gaussian initial states. The non-Gaussianity is characterized by the higher-order cumulants in the characteristic functional $\chi_{0}[J,J^{*}]$ of the initial state. For instance, suppose that there exists a third-order cumulant like

[TABLE]

in the cumulant expansion of $\ln\chi_{0}[J,J^{*}]$. As we did for the second-order correlations, we allow this third-order correlation $K_{\bm{k}_{1}\bm{k}_{2}\bm{k}_{3}}$ to possess a translationally invariant component proportional to $\delta^{D}(\bm{k}_{1}+\bm{k}_{2}-\bm{k}_{3})$; otherwise, it is assumed to be free from singularity. Namely, it assumes the form

[TABLE]

with $\bar{K}_{\bm{k}_{1}\bm{k}_{2}}$ and $\tilde{K}_{\bm{k}_{1}\bm{k}_{2}\bm{k}_{3}}$ being regular functions of the momenta. It physically means that this third-order correlation is allowed to exist everywhere in space, while its correlation lengths are finite. This cumulant evolves in time according to the Hamiltonian (2) as

[TABLE]

This decays in the long-time limit $t\to\infty$ for a generic dispersion relation $\omega_{k}$, according to the Riemann-Lebesgue lemma. In this way, any higher-order cumulants decay in the long-time limit $t\to\infty$ under the assumption of finite correlation lengths mentioned above (they are regular apart from the translationally invariant components), and the system relaxes to the Gaussian GGE as (12), even from a non-Gaussian initial state.

See Appendix A, where explicit formulas for the time evolutions of all the cumulants for Gaussian spectra with a quadratic dispersion relation $\omega_{k}$ are provided. All the cumulants except for the translationally invariant component of the second-order normal correlation decay in the long-time limit $t\to\infty$. The results for the Gaussian spectra are of quite general validity, since by the saddle-point approximation the integrals for the cumulants are approximated by Gaussian integrals for large $t$. For any quadratic dispersion relation $\omega_{k}$, the cumulants exhibit power-law decays as shown in (65) and (66).

If the third- and higher-order correlations exist only locally, with no translationally invariant components, we can show the relaxation to the Gaussian GGE in a simpler way as follows. First, we observe that the Gaussian GGE in (13) is mixing with respect to the Hamiltonian (2) Bratteli1 , namely,

[TABLE]

holds for any local observables $\hat{A}$, $\hat{B}$, and $\hat{C}$, where $\hat{B}(t)=\text{e}^{\text{i}\hat{H}t}\hat{B}\text{e}^{-\text{i}\hat{H}t}$ is the Heisenberg operator corresponding to $\hat{B}$, and $\langle{}\cdots{}\rangle_{\text{GGE}}$ denotes the expectation value in the Gaussian GGE in (13). See Appendix B for a proof. This tells us that any locally perturbed GGE, which is non-Gaussian in general, returns to the Gaussian GGE ref:ReturnEq :

[TABLE]

where $\hat{L}_{j}$ represent the local perturbations, satisfying $\sum_{j}\hat{L}_{j}^{\dagger}\hat{L}_{j}=\hat{\openone}$ for the normalization of $\hat{\rho}_{0}$. Indeed, thanks to the mixing property (17), we get

[TABLE]

for any local observable $\hat{A}$, where $\langle{}\cdots{}\rangle$ denotes the expectation value in the initial state $\hat{\rho}_{0}$ in (18). This proves the relaxation to the Gaussian GGE in (18).

If the third-order cumulant $K_{\bm{k}_{1}\bm{k}_{2}\bm{k}_{3}}$ in (14) in the initial state contains a contribution proportional to $\delta(\omega_{k_{1}}+\omega_{k_{2}}-\omega_{k_{3}})$, it is stationary and survives in the long-time limit $t\to\infty$. This adds integrals of motion like $\delta(\omega_{k_{1}}+\omega_{k_{2}}-\omega_{k_{3}})\hat{b}_{\bm{k}_{1}}\hat{b}_{\bm{k}_{2}}\hat{b}_{\bm{k}_{3}}^{\dagger}$ to the exponent of the GGE in (13). If, however, the initial state is prepared irrespective of the Hamiltonian $\hat{H}$ of the system, e.g. in a quench scenario, it is typically impossible that the initial state is equipped with such a delicate structure in the correlations, matching the characteristics of the Hamiltonian $\hat{H}$. The same applies to the higher-order cumulants. The relaxation to the Gaussian GGE in (13) is therefore quite general for generic initial states.

III Harmonic Oscillator Coupled with Bosonic Reservoirs

In the previous section, we have observed the relaxation of the free bosonic field to a Gaussian GGE. Let us now look at a coupled system. As an interesting and nontrivial example, we consider a harmonic oscillator coupled with multiple bosonic reservoirs. The Hamiltonian is given by ref:OscillatorFriedrichsBarnett ; ref:ExactDiagonalizationCaldeira ; ref:FactorizeModel ; ref:FTExactMonnai ; Monnai1

[TABLE]

with

[TABLE]

The operators $\hat{a}$ and $\hat{b}_{\bm{k}\ell}$ are the canonical operators for the harmonic oscillator and the bosonic reservoirs, respectively, satisfying the canonical commutation relations

[TABLE]

$\Omega>0$ is the frequency of the harmonic oscillator, $\omega_{k\ell}\geq 0$ is the dispersion relation for the $\ell$th reservoir, $u_{\bm{k}\ell}$ is the form factor of the interaction between the harmonic oscillator and the $\ell$th reservoir note:FormFactor , and $\lambda$ is a coupling constant note:BoundState .

It is an exactly solvable quadratic model, and we can diagonalize the Hamiltonian (20) as ref:OscillatorFriedrichsBarnett ; ref:ExactDiagonalizationCaldeira ; ref:FTExactMonnai ; ref:FanoDiagonalization

[TABLE]

where the normal modes are given by

[TABLE]

with

[TABLE]

Here, the self-energy function on the complex $z$ plane is given by

[TABLE]

and the coefficients satisfy the orthogonality

[TABLE]

and the completeness

[TABLE]

The completeness (33) ensures that the operators $\hat{A}_{\bm{k}\ell}$ satisfy the canonical commutation relations

[TABLE]

and the orthogonality (31)–(32) allows us to invert the relation (27) as

[TABLE]

III.1 Factorized Initial State

We first consider the evolution of the system from a factorized initial state

[TABLE]

with no correlations among the harmonic oscillator $\hat{\rho}_{S}$ and the reservoirs $\hat{\rho}_{\ell}$. The characteristic function of the initial state of the harmonic oscillator $\hat{\rho}_{S}$ and the characteristic functional of the initial state of the reservoirs $\hat{\rho}_{B}=\mathop{\bigotimes}_{\ell}\hat{\rho}_{\ell}$ are defined respectively by

[TABLE]

where $\langle{}\cdots{}\rangle$ denotes the expectation value in the initial state $\hat{\rho}_{0}$ in (36). The initial state of the harmonic oscillator $\hat{\rho}_{S}$ is arbitrary, while the initial state of each bosonic reservoir is a (non-Gaussian) state like the one considered in Sec. II, under the assumption of finite correlation lengths. More specifically, in the cumulant expansion of the characteristic functional of the reservoirs $\chi_{B}[\eta,\eta^{*}]$,

[TABLE]

the covariance matrix $V_{\bm{k}\bm{k}^{\prime}}^{(\ell)}$ and the third-order cumulant $K_{\bm{k}_{1}\bm{k}_{2}\bm{k}_{3}}^{(\ell)}$ of the $\ell$th reservoir are endowed with the same structures as those in (9) and (15), respectively, with $f_{\bm{k}}$, $F_{\bm{k}\bm{k}^{\prime}}$, $g_{\bm{k}}$, $G_{\bm{k}\bm{k}^{\prime}}$, $\bar{K}_{\bm{k}_{1}\bm{k}_{2}}$, and $\tilde{K}_{\bm{k}_{1}\bm{k}_{2}\bm{k}_{3}}$ replaced by $f_{\bm{k}}^{(\ell)}$, $F_{\bm{k}\bm{k}^{\prime}}^{(\ell)}$, $g_{\bm{k}}^{(\ell)}$, $G_{\bm{k}\bm{k}^{\prime}}^{(\ell)}$, $\bar{K}_{\bm{k}_{1}\bm{k}_{2}}^{(\ell)}$, and $\tilde{K}_{\bm{k}_{1}\bm{k}_{2}\bm{k}_{3}}^{(\ell)}$, which are all assumed to be regular functions of the momenta. The other (higher-order) cumulants of the reservoirs can also possess translationally invariant components, but otherwise they are assumed to be regular in momenta.

Starting from such a factorized initial state $\hat{\rho}_{0}$ in (36), the characteristic functional of the state of the total system evolves according to the Hamiltonian (26) as

[TABLE]

where $\chi_{S}(\xi,\xi^{*})$ and $\chi_{B}[\eta,\eta^{*}]$ are the characteristic function and the characteristic functional of the initial states of the harmonic oscillator and of the reservoirs given respectively in (37) and (38), and

[TABLE]

In the long-time limit $t\to\infty$, we get

[TABLE]

due to the Riemann-Lebesgue lemma, while

[TABLE]

by recalling

[TABLE]

Therefore, the characteristic functional of the total system (40) behaves asymptotically as

[TABLE]

where we have used the normalization condition $\chi_{S}(0,0)=1$ of the initial state $\hat{\rho}_{S}$ of the harmonic oscillator. Notice that the asymptotic characteristic functional $\chi_{B}[J\text{e}^{\text{i}\omega t},J^{*}\text{e}^{-\text{i}\omega t}]$ in (46) is essentially the same as (10), and the results in Sec. II immediately apply. Each reservoir relaxes as (12), and the characteristic functional of the total system in (46) further relaxes to

[TABLE]

In terms of the density operator, the stationary state is formally given by

[TABLE]

It is a Gaussian state: the total system relaxes to the Gaussian state $\hat{\rho}_{\text{NESS}}$, even from a non-Gaussian initial state $\hat{\rho}_{0}$, under the condition of finite correlation lengths in the reservoirs.

Let us comment on some physical aspects of the stationary state $\hat{\rho}_{\text{NESS}}$.

III.1.1 Nonequilibrium Steady State

In the presence of two or more reservoirs, the stationary state (48) is a nonequilibrium steady state (NESS) ref:NESS-AschbacherPillet-JSP ; ref:NESS-AschbacherJaksicPautratPillet ; ref:NESS-Tasaki ; ref:FTExactMonnai ; Monnai1 , in which steady currents flow among the reservoirs through the harmonic oscillator. That is why we have named the stationary state $\hat{\rho}_{\text{NESS}}$ in (48). For instance, the energy current

[TABLE]

flowing into the $\ell$th reservoir per time is estimated in the stationary state $\hat{\rho}_{\text{NESS}}$ to be

[TABLE]

where

[TABLE]

and we have introduced

[TABLE]

which is normalized as $\int_{0}^{\infty}\text{d}\omega\,|\alpha(\omega)|^{2}=1$, and $|\alpha(\omega)|^{2}\to\delta(\omega-\Omega)$ in the weak-coupling limit $\lambda\to 0$. If $f_{\bm{k}}^{(\ell)}$ is isotropic and is given by $f_{\bm{k}}^{(\ell)}=f_{\ell}(\omega_{k\ell})$ [e.g., in the case of the canonical ensemble, $f_{\ell}(\omega)$ is the Bose distribution function], we have $\mathcal{F}_{\ell}(\omega)=f_{\ell}(\omega)\Gamma_{\ell}(\omega)$, and the current (50) is simplified to

[TABLE]

The steady current flows by the difference between $f_{\ell}(\omega)$ and $f_{\ell^{\prime}}(\omega)$. In the weak-coupling limit $\lambda\to 0$ (more precisely, in the van Hove limit $\lambda\to 0$ keeping the scaled time $\tau=\lambda^{2}t$ finite), it is further simplified to

[TABLE]

which is a standard Landauer formula but with the reservoirs in GGEs.

III.1.2 Equilibration of the Subsystem

The long-time limit (47) shows that the system forgets the initial state $\hat{\rho}_{S}$ of the harmonic oscillator and relaxes to the stationary state $\hat{\rho}_{\text{NESS}}$ independent of $\hat{\rho}_{S}$. Let us look in which state the harmonic oscillator equilibrates. Recalling the inversion formula (35), the characteristic function of the harmonic oscillator $\chi_{t}^{(S)}(\xi,\xi^{*})$ can be extracted from the characteristic functional of the total system $\chi_{t}[J,J^{*}]$ in (40). It relaxes to

[TABLE]

The corresponding equilibrium density operator is given by

[TABLE]

with

[TABLE]

If the harmonic oscillator is immersed in a single reservoir and $f_{\bm{k}}=f(\omega_{k})=1/(\text{e}^{\omega_{k}/k_{B}T}-1)$ is the Bose distribution function, the equilibrium density operator $\hat{\rho}_{\text{NESS}}^{(S)}$ in (58) is reduced to the thermal state $\hat{\rho}_{\text{NESS}}^{(S)}\propto\text{e}^{-\hat{H}_{S}/k_{B}T}$ at the same temperature $T$ as that of the reservoir in the weak-coupling limit $\lambda\to 0$. In general, the reservoirs do not relax to the canonical state, with $f_{\bm{k}}^{(\ell)}$ being different from the Bose distribution functions, but in any case the equilibrium state of the harmonic oscillator $\hat{\rho}_{\text{NESS}}^{(S)}$ in (58) looks like a canonical state with an effective temperature $\Omega/k_{B}\theta$, which depends on the GGE characterized by $f_{\bm{k}}^{(\ell)}$.

III.2 Correlated Initial State

In the previous subsection, we considered factorized initial states (36), with no correlations among the harmonic oscillator and the reservoirs. Even if there are some correlations in the initial state, the system relaxes to the same stationary state $\hat{\rho}_{\text{NESS}}$ as that given in (48), as long as the initial correlations are just local. Namely, we consider, instead of the factorized initial state (36), a correlated initial state

[TABLE]

where the factorized state (36) is perturbed by local operators $\hat{L}_{j}$, which satisfy $\sum_{j}\hat{L}_{j}^{\dagger}\hat{L}_{j}=\openone$ and induce correlations among the harmonic oscillator and the reservoirs. Still, the system relaxes to

[TABLE]

where $\hat{\rho}_{\text{NESS}}$ is the same stationary state as the one presented in (48). A proof is provided in Appendix C.

IV Summary

We have shown a scenario of the relaxation to GGE for integrable models which can be mapped to free bosonic fields in the thermodynamic limit. The unitary transformation to the free bosonic fields would be regarded as a quantum counterpart of the canonical transformation to an assembly of harmonic oscillators for classical systems in the context of the Liouville integrability. Then, the field operators in the characteristic functional just acquire oscillating factors in the evolution as in (10), and the diagonal singularity yields the stationary state. As a result, only the gauge invariant terms such as the occupation number $\hat{I}_{\bm{k}}=\hat{b}_{\bm{k}}^{\dagger}\hat{b}_{\bm{k}}$ survive in the GGE in (12) in the long-time limit $t\rightarrow\infty$.

Moreover, the GGE (12) is a simple Gaussian state: only the quadratic gauge invariant terms $\hat{I}_{\bm{k}}=\hat{b}_{\bm{k}}^{\dagger}\hat{b}_{\bm{k}}$ contribute to the GGE, as long as the initial (not necessarily Gaussian) state fulfills a few physical conditions, where the presence of anomalous correlations is allowed. In contrast to the previous works ref:CramersDawsonEisertOsbornePRL2008 ; ref:CramerEisert-CentralLimitTheorem ; ref:GogolinEisert-GaussificationFermi ; Farrelly1 ; Murthy1 , which proved the Gaussification for large but finite systems, we have directly analyzed the systems in the thermodynamic limit, with continuous spectra. This greatly simplifies the analysis, and in this picture, the mechanism for the equilibration is due to the Riemann-Lebesgue lemma and the diagonal singularity. We have solved the evolutions of the states of the whole systems exactly, and have shown the Gaussification for rather general (non-Gaussian) initial states (particles can be distributed all over the space and correlations can be present everywhere in space, but the correlations should be of finite range), even in the presence of initial correlations among systems.

We stress that the observables are also relevant to the relaxation. In particular, the locality of the observables is important, for the Riemann-Lebesgue lemma to work. It is known that the relaxation times for typical systems and/or typical settings are typically short relaxation1 ; Goldstein3 ; Monnai2 ; relaxation2 ; Goldstein1b ; relaxation3 irrespective of the system characteristics or the observables. For specific (atypical) systems, however, it is not the case, and the relaxation time depends on the choice of the observable as well as the initial state.

In this work, we have studied the relaxation to a GGE for time-independent Hamiltonians. It would be interesting, as a potential future subject, to extend the discussion to the stationary states of periodically driven systems, on the basis of the Floquet theory Lazarides1 ; Lazarides2 ; DAlessio1 ; Haldar1 .

Acknowledgements.

This work was supported by the Grants-in-Aid for Scientific Research (C) (No. 18K03467 and No. 18K03470) and for Fostering Joint International Research (B) (No. 18KK0073) both from the Japan Society for the Promotion of Science (JSPS), and by the Waseda University Grant for Special Research Projects (No. 2018K-262).

Appendix A Decay of Correlations

The Gaussification is a consequence of the decay of the cumulants except for the occupation numbers. We have attributed the mechanism behind the decay to the Riemann-Lebesgue lemma, under the assumptions (i)–(v) listed in Sec. I. Let us here compute the evolutions of the cumulants explicitly for analytically tractable specific forms of spectra to get an idea on how the Riemann-Lebesgue lemma works and how the correlations decay in the long-time limit.

Let us first look at the second-order cumulants in (11). There are four different types. The translationally invariant component of the normal correlation represented by $f_{\bm{k}}$, which is related to the occupation numbers, does not decay and survives in the GGE: $\bar{\mathcal{K}}_{1,1}(t)=\int\text{d}^{D}\bm{k}\,f_{\bm{k}}J_{\bm{k}}J_{\bm{k}}^{*}$ in (11) is independent of time $t$. The remaining (local) component of the normal correlation represented by $F_{\bm{k}\bm{k}^{\prime}}$ evolves in time as $\tilde{\mathcal{K}}_{1,1}(t)=\int\text{d}^{D}\bm{k}\int\text{d}^{D}\bm{k}^{\prime}\,F_{\bm{k}\bm{k}^{\prime}}J_{\bm{k}}J_{\bm{k}^{\prime}}^{*}\text{e}^{\text{i}(\omega_{k}-\omega_{k^{\prime}})t}$. In general, the cumulant involving $m_{+}$ creation operators $\hat{b}_{\bm{k}}^{\dagger}$ and $m_{-}$ annihilation operators $\hat{b}_{\bm{k}}$ evolves in time as

[TABLE]

As we discussed below (13), not only the spectrum $K_{\bm{k}_{1}\ldots\bm{k}_{m_{+}}\bm{k}_{1}^{\prime}\ldots\bm{k}_{m_{-}}^{\prime}}$ of the correlation function but also observables are relevant to the evolution of the cumulant through $J_{\bm{k}}$. To generate the expectation value of an observable, we take the derivatives of the characteristic functional $\chi_{t}[J,J^{*}]$ with respect to $J_{\bm{k}}$ and $J_{\bm{k}}^{*}$, through which $J_{\bm{k}}$ and $J_{\bm{k}}^{*}$ are replaced by the Fourier spectrum of the observable note:Locality . In this sense, the $\bm{k}$ dependence of $J_{\bm{k}}$ represents the spectra of observables of interest. Let us consider a specific form of the overall spectrum for the evolution of the cumulant, i.e. a Gaussian spectrum with a quadratic dispersion relation $\omega_{k}=\frac{1}{2}a^{2}k^{2}$, to facilitate explicit calculation and to get an idea on how the cumulants decay. In the case of a symmetric Gaussian spectrum, the local component evolves as

[TABLE]

The translationally invariant component with such a symmetric Gaussian spectrum, on the other hand, evolves as

[TABLE]

The local component of any $m$th-order cumulant ($m=m_{+}+m_{-}$) decays asymptotically as

[TABLE]

while the translationally invariant component as

[TABLE]

These show that (i) the cumulants decay except for the translationally invariant component $f_{\bm{k}}$ of the normal correlation $(m_{+},m_{-})=(1,1)$. This leads to the Gaussification. In addition, we see that (ii) the translationally invariant components decay more slowly than the local components, and that (iii) the decays of the translationally invariant components of the gauge-invariant cumulants with $m_{+}=m_{-}$ are even slower.

Even if the real spectra for the evolutions of the cumulants are not of such simple Gaussian form, the saddle-point approximation in estimating the integrals for large $t$ yields Gaussian integrals. The integrands after the saddle-point approximation might not be so simple as the symmetric Gaussian considered above, but in any case, for generic (not necessarily symmetric) Gaussian spectra with the quadratic dispersion relation $\omega_{k}$, the cumulants decay asymptotically as (65) and (66). The results (65) and (66) are therefore of quite general validity.

Appendix B Mixing of GGE

Here, we prove the mixing property (17) of the GGE in (13) with the Hamiltonian $\hat{H}$ in (2). Recall the characteristic functional of the GGE in (12),

[TABLE]

where

[TABLE]

and consider

[TABLE]

which is the generating functional for the correlation functions of the type (17). In the GGE characterized by (67), it reads

[TABLE]

by the Riemann-Lebesgue lemma. This shows the factorization (17) for any local observables $\hat{A}$, $\hat{B}$, and $\hat{C}$.

Appendix C Relaxation from Correlated Initial States

Here, we prove the relaxation (61) from the correlated initial state (60). We actually look at the evolution of the characteristic functional of the state of the total system driven by the Hamiltonian $\hat{H}$ in (26),

[TABLE]

where

[TABLE]

and $\langle{}\cdots{}\rangle$ denotes the expectation value in the factorized state $\hat{\rho}_{S}\otimes\hat{\rho}_{B}$ in (36). We introduce the generating functional for the correlation functions in (72),

[TABLE]

It is reduced to

[TABLE]

where $\chi_{S}(\xi,\xi^{*})$ and $\chi_{B}[\eta,\eta^{*}]$ are the characteristic function of the state $\hat{\rho}_{S}$ of the harmonic oscillator in (37) and the characteristic functional of the state $\hat{\rho}_{B}$ of the reservoirs in (38), respectively, with

[TABLE]

For $t\to\infty$, due to the Riemann-Lebesgue lemma, it behaves asymptotically as

[TABLE]

with

[TABLE]

Recall (43) and (44) for the factorized case. Inserting the asymptotic behavior of $\eta_{\bm{k}\ell}(t)$ in the cumulant expansion of $\chi_{B}[\eta,\eta^{*}]$ in (39), we have

[TABLE]

where $\chi_{\text{NESS}}[J,J^{*}]$ is given in (47). Thus, (78) further relaxes to

[TABLE]

This shows that the characteristic functional (72) for the correlated initial state $\hat{\rho}_{0}$ in (60) relaxes to

[TABLE]

for any local perturbations $\hat{L}_{j}$, and proves (61), under the normalization condition $\sum_{j}\langle\hat{L}_{j}^{\dagger}\hat{L}_{j}\rangle=1$ for the correlated initial state $\hat{\rho}_{0}$.