Spectral asymptotics for metastable processes in temperature-dependent traps

The timescale barrier in MD simulations

• Long trajectories are useful to estimate dynamical properties of atomistic systems, and parametrize models on larger scales.
• Length of sampled trajectories is hindered by metastability, the presence of phenomena occuring on widely separated timescales.
• This is due to the presence of free-energetic barriers much higher than \(k T\).
• Accelerated dynamics methods ( Hyperdynamics (Arthur F. Voter (1997)), TAD (Sørensen and Voter (2000))) can help.
• Parallelization in space is “easy”, parallelization in time is much harder. Picture from Danny Perez (LANL).
• We focus on Parallel Replica (A. F. Voter (1998)) which adresses this issue.

The parallel replica algorithm

• Exploits local metastability inside \(\Omega\subset \mathcal X\) to accelerate sampling of the transition (= exit event from \(\Omega\)).
• Idea: if dynamics \(X_t\) stays in \(\Omega\) for a time \(\gg 1\), a local equilibrium is reached, and the way \(X_t\) exits \(\Omega\) becomes independent of how or when it entered.
• This allows to sample the exit event \((T_{\mathrm{exit}}(\Omega),X_{T_{\mathrm{exit}}(\Omega)})\) in parallel, conditionally on reaching local equilibrium.
• Inputs to ParRep: a set of metastable states \((\Omega_i)_{i\geqslant 1}\), and associated decorrelation timescales \(T_{\mathrm{corr}}(\Omega_i)\).

ParRep algorithm

Decorrelation step: when the dynamics \(X_t\) enters a metastable state \(\Omega\) at time \(t=0\), evolve during \(T_{\mathrm{corr}}(\Omega)\).

• If \(X_t\in\Omega\) for all \(0\leqslant t \leqslant T_{\mathrm{corr}}(\Omega)\), \(X_{T_{\mathrm{corr}}}(\Omega)\) is (almost) in local equilibrium \(\nu_\Omega \in \mathcal M_1(\Omega)\).
• If \(X_t\) exits at a time \(T_{\mathrm{exit}}(\Omega) < T_{\mathrm{corr}}(\Omega)\), record the (correlated) exit event and start again.

ParRep algorithm

Dephasing step: spawn \(N\) independent replicas, and evolve for \(T_{\mathrm{corr}}(\Omega)\), in parallel.

• During this step, the physical clock is paused.
• The replicas are asymptotically independent in the limit \(T_{\mathrm{corr}}(\Omega)\to\infty\).

ParRep algorithm

Rejection/conditioning step: discard replicas which exited the state in the dephasing step.

• The remaining \(N_{\mathrm{par}}\) replicas (almost) i.i.d. in local equilibrium.
• More sophisticated variants exist (e.g. Flemming-Viot processes).

ParRep algorithm

Parallel sampling step: let \(X^1_{T_{\mathrm{par}}}\) be the state of the first replica to exit \(\Omega\), after an additional time \(T_{\mathrm{par}}\).

• Conditionally on \(\{X_t \in \Omega\}\) for \(t< T_{\mathrm{corr}}(\Omega)\), the pair \(\left( T_{\mathrm{decorr}} + N_{\mathrm{par}}T_{\mathrm{par}},X^1_{T_{\mathrm{par}}}\right)\) is an (almost) unbiased sample of the exit event.

ParRep efficiency

• The algorithm is efficient if the transition timescale \(\mathbb E_{\nu_\Omega}[T_{\mathrm{exit}}(\Omega)]\) is much larger than the decorrelation timescale \(T_{\mathrm{corr}}(\Omega)\).
• Speedup in wallclock time related to separation of timescales \(\delta(\Omega) = T_{\mathrm{corr}}(\Omega)/\mathbb{E}_{\nu_\Omega}[T_{\mathrm{exit}}(\Omega)]\) and number of available replicas \(N\).
• Standard approach: \(\Omega=\) basin of attraction for steepest descent of the potential, \(T_{\mathrm{corr}}(\Omega)\) given by harmonic approximation around the local minimum.
• In general, this is not optimal nor valid (as we show theoretically).

Mathematical setting: the quasi-stationary distribution (QSD)

