Dynamic principle for ensemble control tools
A. Samoletov B. Vasiev





A.Samoletov@liverpool.ac.uk

Department of Mathematical Sciences, University of Liverpool, Liverpool,
UK

Institute for Physics and Technology, Donetsk, Ukraine


B.Vasiev@liverpool.ac.uk


Department of Mathematical Sciences, University of Liverpool, Liverpool,
UK

 

Dynamical equations describing physical systems in contact with the
thermal bath are commonly extended by mathematical tools called “thermostats”. These tools are designed for sampling ensembles in statistical mechanics.
Here we propose a dynamic principle underlying a range of thermostats
which is derived using fundamental laws of statistical physics and
insures invariance of the canonical measure. 
The principle covers both stochastic and deterministic thermostat schemes.
 Our method has a clear advantage over a range of proposed and widely used thermostat schemes which are based on formal mathematical reasoning. Following the derivation of proposed principle we show its generality and illustrate its  
applications including design of temperature
control tools that differ from the Nosé-Hoover-Langevin scheme. 

Dynamic principle for ensemble control tools
    B. Vasiev
     
============================================





§ INTRODUCTION


Analysis of molecular systems is an essential part of research in
a range of disciplines in natural sciences and in engineering <cit.>.
As molecular systems affected by environmental thermodynamic conditions,
they are studied in the context of statistical physics ensembles.
Methods of dynamical sampling of the corresponding probability measures
are important for applications and they are under extensive study
and development <cit.>.
The traditional application of thermostats is molecular dynamics (MD),
that is sampling of equilibrium systems with known potential energy
functions, V(q), where q is a system's  configuration. However,
the ability to sample equilibrium ensembles at constant temperature
T would also imply the ability to sample arbitrary probability
measures. Indeed, as an alternative to the conventional MD practice,
one may use a probability density σ(q), theoretical or extracted
from experimental data, to define the potential function as V(q)=-k_BT lnσ(q),
where k_B is the Boltzmann constant.

Thermostats embedded into dynamical equations bring in the so-obtained
dynamics rich mathematical content. Such dynamical systems with an
invariant probability measure have become increasingly popular for
mathematical studies in a wide range of applications including investigation
of non-equilibrium phenomena <cit.>,
mathematical biology models <cit.>,
multiscale models <cit.>,
Bayesian statistics and Bayesian machine learning applications <cit.>,
superstatistics <cit.>.

Here, we present a unified approach for derivation of thermostats
sampling the canonical ensemble. The corresponding method is derived using fundamental
physical arguments that facilitate understanding physics of thermostat
schemes in general, and elucidate physics of the Nosé-Hoover (NH)
and the Nosé-Hoover-Langevin (NHL) dynamics in particular. Besides,
our method allows to build a plethora of thermostats, stochastic as
well as deterministic, including those previously proposed. We expect
that it can also be adjusted to arbitrary probability measures.

Classical mechanics and equilibrium statistical physics are adequately
described in terms of the Hamiltonian dynamics. Dynamic thermostat
schemes involve modified Hamiltonian equations of motion where certain
temperature control tools are included. The modified dynamics can
be deterministic as well as stochastic <cit.>.
Recently proposed NHL thermostats <cit.>
combine deterministic dynamics with stochastic perturbations. This
combination ensures ergodicity and allows “gentle” perturbation
of the physical dynamics that is often desired <cit.>.

