Metastability associated with many-body explosion of eigenmode expansion coefficients

Let us consider a system with discrete states $n=1,2,\dots,D$. The energy of the state $n$ is denoted by $E_{n}$. The probability distribution $P_{n}(t)$ evolves following the master equation

where the matrix $G$ is the generator of the stochastic dynamics. For $n\neq m$, $G_{n,m}$ represents the transition rate from the state $m$ to $n$, and $G_{n,n}=-\sum_{m(\neq n)}G_{m,n}$. When the system is coupled to a thermal bath at the inverse temperature $\beta$, $G$ satisfies the detailed-balance condition, $G_{n,m}e^{-\beta E_{m}}=G_{m,n}e^{-\beta E_{n}}$.

The detailed-balance condition ensures that $G$ is diagonalizable with real (and non-negative) eigenvalues since it is made symmetric via the similarity transformation $G\to G^{\prime}=e^{\beta H/2}Ge^{-\beta H/2}$ with $H=\mathrm{diag}(E_{1},E_{2},\dots,E_{D})$. We here assume that the stationary state is unique and eigenvalues are ordered as $0=\Lambda_{0}<\Lambda_{1}\leq\Lambda_{2}\dots\leq\Lambda_{D-1}$. The right and left eigenvectors for an eigenvalue $\Lambda_{\alpha}$ are respectively denoted by $\vec{\Phi}^{(\alpha)}$ and $\vec{\Pi}^{(\alpha)}$: $G\vec{\Phi}^{(\alpha)}=\Lambda_{\alpha}\vec{\Phi}^{(\alpha)}$ and $G^{\mathrm{T}}\vec{\Pi}^{(\alpha)}=\Lambda_{\alpha}\vec{\Pi}^{(\alpha)}$.

Let us normalize right eigenvectors by using the $L^{1}$ norm $\|\cdot\|_{1}$ as follows:

where $\Phi^{(\alpha)}_{n}$ are $n$th component of $\vec{\Phi}^{(\alpha)}$. We will later see that this normalization is convenient for our purpose. Since the dual norm of the $L^{1}$ norm is the $L^{\infty}$ norm $\|\cdot\|_{\infty}$ in the sense that $|\vec{A}\cdot\vec{B}|\leq\|\vec{A}\|_{\infty}\|\vec{B}\|_{1}$ for any set of vectors $\vec{A}$ and $\vec{B}$, we normalize left eigenvectors by using $L^{\infty}$ norm as follows:

where $\Pi^{(\alpha)}_{n}$ are $n$th component of $\vec{\Pi}^{(\alpha)}$.

The mode $\alpha=0$ corresponds to the stationary state, and $\Phi^{(0)}_{n}=e^{-\beta E_{n}}/Z$ with $Z=\sum_{n}e^{-\beta E_{n}}$ and $\Pi^{(0)}_{n}=1$ for all $n$. Without loss of generality, we choose the signs of $\vec{\Phi}^{(\alpha)}$ and $\vec{\Pi}^{(\alpha)}$ so that $\max_{n}\Pi^{(\alpha)}_{n}=1$ and $\vec{\Pi}^{(\alpha)}\cdot\vec{\Phi}^{(\alpha)}>0$.

The probability distribution $\vec{P}(t)=(P_{1}(t),\dots,P_{D}(t))^{\mathrm{T}}$ can be expanded in terms of right eigenvectors of $G$ as $\vec{P}(t)=\sum_{\alpha}C_{\alpha}e^{-\Lambda_{\alpha}t}\vec{\Phi}^{(\alpha)}$, where

is an eigenmode expansion coefficient in the initial state.

Obviously, $\Lambda_{\alpha}$ gives the decay rate of $\alpha$th eigenmode, and hence the slowest decay rate is given by the spectral gap $\Lambda_{1}$ of $G$. This observation indicates that slow relaxation should accompany small $\Lambda_{1}$. Following this idea, in previous works Gaveau and Schulman (1987, 1998); Biroli and Kurchan (2001); Macieszczak et al. (2016), metastability is characterized by the presence of small eigenvalues. This argument implicitly assumes that $C_{\alpha}$ is not too large. As we point out below, however, this assumption is typically not satisfied in many-body systems.