• For theoretical purposes, we consider the overdamped Langevin dynamics \[ \mathrm{d}X_t^\beta = -\nabla V\left(X_t^\beta\right)\,\mathrm{d}t + \sqrt{\frac{2}{\beta}}\,\mathrm{d}W_t. \]
• A quasi-stationary distribution in \(\Omega\) is a probability measure \(\nu_\Omega\in \mathcal M_1(\Omega)\) satisfying \[ \forall\,A\subset \Omega,\qquad \mathbb P_{\nu_\Omega}\left(X_t \in A\,\middle|\,t<T_{\mathrm{exit}}(\Omega)\right) = \nu_\Omega(A). \]
• Under mild conditions on \(V\) and \(\Omega\), the QSD is unique, and related to the principal Dirichlet eigenspace. More precisely, \[ \nu_\Omega(\mathrm{d}x) \propto u_{1,\Omega}(x)\mathrm{e}^{-\beta V(x)}\,\mathrm{d}x,\qquad -\mathcal{L}_\beta u_{1,\Omega}(x) = \lambda_{1,\Omega}u_{1,\Omega}(x), \] where \(\mathcal{L}_\beta\) is the infinitesimal generator of the dynamics with absorbing Dirichlet boundary conditions on \(\partial \Omega\): \[ \mathcal{L}_\beta = -\nabla V\cdot \nabla + \frac1\beta\Delta,\qquad \mathcal D(\mathcal{L}_\beta) = H_0^1 \cap H^2(\Omega; \mathrm{e}^{-\beta V(x)}\,\mathrm{d}x), \] which is self-adjoint with pure point spectrum on \(L^2(\Omega;\mathrm{e}^{-\beta V(x)}\,\mathrm{d}x)\) for bounded \(\Omega\).

Spectral characterization of timescales

• (Theorem ( Le Bris et al. 2012)):

  Under initial distribution \(\nu_\Omega\), \[ T_{\mathrm{exit}}(\Omega)\sim \mathcal E(\color{blue}{\lambda_{1,\Omega}}),\qquad T_{\mathrm{exit}}(\Omega) \perp\!\!\!\!\perp X_{T_{\mathrm{exit}(\Omega)}}. \]

• This justifies the unbiasedness of parallel sampling for the exit event starting from the QSD.

• The exit rate is given by \(\lambda_{1,\Omega}\).

• Denote \(\mu_t = \mathrm{Law}\left(X_t \middle| t<T_{\mathrm{exit}}(\Omega)\right)\) the law of the process conditioned on remaining trapped.

• (Theorem ( Le Bris et al. 2012)):

  Convergence to the QSD and bias on the exit event are exponentially decaying. There exist \(C_1,C_2>0\) independent of \(\mu_0\) such that \[ \left\|\mu_t-\nu_\Omega\right\|_{\mathrm{TV}} \leqslant C_1 \mathrm{e}^{-(\color{red}{\lambda_{2,\Omega}-\lambda_{1,\Omega}})t}, \] \[ \left\|\mathrm{Law}_{\mu_t}(T_{\mathrm{exit}(\Omega)}-t,X_{T_{\mathrm{exit}}(\Omega)}) - \mathrm{Law}_{\nu_\Omega}(T_{\mathrm{exit}(\Omega)},X_{T_{\mathrm{exit}}(\Omega)})\right\|_{\mathrm{TV}}, \leqslant C_2 \mathrm{e}^{-(\color{red}{\lambda_{2,\Omega}-\lambda_{1,\Omega}})t}, \] and the rate is given by the spectral gap \(\lambda_{2,\Omega} - \lambda_{1,\Omega}\), of the Dirichlet generator.

• The decorrelation rate is given by \(\lambda_{2,\Omega}-\lambda_{1,\Omega}\)

• This justifies choosing \(T_{\mathrm{corr}}(\Omega) = n_{\mathrm{corr}}/(\lambda_{2,\Omega}-\lambda_{1,\Omega})\) as a valid decorrelation time, for some fixed tolerance parameter \(n_{\mathrm{corr}}\).

A spectral shape optimization problem