To introduce our scheme, we consider a dynamical system S
consisting of N particles in d-dimensional space (𝒩=dN
degrees of freedom) described by the Hamiltonian function H(x),
where x=(p,q) is a point in the phase space ℳ=ℝ^2dN,
p={𝐩_i∈ℝ^d} _i=1^N are
momentum variables and q={𝐪_i∈ℝ^d} _i=1^N
are position variables. The Hamiltonian dynamics has the form, ẋ=J∇H(x)
in the phase space ℳ, where J
is the symplectic unit. The canonical ensemble describes the system
S in contact with the heat bath Σ (an energy reservoir
permanently staying in the thermal equilibrium with the thermodynamic
temperature T), and S may exchange energy with Σ
only in the form of heat. Thus, the temperature of the system S
is fixed while its energy, E, is allowed to fluctuate. The canonical
distribution has the form ρ_∞(x)∝exp[-β H(x)],
where β=(k_BT)^-1. On average along an ergodic trajectory
⟨ E(t)⟩ =E(T)=const. Rate of energy exchange
between the system S and the thermal bath Σ depends
on the temperature T. Note that Hamiltonian system is unable to
sample the canonical distribution since there is no energy exchange
between the system and the heat bath. To describe the heat transfer,
it is necessary to modify the equations of motion in a way that the
dynamics becomes non-Hamiltonian <cit.>. Suppose ẋ=G(x)
is a modified law of motion and Ḣ(x)=G(x)·∇H(x)
is the rate of energy change (depending on T) such that ⟨G(x)·∇H(x)⟩ =0,
that is the energy is constant on average. Let G(x)·∇H(x)∝ F(x,β)
where the temperature dependence is a key. In order to state the dynamic
principle governing temperature control tools, a few definitions are
required.



§ MICROSCOPIC TEMPERATURE EXPRESSIONS


Consider F(x,β) such that ⟨ F(x,β)⟩ =0
for all β>0. This condition is denoted as F(x,β)∼0
while the function F(x,β) is called the microscopic temperature
expression (TE).

For the system with H(x)=K(p)+V(q) examples of TEs include the
kinetic TE, F_kin(p,β)=2K(p)β-𝒩, and
the configurational TE, F_conf(q,β)=(∇V(q))^2β-Δ V(q)
<cit.>.

Various TEs can be obtained in the following manner. Suppose that
F(x,β) is a polynomial in β, F(x,β)=∑_n=0^2L+1φ_n(x)β^n∼0,
where L∈ℤ_≥0 and functions {φ_n(x)} _n=0^2L+1,
are subject to specification. Rewrite F(x,β) in the form

    F(x,β)=∑_k=0^L(φ_2k(x)+βφ_2k+1(x))β^2k∼0

for all β>0. Thus, from (<ref>) it follows that φ_2k(x)+βφ_2k+1(x)∼0
for all k∈{ 0,1,…,L}. 
To find φ_2k(x) and φ_2k+1(x) satisfying this condition consider the basic expression, F(x,β)=β φ_1(x)+φ_0(x).
Substituting φ(x)∂_iH(x) for φ_1(x),
where φ(x) is an arbitrary function, and utilizing the identity,
∂_ie^-β H(x)=-β∂_iH(x)e^-β H(x)
for all i=1,…,2dN and x∈ℳ, where ∂_i≡∂/∂ x_i,
we get ∂_iφ(x)+φ_0(x)∼0.
Then, excluding φ_0(x) from the basic expression, we arrive at F(x,β)=β φ(x)∂_iH(x)-∂_iφ(x)∼0 (or φ(x)∂_iH(x)-k_BT∂_iφ(x)∼0) for each and every x_i in ℳ provided that   φ(x)exp[-β H(x)]→0
as |x|→∞. This result can be represented in a compact form. Suppose φ_0(x)
is a vector field on ℳ such that φ_0(x)exp[-β H(x)]→0
as |x|→∞. Then 

    F_0(x,T)=φ_0(x)·∇H(x)-k_BT∇·φ_0(x)∼0.

This form of TE was previously discussed<cit.>. More
general TEs are allowed, e.g. vector fields F(x,β)=β ∇H(x)×φ(x)-∇×φ(x)∼0,
and so on. As a further generalization we introduce the notation

    F_l(x,T)=φ_l(x)·∇H(x)-k_BT∇·φ_l(x),

where l=0,1,...,L, φ_0(x)=φ(x),
and {φ_l(x)} _l=0^L
is a set of vector fields such that φ_l(x)exp[-β H(x)]→0
as |x|→∞. Then the general scalar
TE can be represented as

    F_L(x,T)=∑_l=0^LF_l(x,T)(k_BT)^2l∼0

for all L∈ℤ_≥0. A particular example of the use
of such a TE in a limited context (L=1 and φ_l(x)∝(p,0)
leading to the kinetic TE) can be found in the literature<cit.>.
In what follows we focus mainly on F_0(x,T) and only to a certain extent on F_L(x,T) where L≥1.