For a given initial state $\vec{P}(0)=\sum_{\alpha}C_{\alpha}\vec{\Phi}^{(\alpha)}$, let us define the relaxation time $\tau_{\alpha}$ for each eigenmode. The expectation value of a physical quantity $\vec{O}=(O_{1},O_{2},\dots,O_{D})^{\mathrm{T}}$ with $\|\vec{O}\|_{\infty}=1$ is given by

The contribution from $\alpha$th eigenmode is thus bounded as

where the normalization (2) is used. Thus the contribution from $\alpha$th eigenmode is negligible when $|C_{\alpha}|e^{-\Lambda_{\alpha}t}\ll 1$, and hence it is appropriate to define the eigenmode relaxation time $\tau_{\alpha}$ as the time at which $|C_{\alpha}|e^{-\Lambda_{\alpha}t}=\delta$ for a fixed constant $\delta$. Throughout the paper, we fix $\delta=0.5$.

We then have

which implies that large expansion coefficients can cause the delay of the relaxation ²²2If we choose another normalization condition for right eigenvectors, the criterion of the relaxation $|C_{\alpha}|e^{-\Lambda_{\alpha}t}=\delta$ should be modified. For example, if we employ the normalization $\|\vec{\Phi}^{(\alpha)}\|_{2}:=\left(\sum_{n=1}^{D}|\Phi^{(\alpha)}_{n}|^{2}\right)^{1/2}=1$, we have $|\vec{O}\cdot\vec{\Phi}^{(\alpha)}|\sim D^{1/2}$, and hence the criterion of the relaxation should be $|C_{\alpha}|e^{-\Lambda_{\alpha}t}=\delta D^{-1/2}$ with $\delta$ being a fixed constant. Accordingly, the formula (7) should also be modified. In other words, Eq. (7) is appropriate only for the normalization scheme (2).

Which eigenmodes are responsible for metastable states? To answer this question, let us examine the possible largest value of $\tau_{\alpha}$ over all initial distributions. From Eq. (4), we have an upper bound

This upper bound is always realizable: $|C_{\alpha}|=\Psi_{\alpha}$ when $P_{n}(0)$ is nonzero only for $n$ such that $\Pi^{(\alpha)}_{n}=1$. Thus, we obtain

The eigenmode relaxation times $\{\tilde{\tau}_{\alpha}\}$ tell us about which eigenmodes can compose metastable states. In bra-ket notation, the time evolution operator is written as

For a timescale specified by $t$, we consider a subset $\mathcal{M}$ of eigenmodes such that $\tilde{\tau}_{\alpha}\ll t$ for any $\alpha\notin\mathcal{M}$. Since $\Psi_{\alpha}e^{-\Lambda_{\alpha}t}=\delta e^{\Lambda_{\alpha}(\tilde{\tau}_{\alpha}-t)}\ll 1$ for $\alpha\notin\mathcal{M}$, we have $e^{-Gt}\approx\sum_{\alpha\in\mathcal{M}}\Psi_{\alpha}e^{-\Lambda_{\alpha}t}\ket{\Phi^{(\alpha)}}\bra{\Pi^{(\alpha)}}$, which implies that a metastable state consists of eigenmodes in $\mathcal{M}$.