• The goal becomes to maximize the efficiency of ParRep by maximizing the separation of timescales: \[ \underset{\Omega \subset \mathcal X}{\max}\,J(\Omega),\qquad J(\Omega):= \frac{\color{red}{\lambda_{2,\Omega}-\lambda_{1,\Omega}}}{\color{blue}{\lambda_{1,\Omega}}}. \]
• In other words, make the exit rate as small as possible compared to the decorrelation rate, to make \(\Omega\) as locally metastable as possible.
• Direct numerical approaches using shape optimization of eigenvalues and projection to collective variable space can be implemented. Ongoing work in this direction.
• Here, we take a theoretical and asymptotic approach, estimating the Dirichlet spectrum in the low-temperature regime \(\beta\to\infty\) for a temperature-dependent domain \(\Omega_\beta\).
• To get tractable estimates, the shape of \(\Omega_\beta\) is constrained by the value \(\alpha\) of a specific low-dimensional reaction coordinate.
• Goal: derive sharp, quantitative estimates, of the rates \(\lambda_{1,\Omega_\beta},\lambda_{2,\Omega_\beta}\) governing transitions and decorrelation, as a function of \(\alpha\).
• Then, it becomes easy to optimize the asymptotic behavior of \(J(\Omega_\beta)\) with respect to \(\alpha\).

Main geometric assumption

• \(V\) is a Morse function (non-degenerate Hessians at critical points), and \(\Omega_\beta \subset \mathcal K\) with \(\mathrm{Crit} V|_{\mathcal K}=\{z_0,\dots,z_{N-1}\}\).
• Assume the following limit is well-defined in \((-\infty, +\infty]\) for each \(0\leqslant i < N\). \[ \alpha_i = \underset{\beta\to\infty}{\lim}\, \sqrt\beta\sigma_{\Omega_\beta}(z_i),\qquad \sigma_{\Omega_\beta}(x) = \begin{cases} d(x,\partial\Omega_\beta),& x\in \Omega_\beta,\\ -d(x,\partial\Omega_\beta)&x\in\mathcal X\setminus \Omega_\beta \end{cases}. \]
• If \(\alpha_i = +\infty\), \(z_i\) is far from the boundary, and if not, it is close.
• Assume there exists a function \(\delta :{\mathbb{{R}}}_+\to (0,\varepsilon(V))\) such that \[ \sqrt\beta\delta(\beta)\xrightarrow{\beta\to\infty}\,+\infty,\qquad\mathcal{O}_{i}^{-}(\beta) \subseteq B(z_i,\delta(\beta))\cap\Omega_\beta \subseteq \mathcal{O}_{i}^{+}(\beta), \] where the sets \(\mathcal{O}_{i}^{\pm}\) constrain the geometry of the boundary near the \(z_i\) in a specific way.

Main geometric assumption

• In the drawing above, \(\mathcal{O}_{i}^{\pm}\) are the two capped balls, \(\sqrt\beta\gamma(\beta)\to 0\), and \(v_1^{(i)}\) is an eigenvector of \(\nabla^2 V(z_i)\).
• If \(z_i\) is an index 1 saddle points, we assume \(v_1^{(i)}\) is the unique eigenvector for the negative eigenvalue.
• Boundary near a saddle \(z_i \approx\) level set of the position along the minimum energy path through the saddle.
• The standard case \(\Omega_\beta = \Omega\) satisfying some transversality condition is included.

First-order asymptotics: harmonic approximation of the full spectrum