Although the expression (<ref>) implies the existence
of infinite number of TEs, they all are equivalent from the thermodynamic
perspective. However, the time interval required to achieve a specified
accuracy in ⟨ F(x,β)⟩ =0 can differ
for different TEs <cit.>. In general, physical
systems are often distinguished by multimodal distributions and by
existence of metastable states. Their dynamics is characterized by
processes occurring on a number of timescales. We assume that TEs
can be associated with dynamical processes occurring on various time
scales, and thus, they can be combined in multiscale models.



§ DYNAMIC PRINCIPLE


Now we claim the following dynamic principle for ensemble control
tools:  Let   F(x,T) be a TE. Then there exists the
dynamical system, ẋ=G(x), such that 

    ∇H(x)·G(x)∝ F(x,T).

Relationship (<ref>) states that the rates of dynamical
fluctuations in energy and in TE are proportional, both are zero on
average and there is no energy release along a whole trajectory in
ℳ. It is a necessary condition for any thermostat. In
what follows, with implication of the fundamental requirements of
statistical physics, we show that the relationship (<ref>)
leads to a general method for obtaining stochastic and deterministic
thermostats.

Let us consider the exchange of energy between the system S
and the thermal bath Σ. Any system placed in the heat bath
should in some extent perturb it and be affected by backward influence
of this perturbation. There exists a subsystem S_ad
of Σ such that S_ad is involved in a joint dynamics
with S. The rest of the heat bath is assumed to be unperturbed,
permanently staying in thermal equilibrium. This is an approximation
that is based on separation of relevant time scales. For instance,
Brownian dynamics assumes that characteristic time scales of S
and Σ are well separated and the system S does
not perturb Σ. If the time scale is refined (which is of particular
importance for small systems) then we have to take into account joint
dynamics of S and S_ad. We will show that
this case is closely related to NHL <cit.>
and NH dynamics <cit.>.

Thus, we have two cases: () the system S
doesn't perturb the thermal bath and there are no new dynamic variables.
The thermal bath in this case can only be taken into account implicitly
via stochastic perturbations (similar to the Langevin dynamics); ()
the system S perturbs a part (S_ad) of the
thermal bath Σ, while the rest of the thermal bath remains
unperturbed. We assume that there is no
direct energy exchange between S and Σ. 
Fundamentals of the statistical mechanics require that
the systems S and S_ad are statistically
independent at thermal equilibrium. 
Let us consider cases  and 
in detail.



 §.§ Stochastic dynamics


Suppose  ∇H(x)·ẋ=λ F_0(x,T),
where λ is a constant. Without loss of generality, we can
consider modified Hamiltonian dynamics in the form, ẋ=J∇H(x)+ψ(x,λ),
and consequently: 

    ∇H(x)·ψ(x,λ)=λ F_0(x,T),

where the vector field ψ(x,λ)
is to be found. Since the thermal bath does not appear in equation
(<ref>) explicitly, only stochastic thermal noise may
be involved in the dynamics. To find ψ,
we introduce 2𝒩-vector of independent thermal white noises,
ξ(t), such that ⟨ξ(t)⟩ =0,
⟨ξ_i(t)ξ_j(t')⟩ =2λ k_BTδ_ijδ(t-t'),
and the vector field, Φ(x), such that

    ⟨ξ(t)·Φ(x)⟩ =λ k_BT ⟨∇·φ(x)⟩,

where ⟨⋯⟩ is the Gaussian average
over all realizations of ξ(t). Using Novikov's formula
<cit.>, we get

    ⟨ξ(t)·Φ(x)⟩ = ∑_i,k⟨∂Φ_k/∂ x_iδ x_i(t)/δξ_k(t)⟩λ k_BT.

Suppose δ x_i(t)/δξ_k(t)=ζ_i(x)δ_ik,
where the vector field ζ(x) is such that each
component ζ_i(x) does not depend on x_i, that is

    ∇∘ζ(x)=0,