Let us begin with the simplest model, i.e., $N$ independent particles in a double-well potential under the two-state approximation (see Fig. 1). Each particle $i$ ($i=1,2,\dots,N$) is either in the left or right well, which is expressed by a spin variable $\sigma_{i}=\pm 1$ ($\sigma_{i}=1$ corresponds to the left well). Then $N$-particle state is specified by a set of spin variables $(\sigma_{1},\sigma_{2},\dots,\sigma_{N})$, and there are $D=2^{N}$ different states. The state $\sigma_{i}=+1$ (left well) has the energy $\varepsilon>0$, and the energy barrier undergoing in the transition from $\sigma_{i}=+1$ to $-1$ is given by $E_{B}>0$. By choosing the temperature as the unit of energy (i.e., $\beta=1$), transitions of $i$th particle from $\sigma_{i}=1$ to $-1$ and vice versa happen at rates given by the Arrhenius law $\tau_{0}^{-1}e^{-E_{B}}$ and $\tau_{0}^{-1}e^{-(E_{B}+\varepsilon)}$, respectively. Here, $\tau_{0}$ is a certain microscopic timescale, which is chosen as the unit of time, i.e., $\tau_{0}=1$. This transition rate satisfies the detailed-balance condition for the Hamiltonian $H=\varepsilon N_{+}:=\varepsilon\sum_{i=1}^{N}(\sigma_{i}+1)/2$, where the number of particles in the left well is denoted by $N_{+}$. Similarly, we define $N_{-}=N-N_{+}$ as that in the right well.

In the single-particle problem, the generator is a two-by-two matrix, and it has two eigenvalues $\lambda_{0}=0$ and $\lambda_{1}=e^{-E_{B}}(1+e^{-\varepsilon})$. The right and left eigenvectors are denoted by $\vec{\phi}^{(\alpha)}$ and $\vec{\pi}^{(\alpha)}$ ($\alpha=0$ or 1), respectively. In the $N$-particle case, the eigenvalues $\{\Lambda_{\alpha}\}$ are given by the sum of single-particle eigenvalues as $\Lambda_{\alpha}=\sum_{i=1}^{N}\alpha_{i}\lambda_{1}$, where $\alpha_{i}=0,1$ specifies which eigenmode is occupied by $i$th particle. The corresponding eigenvectors are given by the product of single-particle eigenvectors: $\vec{\Phi}^{(\alpha)}=\vec{\phi}^{(\alpha_{1})}\otimes\vec{\phi}^{(\alpha_{2})}\dots\otimes\vec{\phi}^{(\alpha_{N})}$ and $\vec{\Pi}^{(\alpha)}=\vec{\pi}^{(\alpha_{1})}\otimes\vec{\pi}^{(\alpha_{2})}\dots\otimes\vec{\pi}^{(\alpha_{N})}$.

The maximum expansion coefficient $\Psi_{\alpha}$ in Eq. (8) for each $\alpha$ is explicitly calculated as follows:

When $\sum_{i=1}^{N}\alpha_{i}=O(N)$, $\Psi_{\alpha}=e^{O(N)}$ for $\varepsilon>0$. Expansion coefficients diverge exponentially in $N$. Although the corresponding eigenvalue is extensively large $\Lambda_{\alpha}=O(N)$, it is not correct that this eigenmode decays in a timescale $t\sim\Lambda_{\alpha}^{-1}=O(N^{-1})$. The eigenmode relaxation time $\tilde{\tau}_{\alpha}$ is evaluated as

whose second term is independent of the system size and $\Lambda_{\alpha}$.

In this model, the “all-left” state, i.e., $\sigma_{1}=\sigma_{2}=\dots=\sigma_{N}=1$, is considered to be a many-body metastable state for large $E_{B}$. Expansion coefficients for this state are explicitly calculated, and it is found that $C_{\alpha}=\Psi_{\alpha}$ and $\tau_{\alpha}=\tilde{\tau}_{\alpha}$ for all $\alpha$, i.e., this state gives the maximum expansion coefficients for all eigenmodes.

Even in this simplest case, explosive growth of expansion coefficients occurs. However, $\Psi_{\alpha}$ does not depend on $E_{B}$, and hence large $\Psi_{\alpha}$ is not related to metastability in this non-interacting model.

