On Ergodicity for Multidimensional Harmonic Oscillator Systems with Nosé – Hoover-type Thermostat

Abstract

A simple proof and a detailed analysis of the nonergodicity for multidimensional harmonic oscillator systems with the Nosé – Hoover-type thermostat are presented. The origin of the nonergodicity is the symmetries in the multidimensional target physical system, and it differs from that in the Nosé – Hoover thermostat with the one-dimensional harmonic oscillator. A new and simple deterministic method to recover the ergodicity is also proposed. An individual thermostat variable is attached to each degree of freedom, and all of these variables act on a friction coefficient for each degree of freedom. This action is linear and controlled by a Nosé mass matrix $\mathbf{Q}$, which is the matrix analogue of the scalar Nosé mass. The matrix $\mathbf{Q}$ can break the symmetry and contribute to the ergodicity.

Appendices

APPENDIX A

We will state that a linear symmetry $S:\mathbb{R}^{n}\rightarrow\mathbb{R}^{n}$ acts on the ODE (2.1) if it satisfies the following:

Definition 1

$S\in\mathrm{End}\mathbb{R}^{n}$ preserves the functions $\lambda$ and $\Lambda$ as well as the domain $D$, such that $\lambda\left(S(x),S(p),\zeta\right)=\lambda\left(x,p,\zeta\right)$ and $\Lambda\left(S(x),S(p),\zeta\right)=\Lambda\left(x,p,\zeta\right)$ hold for all $\left(x,p,\zeta\right)\in\Omega$ and $S(D)\subset D$. Commutativities also hold: $\mathbf{M}^{-1}\circ S=S\circ\mathbf{M}^{-1}$ and $F\circ S=S\circ F$. 

We denote by $\mathcal{S}$ the set of all $S\in\mathrm{End}\mathbb{R}^{n}$ that act on ODE (2.1).

Lemma 2

For any $S\in\mathcal{S}$ , we have $S(x(t))=x(t)$ and $S(p(t))=p(t)$  for all $t$ in an interval $J\subset\mathbb{R}$ if $\varphi:J\rightarrow\Omega,t\mapsto\left(x(t),p(t),\zeta(t)\right)$ is a solution to ODE (2.1) with an initial condition satisfying $S\big{(}x(0)\big{)}=x(0)$ and $S\big{(}p(0)\big{)}=p(0)$ .

Proof

It follows from Definition 1 that $\hat{\varphi}:\mathbb{R\supset}J\rightarrow\Omega,t\overset{\mathrm{d}}{\mapsto}\left(S(x(t)),S(p(t)),\zeta(t)\right)$ also becomes a solution to the $C^{1}$ ODE and $\hat{\varphi}(0)=\varphi(0)$ holds. Thus, the uniqueness of the initial value problem ensures that $\hat{\varphi}=\varphi$, so that $S(x(t))=x(t)$ and $S(p(t))=p(t)$ for all $t\in J$.

Thus, for every $S\in\mathcal{S}$, we have an invariant set

$$\Omega_{S}:=\Gamma_{s}^{D}\times\Gamma_{s}\times\mathbb{R}^{m},$$

with

$$\Gamma_{s}\equiv\{p\in\mathbb{R}^{n}|S(p)=p\}$$

and $\Gamma_{s}^{D}\equiv\Gamma_{s}\cap D\equiv\{x\in D$ $|$ $S(x)=x\}$, indicating that the symmetry $\big{(}S(x),S(p)\big{)}=(x,p)$ is maintained in the dynamics or is compatible with the ODE. To analyze the nonergodic condition (3.5), we use the approach of finding an invariant set $A$, which either meets condition (3.5) itself or separates the total phase space into three invariant sets,

$$\Omega=A\sqcup A^{[+]}\sqcup A^{[-]},$$

(A.1)

where $A^{[+]}$ satisfies condition (3.5). In both cases, $A$ should be “large” (when selecting $A=\Omega_{S}$, there is an extremely small case of $\Gamma_{s}=\{0\}$ or a case of a low-dimensional subspace). Thus, our target for $A$ is an invariant set that is these $\Omega_{S}$ summed in a certain manner:

$$A_{\mathcal{R}}\equiv\bigcup_{S\in\mathcal{R}}\Omega_{S}=\bigcup_{S\in\mathcal{R}}(\Gamma_{s}^{D}\times\Gamma_{s})\times\mathbb{R}^{m},$$

(A.2)

where $\mathcal{R}$ is a certain subset of $\mathcal{S}$ such that it is sufficiently but not excessively large. For example, for the latter condition, we should remove the case in which $S$ is the identity $\mathrm{id}_{\mathbb{R}^{n}}$; otherwise, $A_{\mathcal{R}}$ becomes “too large” ($S=\mathrm{id}_{\mathbb{R}^{n}}$ provides $\Omega_{\mathrm{id}_{\mathbb{R}^{n}}}=\Omega$, and thus yields $\Omega\backslash A_{\mathcal{R}}=\varnothing$, which does not contribute to the nonergodic condition (3.5) for $A\equiv A_{\mathcal{R}}$).  

We will show that (A.1) holds with $A^{[\pm]}\equiv\Upsilon_{ij}^{\pm}$ if $A\equiv A_{\mathcal{R}}$, in a special case of the isotropic $n$HO. Here $\Upsilon_{ij}^{\pm}$ are defined in Lemma 1, and this fact can explain the route whereby $\Upsilon_{ij}^{\pm}$ arise. That is, they arise as a complementary set to a sum, within a certain range $\mathcal{R}$, of the invariant set $\Omega_{S}$ based on the symmetry $S$ that acts on the ODE. To illustrate the issue, we restrict the condition such that the dependence of $x,p$ in the functions $\lambda$ and $\Lambda$ is only through the potential and kinetic energies, as follows:

Condition 1

There exist $C^{1}$ functions $\tilde{\lambda},\tilde{\Lambda}\colon\mathbb{R}\!\times\!\mathbb{R}\!\times\!\mathbb{R}^{m}\!\rightarrow\!\mathbb{R}\$such that$\ \lambda\left(x,p,\zeta\right)\!=\!\tilde{\lambda}\left(U(x),K(p),\zeta\right)$ and $\Lambda\left(x,p,\zeta\right)\!=\!\tilde{\Lambda}\left(U(x),K(p),\zeta\right)$ for all $\left(x,p,\zeta\right)\in\Omega$.

This is not a special condition and it is satisfied by Examples 1–4. For the HO system described by Eq. (3.6) under Condition 1, orthogonal transforms of $\mathbb{R}^{n}$ actually act on the ODE (2.1):

Lemma 3

$\mathcal{S\supset}O(n)\equiv\{S\in\mathrm{End}\mathbb{R}^{n}$ $|$ ${}^{\mathrm{T}}SS=\mathrm{id}_{\mathbb{R}^{n}}\}$ holds for the HO system.

Proof

Consider any $S\in O(n)$. Because $U\big{(}S(x)\big{)}=U(x)$ and $K\big{(}S(p)\big{)}=K(p)$ hold for any $(x,p)$, and because $\mathbf{M}^{-1}$ and $F$ become diagonal, the conditions in Definition 1 are valid, indicating that $O(n)\subset\mathcal{S}$.

We then put

$$\mathcal{R}\equiv O(n)\backslash\{\mathrm{id}_{\mathbb{R}^{n}}\}=:O(n)_{\times},$$

namely, we sum $\Omega_{S}$ to obtain $A_{\mathcal{R}}$ for all $S$ that are non-identical orthogonal transforms of $\mathbb{R}^{n}$. We also restrict our consideration to $n=2$ for ease, wherein the discussion would be extended to a larger $n$. The following proposition demonstrates how $\Upsilon_{12}^{\pm}$ arises from $\mathcal{R}$.

Proposition 2

The equation $A_{\mathcal{R}}=\Upsilon_{12}^{0}$ holds, and decomposition (A.1) holds for $A=A_{\mathcal{R}}$ and ${A^{[\pm]}=\Upsilon_{12}^{\pm}}$ .

Proof

By using the fact that $O(2)$ is bijectively parameterized by $S^{1}$ and the signature