where ∘ denotes the component-wise (Hadamard) product of two
vectors and 0 is the null vector. Then ∇·φ(x)=∇·(ζ(x)∘Φ(x)).
Thus, we get φ(x)=ζ(x)∘Φ(x)
and it follows that Φ(x)=ζ^-1(x)∘φ(x),
where ζ^-1(x) is the vector field such that
ζ^-1(x)∘ζ(x)=1.
Assuming φ(x)=η(x)∘∇H(x),
where η(x)≡ζ(x)∘ζ(x),
we get 

    ψ(x,λ)=-λη(x)∘∇H(x)+ζ(x)∘ξ(t)

and the modified Hamiltonian dynamics takes the form of stochastic
differential equation (SDE): 

    ẋ=J∇H(x)-λη(x)∘∇H(x)+ζ(x)∘ξ(t).


The Fokker-Planck equation (FPE) corresponding to SDE (<ref>)
has the form ∂_tρ=ℱ^*ρ,
where 

    ℱ^*ρ=-J∇H(x)·∇ρ+λ∇·[η(x)∘∇H(x) ρ]
    
    +λ k_BT∇·[η(x)∘∇ρ].

Note that the last term here was found using the following specific relationship for the vector field
ζ(x):

    (ζ(x)∘∇)·(ζ(x)∘∇ρ)=∇·[η(x)∘∇ρ].

 
Invariant probability density for dynamics (<ref>)
is determined by the equation ℱ^*ρ=0.
It is expected that this is a unique invariant density <cit.>. 

We claim that for the defined above vector field ζ(y)
the canonical density, ρ_∞∝exp[-β H(x)],
is invariant for the stochastic dynamics given by (<ref>),
that is ℱ^*ρ_∞=0.  The proof
is by direct calculation.

The Langevin equation is a particular case of (<ref>).
For example, for the system with H(x)=p^2/2m+V(q),
where x=(p,q)∈ℝ^2 we have:


    if ζ=(1,0), then
 ṗ=-V'(q)-λp/m+ξ(t),
	q̇=p/m;

  if ζ=(0,1), then
 ṗ=-V'(q), q̇=p/m-λ V'(q)+ξ(t).


The procedure for obtaining stochastic dynamics (<ref>) is essentially a general and can be a quite straightforwardly extended to other TEs, for example, the general scalar TE (<ref>).
Indeed,
let us introduce the set of 2𝒩-vectors of independent
thermal white noises, {ξ(l;t)} _l=0^L, L∈ℤ_≥0,
such that ⟨ξ(l;t)⟩ =0,
⟨ξ_i(l;t)ξ_j(l';t')⟩ =2λ_lk_BTδ_ijδ_ll'δ(t-t'),
and the set of vector fields, {ζ(l;x)} _l=0^L, L∈ℤ_≥0,
such that ∇∘ζ(l;x)=0
for any l≥0, where ∘ denotes the component-wise (Hadamard)
product of two vectors and 0 is the null vector. Starting
from the relationship,

    ∇H(x)·ψ(x,λ)=∑_l=0^Lλ_lF_l(x,T) (k_BT)^2l,

and then strictly following arguments as stated above, we get

    ψ(x,λ)=-∑_l=0^Lλ_lη(l;x)∘∇H(x) (k_BT)^2l
    
    +∑_l=0^Lζ(l;x)∘ξ(l;t) (k_BT)^l.

 where η(x)≡ζ(x)∘ζ(x).
Thus, we arrive at the following stochastic dynamics

    ẋ=J∇H(x)-∑_l=0^Lλ_lη(l;x)∘∇H(x) (k_BT)^2l
    
    +∑_l=0^Lζ(l;x)∘ξ(l;t) (k_BT)^l.

One can verify that the canonical measure is invariant for this stochastic
equation of motion. Generally speaking, the dynamics 
(<ref>)
includes 2𝒩(L + 1) independent white noise processes. This seems impractical. However,
we can point out that (<ref>) potentially
useful for multi-timescale stochastic simulations. As a simple example,
let H(x)=p^2/2m+V(q), L=1, ζ(0;x)=(1,0),
and ζ(1;x)=(0,1), then we arrive at the stochastic
dynamics with two timescales involved,

    ṗ=-V'(q)-λ_0p/m+ξ_p(0;t),
    q̇=p/m-λ_1(k_BT)^2 V'(q)+k_BT ξ_q(1;t),