Let us introduce interactions. We here consider a simple model in which two particles in the same potential well equally interact with each other. The interaction energy is given by $E_{\mathrm{int}}=-(g/N)[N_{+}(N_{+}-1)/2+N_{-}(N_{-}-1)/2]$, where $g>0$ denotes the interaction strength (the scaling $1/N$ ensures the extensivity of the energy). The interaction changes the effective potential felt by a single particle: the energy is written as $E(N_{+})=\varepsilon_{+}N_{+}+\varepsilon_{-}N_{-}$ with $\varepsilon_{+}=\varepsilon-g(N_{+}-1)/N$ and $\varepsilon_{-}=-g(N_{-}-1)/N$. The effective energy barrier from $\sigma_{i}=\pm 1$ to $\mp 1$ is thus given by $E_{B}^{(\pm)}=E_{B}+\varepsilon-\varepsilon_{\pm}$. The model with interactions is obtained by replacing $(\varepsilon,E_{B})$ in the non-interacting model by $(\varepsilon_{+}-\varepsilon_{-},E_{B}^{(+)})$. Obviously, interactions $g>0$ effectively increases the energy barrier and causes slower relaxation, which we call interaction-induced metastability.

The master equation in this model is generated by a $2^{N}$-dimensional generator. Instead of considering the full generator, let us employ the simplified description by focusing on the dynamics of a collective variable $N_{+}\in\{0,1,\dots,N\}$ only. Since the energy only depends on $N_{+}$ in our model, we can exactly obtain the dynamics of $N_{+}$. Let the probability distribution of $N_{+}$ at time $t$ be denoted by $P_{N_{+}}(t)$. Then the master equation for $P_{N_{+}}(t)$ is written as

where $G_{N_{+}\mp 1,N_{+}}=-N_{\pm}e^{-E_{B}^{(\pm)}}$. This generator satisfies the detailed-balance condition of the form