$$S^{1}\times\{\pm 1\}\rightarrow O(2),\text{ }(\theta,\sigma)\mapsto\left[\begin{array}[c]{cc}\cos\theta&-\sigma\sin\theta\\ \sin\theta&\sigma\cos\theta\end{array}\right]=:s_{\theta}^{\sigma},$$

and by solving an eigenvalue problem $s_{\theta}^{\sigma}(p)=p$, we observe that

$$\displaystyle\Gamma_{s_{\theta}^{+1}}=\left\{\begin{array}[c]{cc}\mathbb{R}^{2}&\text{if }\theta=0\\ \{0,0\}&\text{otherwise}\end{array}\right\},$$

$$\displaystyle\Gamma_{s_{\theta}^{-1}}=\left\{\begin{array}[c]{lc}\mathbb{R}\times\{0\}&\text{if }\theta=0\\ \{0\}\times\mathbb{R}&\text{if }\theta=\pi\\ \left\{(p_{1},p_{2})\in\mathbb{R}^{2}\text{ }|\text{ }p_{2}=\frac{1-\cos\theta}{\sin\theta}p_{1}\right\}&\text{otherwise}\end{array}\right\}.$$

Note that $\Gamma_{s_{\theta}^{-1}}$ is a line through the origin of $\mathbb{R}^{2}$ with a gradient $k=k_{\theta}$ for $\theta\in S^{1}\backslash\{0,\pi\}$, where $k_{\theta}$ can take any value in $\mathbb{R}_{\times}$, and that $\Gamma_{s_{0}^{-1}}=\mathbb{R}\times\{0\}$ and $\Gamma_{s_{\pi}^{-1}}=\{0\}\times\mathbb{R}$ are also lines with gradients of $0$ and $\infty$, respectively. Thus $\Gamma_{s_{\theta}^{-1}}=L_{k}$, a line through the origin of $\mathbb{R}^{2}$ with a gradient $k\in(\mathbb{-\infty},\mathbb{\infty]}$. Hence, by applying $\mathcal{R=}O(2)_{\times}=\{s_{\theta}^{+1}$ $|$ $\theta\in S_{\times}^{1}\}\cup\{s_{\theta}^{-1}$ $|$ $\theta\in S^{1}\}$, we get

$$\displaystyle\bigcup_{S\in\mathcal{R}}(\Gamma_{s}^{D}\times\Gamma_{s})$$

$$\displaystyle=\bigcup_{\theta\in S_{\times}^{1}}(\Gamma_{s_{\theta}^{+1}}\times\Gamma_{s_{\theta}^{+1}})\cup\bigcup_{\theta\in S^{1}}(\Gamma_{s_{\theta}^{-1}}\times\Gamma_{s_{\theta}^{-1}})$$

$$\displaystyle=\bigcup_{\theta\in S^{1}}(\Gamma_{s_{\theta}^{-1}}\times\Gamma_{s_{\theta}^{-1}})$$

$$\displaystyle=\bigcup_{k\in(\mathbb{-\infty},\mathbb{\infty]}}(L_{k}\times L_{k})$$

$$\displaystyle=\{(x_{1},kx_{1},p_{1},kp_{1})|\text{ }x_{1},p_{1},k\in\mathbb{R}\}$$

$$\displaystyle\cup\text{ }\{0\}\times\mathbb{R}\times\{0\}\times\mathbb{R}$$

$$\displaystyle=\{(x_{1},x_{2},p_{1},p_{2})\in\mathbb{R}^{4}\text{ }|\text{ }x_{1}p_{2}-x_{2}p_{1}=0\}.$$

Thus, $A_{\mathcal{R}}=\bigcup_{S\in\mathcal{R}}(\Gamma_{s}^{D}\times\Gamma_{s})\times\mathbb{R}^{m}=\{(x,p)\in\mathbb{R}^{4}|$ $x_{1}p_{2}-x_{2}p_{1}=0\}\times\mathbb{R}^{m}=\{\omega\in\Omega$ $|$ $\gamma_{12}(\omega)=0\}=\Upsilon_{12}^{0}$. Therefore, decomposition (A.1) holds for $\Omega=\Upsilon_{12}^{0}\sqcup\Upsilon_{12}^{+}\sqcup\Upsilon_{12}^{-}$.