where ⟨ξ_p(0;t)⟩ =0, ⟨ξ_q(1;t)⟩ =0,
⟨ξ_p(0;t)ξ_q(1;t)⟩ =0, ⟨ξ_p(0;t)ξ_p(0;t')⟩ =2λ_0k_BTδ(t-t'),
⟨ξ_q(1;t)ξ_q(1;t')⟩ =2λ_1k_BTδ(t-t'),
as specified above. 
Analysis of p- and q-dynamics can be performed in reduced systems following the separation of these variables according to their time scales
<cit.>.



 §.§ Deterministic and stochastic dynamics


Let S_ad be associated with an even-dimensional
phase space ℳ_ad, the Hamiltonian function
h(y), y∈ℳ_ad, and the Hamiltonian dynamics,
ẏ=J_y∇_yh(y), where J_y
is the symplectic unit. Without loss of generality, we can assume
that the modified Hamiltonian dynamics of the system composed by S
and S_ad has the form, 

    ẋ=J_x∇_xH(x)+ψ(x,y),
    ẏ=J_y∇_yh(y)+ψ^∗(y,x),

where ψ(x,y) and ψ^*(y,x)
are vector fields on ℳ and ℳ_ad
correspondingly. To derive deterministic dynamics, let us temporarily
ignore the heat exchange between S_ad and Σ, that is, ∇_yh(y)·ψ^∗(y,x)=λ^∗F^∗(y,T)
and ∇_xH(x)·ψ(x,y)=λF(x,T).
As discussed above, these  relationships lead to the stochastic dynamics.
Systems
S and S_ad must be statistically independent
in the thermal equilibrium, so that ∇_xH(x)·ẋ∼0
and ∇_yh(y)·ẏ∼0 are satisfied
simultaneously. Thus, we assume that 

    ∇_xH(x)·ψ(x,y)     =g(x)F_0^∗(y,T),
    ∇_yh(y)·ψ^∗(y,x)     =-g^∗(y)F_0(x,T),

where g(x) and g^∗(y) are some vague functions, and

    F_0(x,T)=φ(x)·∇_xH(x)-k_BT ∇_x·φ(x),
    
    F_0^∗(y,T)=Q(y)·∇_yh(y)-k_BT ∇_y·Q(y),

are TEs for the systems S and S_ad correspondingly.
These relationships are valid for any H(x) and h(y). To specify
ψ(x,y) and ψ^*(y,x), we assume
that g(x)=a(x)·∇_xH(x),
g^∗(y)=b(y)·∇_yh(y),
where a(x) and b(y) are vector fields
on ℳ and ℳ_ad, respectively.
It follows that

    ψ(x,y)=a(x) F_0^∗(y,T),  ψ^∗(y,x)=b(y) F_0(x,T).

To determine the relationship between the vector fields a(x),
b(y) and TEs F_0(x,T), F_0^∗(y,T),
recall that if the combined system S+S_ad
is isolated, then Ḣ(x)=-ḣ(y); and if T≠0, then Ḣ(x)+ḣ(y)∼0. 
Straightforward calculations show that

    a(x)=φ(x),  b(y)=Q(y),

provided that b(y)exp[-β h(y)]→0
as |y|→∞ and a(x)exp[-β H(x)]→0
as |x|→∞. As a result, we have the equations of motion 

    ẋ=J_x∇_xH(x)+F_0^∗(y,T)φ(x),
    ẏ=J_y∇_yh(y)-F_0(x,T)Q(y),

which are generalized NH equations. 

The Liouville equation associated
with the system (<ref>) has the form ∂_tρ=-ℒ^*ρ,
where ℒ^*ρ=∇_x·(ẋρ)+∇_y·(ẏρ).
Invariant probability densities are determined
by the equation ℒ^*ρ=0. 
We claim that if Q(y) and φ(x)
are the defined above vector fields, then the canonical density ρ_∞∝exp[-β H(x)]·exp[-β h(y)]
is invariant for dynamics (<ref>), that is ℒ^*ρ_∞=0.
The proof is by direct calculation.