• Theorem (Blassel, Lelièvre, and Stoltz (2024), in preparation):
  

  For any \(k\), it holds \[ \underset{\beta\to\infty}{\lim}\, \lambda_{k,\Omega_\beta} = \lambda_k^{\mathrm H}, \] where \(\lambda_k^{\mathrm H}\) is the \(k\)-th (temperature independent) eigenvalue of a direct sum of local harmonic models \[ H^\alpha_\beta = \bigoplus_{i=0}^{N-1} H_\beta^{(i),\alpha_i},\qquad H_\beta^{(i),\alpha_i} = \beta(x-z_i)^\top\nabla V^2(z_i)(x-z_i)/4-\Delta V(z_i)/2 - \frac1\beta \Delta, \] acting on the half-space \(\left\{ (x-z_i)^\top v_1^{(i)} < \alpha^{(i)}/\beta\right\}\) with Dirichlet boundary conditions.

  • Example: Assume only one minimum \(z_0\) and index 1 saddle points \(z_1,\dots,z_{N-1}\). \[ \lambda_{1,\Omega_\beta}\xrightarrow{\beta\to\infty}0,\quad \lambda_{2,\Omega_\beta}\xrightarrow{\beta\to\infty} \min\left[\nu^{(0)}_{1},\underset{0<i<N}{\min}|\nu^{(i)}_{1}|\left(\mu_{0,\alpha^{(i)}\sqrt{|\nu^{(i)}_{1}|/2}}+\frac12\right)\right], \] where \(\mu_{0,\theta}=\)ground-state energy of \(\frac12(x^2-\partial_x^2)\) acting on \((-\infty,\theta)\) with Dirichlet boundary conditions.
    Here, \(\nu^{(i)}_{k}=\) \(k\)-th eigenvalue of \(\nabla^2 V(z_i)\).
  • Result provides quantitative estimate for the decorrelation rate \(\lambda_{2,\Omega_\beta}-\lambda_{1,\Omega_\beta}\), and shows bottom-of-the-well approximation cannot be valid in general.
  • Proof adapts a method of Simon (1983) from semiclassical analysis to our setting of metastable diffusions with temperature-dependent boundary.

Additional assumptions

• We now restrict the setting to the case of only one minimum far from the boundary \(z_0\), and any other critical points.
• Harmonic approximation implies \(\lambda_{1,\Omega_\beta}\to 0\), \(\lambda_{2,\Omega_\beta}\to c(\alpha)>0\).
• Define \(V^* = \min\left\{V(z),\,z\in\partial \left[{\mathbb{{R}}}^d\setminus \mathcal A(z_0)\right]\right\}\), where \(\mathcal A(z_0)=\) basin of attraction of \(z_0\) for \(X'(t)=-\nabla V(X(t))\).
• The minimizers are index one saddle points \(\left\{z_i,i\in I_{\min}\right\}\).
• Assume \(\mathcal A(z_0)\cap \{V<V^*\} \subset \Omega_\beta\) for all \(\beta\).

Modified Eyring-Kramers formula

• Theorem (Blassel, Lelièvre, and Stoltz (2024), in preparation):
  

  In the limit \(\beta\to\infty\), it holds \[ \lambda_{1,\Omega_\beta} = \mathrm{e}^{-\beta(V^*-V(z_0))} \left[\sum_{i\in I_{\min}} \frac{|\nu^{(i)}_{1}|}{2\pi\color{blue}{\Phi\left(|\nu^{(i)}_{1}|^{\frac12}\alpha^{(i)}\right)}}\sqrt{\frac{\det \nabla^2 V(z_0)}{\left|\det \nabla ^2 V(z_i)\right|}} \left(1 +\mathcal{O}(1/\sqrt\beta)\right)\right], \] where \[ \Phi(x) = \frac1{\sqrt{2\pi}}\int_{-\infty}^x \mathrm{e}^{-\frac{t^2}2}\,\mathrm{d}t. \]

• Result provides sharp quantitative estimate for the metastable exit rate \(\lambda_{1,\Omega_\beta}\) Dependence on \(\alpha\) is only in the prefactor \(\color{blue}{P}\). Generalizes known Eyring-Kramers (\(\alpha=+\infty\), \(P=1\)) and (\(\alpha=0\), \(P=1/2\)) formulae.
• Can also be understood as providing finer description (compared to previously known results) of the transition in the first eigenvalue as the boundary of a domain crosses a low-energy saddle point.
• Proof relies on the construction of a very precise quasimode for \(u_{1,\Omega_\beta}\), inspired by works in probability by Bovier, Gayrard, and Klein (2005) and semiclassical analysis of Witten Laplacians by Le Peutrec and Nectoux (2021).

Conclusion and outlook

Transition for the prefactor \(\mathrm{e}^{-\beta (V^*-V(z_0))} (\lambda_{2,\Omega_\beta}-\lambda_{1,\Omega_\beta})/\lambda_{1,\Omega_\beta}\) in the case of a single saddle point.

• Our results provide sharp estimates of crucial timescales for ParRep in the low-temperature regime.
• Lead to practical prescriptions on how to define better metastable states and decorrelation times.
• However, limited to setting of energetic barriers. More work needed to take entropic effects into account.
• Preprint soon to appear on ar\(\chi\)iv.