Note that the explicit form of $\Gamma_{s}$ in the proof directly indicates nontrivial examples to explain that $\mathcal{R}$ should be sufficiently large. For example, if we take $\mathcal{R}$ as a one-point set, $A_{\mathcal{R}}=\bigcup_{S\in\mathcal{R}}(\Gamma_{s}^{D}\times\Gamma_{s})\times\mathbb{R}^{m}$ does not separate the phase space into three spaces as described in (A.1): if we set $\mathcal{R}=\{s_{\theta}^{+1}\}$ with $\theta\neq 0$ or $\mathcal{R}=\{s_{\theta}^{-1}\}$, then $\bigcup_{S\in\mathcal{R}}(\Gamma_{s}^{D}\times\Gamma_{s})=\Gamma_{s_{\theta}^{+1}}\times\Gamma_{s_{\theta}^{+1}}=\{0,0,0,0\}\subset\mathbb{R}^{4}$ or $\bigcup_{S\in\mathcal{R}}(\Gamma_{s}^{D}\times\Gamma_{s})=\Gamma_{s_{\theta}^{-1}}\times\Gamma_{s_{\theta}^{-1}}=$line$\times$line$\subset\mathbb{R}^{4}$, respectively, clearly induces no separation. The same applies to, for example, any finite set $\mathcal{R}=\{s_{\theta_{1}}^{\pm 1},\ldots,s_{\theta_{q}}^{\pm 1}\}$.

APPENDIX B

We also examined the kinetic moments (KM) method (Example 3) as a conventional method, with the same simulation conditions, including the masses of $Q_{1}=Q_{2}=1$, as those used for the NHC method (Example 2), except for the use of a smaller unit time of $h=2\times 10^{-4}$ and longer time steps of $5\times 10^{8}$. This is because the feature of the higher-order terms with respect to $p$, namely, $\hat{K}(p)p$ and $\hat{K}(p)\hat{K}(p)$, makes the EOM stiff. The explicit second-order scheme described in Ref. [10] with special care of the local refinement [24, 15], along with the small unit time, made the conservation of the invariant function value comparable with that in the NHC method. Four initial conditions, namely, Cases (A), (A’), (B), and (C) were considered.

Whereas Case (A) maintained the interchange symmetry and trajectory confinement in the invariant set $A=A_{\mathcal{R}}$ (which means that $x_{1}(t)=x_{2}(t)$, $p_{1}(t)=p_{2}(t)$ and $\gamma_{12}(\omega(t))=0$, respectively, hold for all times $t$), the trajectory of $(x_{1},p_{1})$ and distribution of $x_{1}$ implied nonergodic behavior, as illustrated in Fig. 7. The trajectory of $(\zeta_{1},\zeta_{2})$ also exhibited a deviation from the theoretical Gaussian distribution $\varpropto\exp\left[-\sum_{i=1}^{2}\zeta_{i}^{2}/Q_{j}\right]d\zeta$. For Case (A’), the trajectories and distributions exhibited similar nonergodic behavior (not shown). These features are essentially the same as those obtained in the NHC method. The result of Case (B) differed from that of the NHC method (Fig. 2), where the KM method yielded relatively good results for the trajectories and distributions (not shown). A key feature is $\gamma_{12}$, which yielded the expected result on the signature ($\gamma_{12}\big{(}\omega(t)\big{)}>0$ for all $t$) as in the case of the NHC, but generated fluctuations that were one order of magnitude greater on $\gamma_{12}\big{(}\omega(t)\big{)}$ than those of the NHC (Fig. 3b). This would be the reason for the good results, as stated in the text, such that the trajectories that can be distant from $A$ exhibit relatively good sampling. A more detailed and analytical aspect of the reason is discussed in Appendix C in terms of numerical error-induced ergodic-like behavior. In Case (C), although unignorable errors were observed in the distributions, practically sufficient results were obtained, as in the NHC method. This is probably because the magnitude of $\gamma_{12}\big{(}\omega(t)\big{)}$ was comparable to that of the NHC. The stiffness in the EOM of the KM method should contribute to moving away from the invariant set $A$ and facilitate the numerical error-induced ergodic-like behavior.