As a particular case, let Q(y) be an incompressible
vector field (i.e. ∇_y·Q(y)=0
for all y∈ℳ_ad). Then we arrive at the NH
equations 

    ẋ=J_x∇_xH(x)+(Q(y)·∇_yh(y))φ(x),
    ẏ=J_y∇_yh(y)-F_0(x,T)Q(y).


Now we include into our consideration the effect of the
thermal bath Σ on S_ad dynamics, that is the
relationship ∇_yh(y)·ψ^∗=λ F_0^∗(y,T).
Following the arguments and notations used to derive SDE (<ref>),
we arrive at the stochastic dynamics: 

    ẋ=J_x∇_xH(x)+F_0^∗(y,T) φ(x),
    ẏ=J_y∇_yh(y)-F_0(x,T) Q(y)-λη(y)∘∇_yh(y)
    
    +ζ(y)∘ξ(t),

which are generalized NHL equations<cit.>.
In the particular case of an incompressible vector field Q(y)
we get the NHL equations: 

    ẋ=J_x∇_xH(x)+(Q(y)·∇_yh(y)) φ(x),
    ẏ=J_y∇_yh(y)-F_0(x,T) Q(y)-λη(y)∘∇_yh(y)
    
    +ζ(y)∘ξ(t),


FPE corresponding to (<ref>) has the form ∂_tρ=ℱ^*ρ,
where 

    ℱ^*ρ=-J_x∇_xH(x)·∇_xρ-J_y∇_yh(y)·∇_yρ
    
    -F_0^∗(y,T)∇_x·[φ(x)ρ]+F_0(x,T)∇_y·[Q(y)ρ]
    
    +λ k_BT∇_y·[η(y)∘∇_yρ]+λ∇_y·[η(y)∘∇_yh(y)ρ].

Invariant probability density for the SDE (<ref>) is determined
by the equation ℱ^*ρ=0. 

We claim that if Q(y), φ(x),
and ζ(y) are the defined above vector fields ,
then the canonical density, ρ_∞∝exp[-β H(x)]·exp[-β h(y)],
is invariant for the NHL dynamics (<ref>), that is ℱ^*ρ_∞=0.
The proof is by direct calculation. 

Besides, we expect that this dynamics is ergodic <cit.>.

Commonly used NH <cit.> and NHL <cit.> thermostats are particular cases of thermostats given by (<ref>)
and (<ref>) correspondingly. For example, by substituting
ζ^2/2Q for h(y), y=(ζ,η)∈ℝ^2,
(-Q,0) for Q(y) and (𝐩,0) for φ(x)
in (<ref>) we get classical NH equations <cit.>.

It is worth to note that the case of the general TE can be considered
straightforwardly following the method of dynamic principle, as developed
above. Assume that

    ∇_xH(x)·ψ(x,y)     =∑_l=0^Lg_l(x)F_l^∗(y,T) (k_BT)^2l,
    ∇_yh(y)·ψ^∗(y,x)     =-∑_l=0^Lg_l^∗(y)F_l(x,T) (k_BT)^2l.

These relationships must be valid for any H(x) and h(y). To specify
ψ(x,y) and ψ^*(y,x), we set
g_l(x)=a_l(x)·∇_xH(x),
g_l^∗(y)=b_l(y)·∇_yh(y),
from what follows that
a_l(x)=φ_l(x)
and b_l(y)=Q_l(y). Thus,

    ψ(x,y)=∑_l=0^LF_l^∗(y,T) (k_BT)^2lφ_l(x),
    ψ^∗(y,x)=-∑_l=0^LF_l(x,T) (k_BT)^2lQ_l(y).

Finally, we arrive at the deterministic equations of motion (modified Hamiltonian dynamics),

    ẋ=     J_x∇_xH(x)+∑_l=0^LF_l^∗(y,T) (k_BT)^2lφ_l(x),
    ẏ=    J_y∇_yh(y)-∑_l=0^LF_l(x,T) (k_BT)^2lQ_l(y).