where $F(N_{+})=E(N_{+})-S(N_{+})$ denotes the free energy of a state specified by $N_{+}$. The entropy $S(N_{+})$ is given by $S(N_{+})=\ln(N!/(N_{+}!N_{-}!)$.

The $(N+1)$-dimensional generator in Eq. (13) is a diagonal block of the original $2^{N}$-dimensional full generator. Thus, the full generator has other $2^{N}-(N+1)$ eigenmodes, which are dropped in this description. However, it does not affect the conclusion (see Appendix A).

We now fix $E_{B}=1$ and $\varepsilon=5$. In our model, there is the non-ergodic phase for $g\gtrsim 8.8$, in which we have an extensive free energy barrier $\Delta F\propto N$ and the relaxation time $\tau_{\mathrm{rel}}$ diverges in the thermodynamic limit as $\tau_{\mathrm{rel}}=e^{\Delta F}=e^{O(N)}$. We discuss the non-ergodic phase in Appendix B, and here we focus on the ergodic phase at which the free energy $F(N_{+})$ is a monotonic function of $N_{+}$.

First let us investigate the $g$-dependence of the spectral gap of $G$. Numerical results for different $g$ and $N$ are presented in Fig. 2. It clearly shows that the spectral gap of $G$ is almost independent of the value of $g$ for large $N$, although the finite-size effect is stronger for larger $g$. It means that the spectrum gap $\Lambda_{1}$ does not reflect the increase of the relaxation time due to interactions.

Next, we compute $\{\tilde{\tau}_{\alpha}\}$ at $g=2$ for different $N$ and show the result in Fig. 3. Although the typical value of $\Lambda_{\alpha}$ linearly grows with $N$, that of $\tilde{\tau}_{\alpha}$ looks convergent for large $N$. According to Eq. (9), it means that expansion coefficients typically grows exponentially in $N$, which is already seen in the non-interacting model. In Fig. 3, $\tau_{\alpha}$ for the all-left initial state is also plotted by squares with a dashed line in Fig. 3, which completely agrees with $\tilde{\tau}_{\alpha}$ as in non-interacting case. The all-left state has the maximum expansion coefficients for all eigenmodes, $|C_{\alpha}|=\Psi_{\alpha}$.

It turns out that the eigenmode $\alpha^{*}$ that gives the largest relaxation time $\tau_{\mathrm{max}}=\max_{\alpha}\tilde{\tau}_{\alpha}=\tilde{\tau}_{\alpha^{*}}$ appears in the middle of the spectrum, $\Lambda_{\alpha^{*}}\propto N$, not at the spectrum edge. Since the relaxation particularly slows down far from equilibrium with large $N_{+}$, an extensive number of particles simultaneously undergo dissipation in a many-body metastable state, which would be the reason why eigenmodes with extensively large eigenvalues are important here.

We plot $\tau_{\mathrm{max}}$ against $g$ in the inset of Fig. 3. We find that $\tau_{\mathrm{max}}$ increases exponentially in $g$, which is consistent with the intuitive picture that interactions effectively increase the energy barrier. This increase of the relaxation time stems from rapid growth of expansion coefficients with $g$.

In order to compare the physical relaxation time with $\tau_{\mathrm{max}}$, we numerically compute the dynamics of $\braket{N_{+}}=\sum_{N_{+}}N_{+}P_{N_{+}}(t)$ starting from the all-left initial state. The numerical result for different $g$ at $N=30$ is shown in Fig. 4. The quantity $\braket{N_{+}}$ reaches the stationary value roughly at $t=\tau_{\mathrm{max}}$, which is shown by vertical dashed lines in Fig. 4. It is thus confirmed that $\tau_{\mathrm{max}}$ coincides with the physical relaxation time of the system.

In the model considered so far, interaction-induced metastability is explained by the many-body explosion of expansion coefficients. Here we examine the FA model Fredrickson and Andersen (1984) as a prototypical example of kinetically constrained models Ritort and Sollich (2003) studied in glass and jamming transitions Berthier and Biroli (2011). Metastability in a kinetically constrained model is not due to energy barriers but due to dynamical rules. As we see below, the many-body explosion of expansion coefficients also explains this kind of metastability.

Suppose that each site $i\in\{1,2,\dots,L\}$ is occupied ($n_{i}=1$) or empty ($n_{i}=0$), and the transition rates from $n_{i}=0$ to 1 and that from $n_{i}=1$ to 0 are given by $C_{i}\rho$ and $C_{i}(1-\rho)$, respectively, where $\rho$ corresponds to the average density at equilibrium and $C_{i}$ expresses dynamical constraints. In the FA model Fredrickson and Andersen (1984), $C_{i}=1-n_{i-1}n_{i+1}$, which means that the site $i$ can change its state only if at least one of the two neighbors $i-1$ and $i+1$ is empty. Since this transition rate satisfies the detailed balance condition with respect to the non-interacting Hamiltonian $H=-\mu\sum_{i=1}^{L}n_{i}$, where $\mu=\ln[\rho/(1-\rho)]$ denotes the chemical potential, the equilibrium state is trivial but its dynamics exhibits interesting properties, which is a common characteristic of kinetically constrained models.

In the FA model, it is rigorously proved that the spectral gap of $G$ remains finite in the thermodynamic limit for any $0\leq\rho<1$ Cancrini et al. (2008). However, the maximum relaxation time exhibited by this model diverges in the thermodynamic limit for any $\rho$. Because of the kinetic constraint, a cluster of particles can decay or grow only from its edges. Therefore, if the initial state has a macroscopically large cluster of size $O(L)$, the relaxation time is at least of $O(L)$. A state with macroscopically large clusters is therefore regarded as a metastable state. This increase of the maximum relaxation time should be due to an exponential growth of $\Psi_{\alpha}$ ³³3At $\rho=1/2$, the Hamiltonian vanishes and $G$ is a symmetric matrix. Even in this case, the many-body explosion of expansion coefficients occurs..

Under the periodic boundary condition, $\{\tilde{\tau}_{\alpha}\}$ for $\rho=0.1$ are computed and plotted in Fig. 5. We see that in the FA model $\Psi_{\alpha}=e^{O(L)}$ for $\alpha$ with $\Lambda_{\alpha}=O(1)$, which results in $\tau_{\mathrm{max}}\propto L$ even though the spectral gap remains finite.