Fig. 7

Simulation results obtained by the conventional thermostat method (KM method) for $2$HO using the initial condition $x_{1}(0)=x_{2}(0)=0,$ $p_{1}(0)=p_{2}(0)=1$: (a) trajectory (scatter plot) of $(x_{1},p_{1})$ and marginal distributions for (b) $x_{1}$ (with theoretical distribution and discrepancies) and (c) trajectory of $(\zeta_{1},\zeta_{2})$.

APPENDIX C

For a further study of the $n=2$ isotropic HO system, we also implemented an initial condition in the NHC method, that is, $x_{1}(0)=x_{2}(0)=0,$ $p_{1}(0)=1+10^{-9}$, and $p_{2}(0)=1$, which is obtained by adding a small perturbation to $p_{1}(0)$ in Case (A) described in the text. All other simulation conditions were the same as those of Case (A). Note that this initial condition, or the condition that $x(0)=0$ in general, does not lose the symmetry, but rather retains the symmetry. This is because for any $p(0)\in\mathbb{R}^{2}$ there exists $\theta\in S^{1}$ such that $s_{\theta}^{-1}p(0)=p(0)$, where $S=s_{\theta}^{-1}=\left[\begin{array}[c]{cc}\cos\theta&\sin\theta\\ \sin\theta&-\cos\theta\end{array}\right]$ (see the proof of Proposition 2).

The distribution of $x_{1}$ obtained in the simulation is presented in Fig. 8c. This distribution agrees well with the theoretical distribution, and thus, the system appears ergodic, which is in contrast to the nonergodic distribution of Case (A) depicted in Fig. 1b. This contrast may appear contradictory, considering our theoretical conclusion that the isotropic HO system with $n>1$ thermostatted by the EOM (2.1) is nonergodic (Proposition 1). However, this good distribution obtained in this case should be owing to ergodic-like behavior derived from numerical errors, which we refer to as numerical error-induced ergodic-like behaviour [11] (which should be essentially the same context as that discussed in Ref. [5]). Namely, although the system is intrinsically nonergodic, it may often appear ergodic from a certain aspect owing to the numerical errors, which is due to the numerical integrator (time-discretization scheme of the ODE) and would also be due to computational round-off errors.

Fig. 8

(a) Trajectories of $\gamma_{12}(\omega(t))$ for several values of unit time $h$ obtained by the NHC method with the initial condition $x_{1}(0)=x_{2}(0)=0$, $p_{1}(0)=1+10^{-9}$, and $p_{2}(0)=1$. The arrow indicates the time ($t_{\text{break}}$) when $|\gamma_{12}(\omega(t))|>\epsilon$ occurs, which shows the theoretical condition breaking, as measured by the threshold $\epsilon=0.1$. (b) Standard deviations of values of the invariant function $I$ during the total simulation time of $10^{5}$, which is common for all cases, with varying $h$. (c) Marginal distributions of $x_{1}$ obtained for total time $10^{5}$ with unit time step $h=10^{-3}$, and (d) that with $h=10^{-6}$. (e) Marginal distributions of $x_{1}$ obtained by an alternative integrator [33] with the total time of $10^{5}$ and $h=10^{-6}$. (f) Schematic figure of phase space $\Omega$, which is not arcwise connected: an invariant subset $A=A_{\mathcal{R}}$ separates $\Omega$ into the two positive-measure subsets $A^{[+]}$ and $A^{[-]}$. This indicates the nonergodicity, where the trajectories starting at $A$, $A^{[+]}$ and $A^{[-]}$ are confined therein, respectively, for the entire time; namely, transitions among distinct sets $A$, $A^{[+]}$, and $A^{[-]}$ are forbidden. However, numerical errors often induce transitions, such as that from $A$ to $A^{[+]}$, as indicated by the trajectory (green curve), which clearly appears in the behavior of $\gamma_{12}(\omega(t))$ in (a).