We will not discuss the equations (<ref>) in detail and only note that the canonical measure is invariant for this dynamics, and a generalization to stochastic NHL type dynamics can be obtained. Strictly speaking, such a generalization is important since it simulates an equilibrium reservoir of the energy and ensures the ergodicity of dynamics. To outline a connection
between equations of motion (<ref>) and known deterministic thermostats<cit.>,
we provide the following simple example. Let L=1, H(x)=p^2/2m+V(q),
h(y)=η_0^2/2Q_0+η_1^2/2Q_1,
φ_0(x)=(p,0), φ_1(x)=(p^3,0),
Q_0(y)=(-Q_0,0,0,0), and Q_1(y)=(0,-Q_1,0,0),
then

    ṗ=-V'(q)-η_0p-η_1k_BT p^3,
    q̇=p/m,
    η̇_0=Q_0(p^2/m-k_BT),
    η̇_1=Q_1(p^4/m-3k_BTp^2)(k_BT)^2,

the  dynamic equations equipped with the control of first two moments of the equilibrium kinetic energy<cit.>. Similarly, we can obtain dynamic equations that control the configurational temperature moments. 



§ REDESIGN OF NHL THERMOSTAT


In this Section we consider an alternative to the conventional NH and NHL thermostat schemes. 
This alternative (seen as a particular case of dynamical equations (<ref>)) 
is based on the consideration of physically reasonable chain of interactions, 
S↭ S_ad↭Σ, that is, the system S_ad is a buffer between the physical system S and the infinite energy reservoir Σ. 

Consider the dynamical equations (<ref>) and (<ref>),
and assume that ∇_x·φ(x)=0, ∇_y·Q(y)≠0.
Note, that these assumptions are opposite to the requirements for the NH and NHL dynamics, where ∇_x·φ(x)≠0, ∇_y·Q(y)=0.
We get


ẋ     =J_x∇_xH(x)

     +[Q(y)·∇_yh(y)-k_BT ∇_y·Q(y)] φ(x),

ẏ     =J_y∇_yh(y)-(φ(x)·∇_xH(x))Q(y),

and

ẋ     =J_x∇_xH(x)

     +[Q(y)·∇_yh(y)-k_BT ∇_y·Q(y)] φ(x),

ẏ     =J_y∇_yh(y)-(φ(x)·∇_xH(x)) Q(y)

     -λη(y)∘∇_yh(y)+ζ(y)∘ξ(t),

where vector fields involved are such as indicated above. Thus, there is plenty of freedom in specification of particular thermostat equations of motion. 

To illustrate the redesigned NH and NHL thermostat dynamical systems (described by the equations (<ref>) and (<ref>) correspondingly)
let us consider system S with the Hamiltonian function
H(p,q),

    H(p,q)=p^2/2m+1/2mω^2q^2,   x=(p,q)∈ℝ×ℝ,

that is a harmonic oscillator of mass m and frequency ω,
and system S_ad with the Hamiltonian function h(v,u),

    h(v,u)=v^2/2μ,   y=(v,u)∈ℝ×ℝ,

that is a free particle of mass μ. Harmonic oscillators are among
central instruments in analysis of many physical problems, classical
as well as quantum mechanical. It is known that generating the canonical
statistics for a harmonic oscillator is a hard problem. For example,
the NH scheme is proven to be non-ergodic<cit.> and the NHL scheme<cit.>, and earlier the NHC scheme<cit.>,
was proposed to overcome this difficulty. Anyway, it is important for any dynamic thermostat
to correctly generate the canonical statistics for a harmonic oscillator.

The deterministic thermostat dynamics (<ref>) as
well as stochastic dynamics (<ref>) allow a plethora
of further specifications. To be as close as possible to redesign
of original NH dynamics<cit.>, we set Q(y)=(v,0),
∇·Q=1, and φ(x)=(γ,0),
where γ is a dimensional parameter, ∇·φ=0.
Thus, we arrive at the following equations of motion: 

ṗ     =-mω^2q+γ[v^2/μ-k_BT],

q̇     =1/mp,

v̇     =-γp/mv,

u̇     =v/μ;

and

ṗ     =-mω^2q+γ[v^2/μ-k_BT],

q̇     =1/mp,

v̇     =-γp/mv-λv/μ+ξ(t),

u̇     =v/μ;

where ζ=(1,0) and ⟨ξ(t)ξ(t')⟩ =2λ k_BTδ(t-t'). Note, that  equations (<ref>) and (<ref>) are redesign
of NH (denote RNH) and NHL (RNHL) thermostats correspondingly.

System (<ref>) has two integrals of motion, that is

    I_1=vexp(γ q)=const,
    
    I_2=p^2/2m+1/2mω^2q^2+v^2/2μ+γ k_BT/mq=const,

indicating the lack of ergodicity. For example, if all parameters of the system
(<ref>) are set equal to unity, m=1,ω=1,μ=1,γ=1,k_BT=1,
and initial conditions are p=1,q=0,v=1, then the phase trajectory is represented by the closed curve
and the Poincaré section (p,q) shown on Figure <ref>. This is expected from the existence
of two integrals of motion, that is I_1 and I_2. It is clear
that the trajectory does not explore the phase space available for
the harmonic oscillator. This ergodicity problem is not surprising,
the convenient NH dynamic suffer from the same problem. It is questionable  that the situation can be improved with a more complex φ and Q,
for example,
φ=(γ_1,γ_2), φ=(γ_1mω^2q,γ_21mp),
Q=(v,u), and so on. If φ=(γ_1,γ_2),
then we get

ṗ     =-mω^2q+γ_1[v^2/μ-k_BT],

q̇     =p/m+γ_2[v^2/μ-k_BT],

v̇     =-(γ_1p/m+γ_2μω^2q)v,

u̇     =v/μ,

and it is easy to show that this dynamics is not ergodic.



Our next illustration will be devoted to the system described by thermostat dynamical equations (<ref>). We will show, by means of numerical simulations, that a certain realization of the whole length chain of physically
reasonable interactions, that is S↭ S_ad↭Σ,
generates the correct statistics.

Let us consider the case when all parameters of the system (<ref>) are set equal to unity, m=1,ω=1,μ=1,γ=1,k_BT=1,λ=1,
and the initial conditions are: p=0,q=0,v=0. Phase trajectories of length 10^6 are generated using the Euler method with a time step of Δ t=0.0005. We have repeated simulations using the fourth-order Runge-Kutta method with a random contribution held once for the entire interval from t to t+Δ t, and arrive at the same result.

Figure <ref>(a) shows the Poincaré section (p,q) for a harmonic oscillator equipped with the temperature control tool (<ref>). This figure demonstrates that the trajectory generates proper sampling of the full phase space of the harmonic oscillator. Figures <ref>(b) and <ref>(c) show the momentum and position distribution functions from simulations as 
compared with the exact analytical expressions. In both cases, the
Gaussian distribution is generated in agreement with the theoretical prediction. Presented results serve as an evidence of ergodic sampling the canonical statistics. 



A key difference between the NHL and RNHL schemes is that the latter relates
the temperature control tool to the system S_ad rather than to the system S, and the corresponding variable, v, must be Gaussian, according to the equations (<ref>). Thus, it is important that the RNHL dynamical equations properly generate the Gaussian
statistics of v variable. Figure <ref> shows the v-distribution function from simulations as compared with the exact analytical solution and indicates a good agreement between them. 





§ CONCLUSION


In conclusion, we emphasize that the method proposed in this work
is based on the fundamental laws of statistical physics and offers
a unified approach in developing stochastic and deterministic thermostats. For clarity of presentation we have illustrated our method using a few simple TEs and restricted our consideration by Markov dynamics. The presented method allowed us to obtain a wide spectrum of stochastic and deterministic dynamical systems with the invariant canonical measure. We note that the idea
of presented method is general and adaptable to a variety of TEs so
that it can be used to produce thermostats of novel types. For example the thermostat for the system with non-Markov dynamics, i.e. the one described by the equation ∇H(x(t))·ẋ(t)∝∫_0^tdt'G(t-t')F(x(t'),T). As a second example of new type of thermostats we can mention the one for the gradient dynamical system. 

We realize that non-trivial new thermostats should be verified by
test simulations. In our follow up work we will focus on these and other applications of the presented method. 

This work has been supported by the BBSRC grant BB/K002430/1 to BV. 


apsrev4-1