To confirm this, we investigated the system behaviors by reducing the integration unit time $h$, from the default value of $10^{-3}$ to the smallest value of $10^{-6}$. Figure 8a presents the values of $\gamma_{12}(\omega(t))$, which should theoretically be zero under this initial condition $\omega_{0}\in A$, as stated in the text. However, it deviates from zero after a certain time $t_{\text{break}}$, as indicated by the arrow, and thereafter, it moves significantly away from zero and begins to oscillate. With decreasing $h$, the time $t_{\text{break}}$ is delayed, but $\gamma_{12}\big{(}\omega(t)\big{)}$ is eventually away from zero and begins to oscillate, while the signatures are maintained (positive for $h=10^{-5}$, or negative for other $h$). The time delay corresponds to the fact that the numerical integration errors are smaller with decreasing $h$, which is explicitly exhibited in Fig. 8b in terms of the error of the invariant function defined on the extended space [16, 32] (where the total simulation time was $10^{5}$, which is common for every $h$ case). In fact, $\gamma_{12}(\omega(t))$ remains zero for the longest duration in the case of $h=10^{-6}$, and the marginal distribution shown in Fig. 8d exhibits unignorable errors around $x_{1}=0$. Namely, the distribution with $h=10^{-6}$, which should be the most accurate, is worse than the distribution with $h=10^{-3}$, which should be the most inaccurate. Thus, the distribution that deviates from the exact BG distribution, as illustrated in Fig. 8d, almost reflects the intrinsic property of the NHC EOM, and the distribution presented in Fig. 8c should be due to the numerical error-induced ergodic-like behavior, which is far from the feature of the exact EOM.

Even if strictly keeping $\gamma_{12}(\omega(t))=0$ would be difficult, a weaker condition represented by the long-time average, $\bar{\gamma}_{12}=0$, could be expected to hold. However, apparent nonzero averages were confirmed numerically such that $\gamma_{12}(\omega(t))$ was always positive or always negative for all cases shown in Fig. 8a, as well as Figs. 3b and 3c, and for cases for longer simulations (not shown). The numerical observations of $\bar{\gamma}_{12}\neq 0$ are reflected from the nonergodicity, as stated in the text. Nonzero averages after breaking the condition $\gamma_{12}=0$ were also observed when replacing the numerical integrator from the current one (the P2S1 integrator with the simplest evaluation for the thermostat part) to an alternative one [33] (which performs an enhanced evaluation for the thermostat part by combining a symmetrizing scheme, Yoshida’s order-increasing protocol, and simple iterations; we set $n_{\text{ys}}=n_{\text{c}}=3$), keeping other simulation conditions the same. As in the current integrator case (Figs. 8c and 8d), the distribution obtained in the alternative integrator for $h=10^{-6}$ shown in Fig. 8e exhibits an erroneous feature around $x_{1}=0$, while that for $h=10^{-3}$ (not shown) was relatively good.

Our remaining question is why we can obtain the distribution that is compatible with the exact BG distribution by numerical errors, although the system is nonergodic. We conjecture that this can be explained by the ergodic decomposition theorem [44]: (dynamics) a trajectory that is wandering in $A$ after starting at $\omega_{0}\in A$ is perturbed via noise to move out of $A$ and reaches $A^{[+]}$ or $A^{[-]}$, as illustrated in Fig. 8f; (measure) it then permits a description with a certain invariant measure (such as the Kolmogorov measure [5]) that would exist in (some stable subset of) $A^{[+]}$ or $A^{[-]}$. This should be related to the fact that the distribution in Case (C) was close to the exact one, as noted in the text. Furthermore, it may be applicable to other MD EOMs or other systems. For example, in the KM method (Example 3), we observed the erroneous behavior of $\gamma_{12}(\omega)$ from zero to nonzero values, as well as relatively good distributions for Cases (B) and (C), as stated in Appendix B. We believe that these conjectures will facilitate research on the measure theoretical study of a stable subset in an invariant space (a type of ergodicity in a local region), or a manner of handling numerical errors to lead to ergodicity (initiatively capturing the positive aspect of the numerical errors; stochastic methods can be viewed as such an approach). This should also contribute to revealing which initial conditions will result in error-induced ergodic-like behavior. The error-induced ergodic-like behavior is inclined to be compatible with stiff EOMs, as implied from the comparison between the NHC and KM methods, and may be related to Hörmander’s condition for certainly formulated noises [2, 29].