Double density dynamics: realizing a joint distribution of a physical system and a parameter system

To perform a variety of types of molecular dynamics simulations, we created a deterministic method termed ‘double density dynamics’ (DDD), which realizes an arbitrary distribution for both physical variables and their associated parameters simultaneously. Specifically, we constructed an ordinary differential equation that has an invariant density relating to a joint distribution of the physical system and the parameter system. A generalized density function leads to a physical system that develops under nonequilibrium environment-describing superstatistics. The joint distribution density of the physical system and the parameter system appears as the Radon–Nikodym derivative of a distribution that is created by a scaled long-time average, generated from the flow of the differential equation under an ergodic assumption. The general mathematical framework is fully discussed to address the theoretical possibility of our method, and a numerical example representing a 1D harmonic oscillator is provided to validate the method being applied to the temperature parameters.


Introduction
Numerical simulations of physical systems are important tools for understanding the structures, phases, and mechanisms of the systems in a systematic manner. Molecular dynamics (MD) simulations are used for classical/quantum Hamiltonian/non-Hamiltonian systems to help understand the characteristics of the systems in terms of their microscopic descriptions [1].
The physical systems, or rather the equations of motion in these methods, include some parameters that may be key for performing the simulations. For example, intensive-quantity parameters such as temperature and pressure are important for understanding equilibrium and nonequilibrium states. There will be cases in which we are motivated to fluctuate these values to ascertain the effects of these parameters. A typical situation is to change or fluctuate these intensive quantity parameters, which is relevant for nonequilibrium simulations [2,3]. It may also be useful to fluctuate the temperature to obtain a good sampling method [4,5]. In other cases, there might be a parameter for which the best way to handle it is unclear. For example, the optimal values of the parameters are unknown a priori, such as the 'mass' parameters for the thermostat or barostat in extended system methods [6,7]. In addition, the optimal parameter values in the sampling method are often system-dependent, and efforts are required for their identification. For example, the optimal value of the extensive parameter q in the Tsallisdistribution sampling [8,9] is unknown a priori.
In these situations, what should we do? Ad hoc methods for varying the parameter values do not provide important information, such as the probability distribution that the physical system obeys, and they thus preclude an easy interpretation to compare the results with experiments. The underlying motivation of the current study is to construct a route to solve this problem. A possible approach, which is the topic of this paper, is to control the dynamical variations of the parameter values by a certain protocol, based on a probabilistic description. Specifically, along the contexts introduced above, we consider performing the following three types of MD simulations. First, (I) temperature is fluctuated, obeying a certain probability distribution (for a nonequilibrium simulation or for good sampling). Second, taking account of a more general situation, where the target is not only the temperature parameter but also an arbitrary parameter, we consider that (II) any parameter is fluctuated, obeying a certain probability distribution (for varying an unknown parameter value). Third, considering the fact that such a fluctuating parameter yields a certain system other than the originally given physical system, a more general situation is that (III) any other system is fluctuated, obeying a certain probability distribution (for coupling two systems, e.g. all-atom model and the associated course-grained model [10]). Namely, our purpose is to produce a general deterministic equation of motion (EOM) in order to perform many kinds of MD simulations represented by these three categories (I)-(III), which are sorted in ascending order of generality.
How can we realize these methods?
The key is to focus on rigorously connecting the probability densities and the dynamics that generate the densities. In this regard, we are fortunate to have a rich background to draw from. To create the Boltzmann-Gibbs (BG) (or canonical) distribution density described with the coordinate variables x x x , , n n 1 via a deterministic EOM, the Gaussian isokinetic dynamics [11,12] can be used. Here U x ( ) is the potential energy of the physical system, and β is the inverse temperature (Boltzmannʼs constant is unity). Similarly, to create the BG distribution density described with both x and momentum variables p p p , , n n 1 via a deterministic EOM, the Nosé-Hoover (NH) dynamics [13,14] can be used. This is represented by where ζ is a variable to control the kinetic energy, K p ( ), by mimicking a function of a heat bath with a temperature of 1 b, and γ is an associated parameter. The NH scheme has been actively utilized in revisiting thermostatting methods as detailed in a number of reviews [15,16].
Based on this background, we will create an arbitrary distribution density of the form, Phys ( ) ( ) ( ) ( ) r r ¼ µ via a deterministic EOM, which is sought for. Here x p Z , ; Phys ( ) r stands for a density value of a conditional probability for the physical system, such that the physical system variable x p , ( ) emerges under the condition that 'parameter' Z takes a certain value. f Z ( ) stands for a density value of the probability for the parameter itself. The left-hand side x p Z , , , ( ) r ¼ thus represents a joint distribution density defined in 'physical-space × parameter-space. ' In the context of the three issues (I)-(III), this concept is specifically demonstrated, as follows. The situation that (I) 'temperature 1 b is fluctuated obeying a certain probability distribution f ' is described as a density, x p x p f , , , , ; . 5 x p , ; Phys ( ) r b corresponds to the rhs of equation (2). However, β in equation (5) is not a constant parameter but a dynamical variable obeying statistics characterized by the density, f. Namely, we can imagine a physical view such that the target physical system interacts with a heat bath with a temperature that is instantaneously 1 b, and this temperature itself fluctuates throughout the duration, according to the probability distribution density f. The situation that (II) 'any parameter q is fluctuated obeying a certain probability distribution f' is described as a density such as x p q x p q f q , , , , ; . 6 The situation that (III) 'any other system described by X is fluctuated obeying a certain probability distribution f' is described as a density such as x p X x p X f X , , , , ; . 7 In this case, X may no longer be a 'parameter,' but it can be viewed as a generalized parameter concerned with the original physical system. Below, for simplicity, instead of using Z in the most general expression (4), we just use β as in the case of (I), and express the density as This is also because we shall mainly emphasize the temperature parameter in this paper, since temperature fluctuations are one of the most interesting and important issues in a variety of research areas. Namely, we realize a joint distribution density of the physical system and the parameter system, represented as the form (8). We derive an ordinary differential equation (ODE) that produces the density, equation (8), with dynamical variables including x p , , and β, by providing an invariant density (a density of an invariant measure) that is related to equation (8).
In terms of the superstatistics, our method focusing on situation (I) leads to a deterministic, time-reversible method which generates a superstatistics distribution in a dynamical manner. Generating density (8), we observe that the superstatistics distribution is realized by focusing only on the physical variables x p , ( ) in our method. Our approach is the dynamical realization of superstatistics distributions, namely, the realization of the integrand in equation (9), and is not limited to the realization of the integral in equation (9) ('dynamical realization' thus does not only mean realizing the superstatistics in a certain dynamics, but also realizing the density (8), which generates x p , R ( ) r , by dynamically fluctuating β). With this realization, a physical system in a nonequilibrium environment, which is yet under control in light of the probability distribution (8), can be simulated. It must provide new gains, due to the freedom of choosing Phys r and f. Approaches to realize the marginal distribution, x p , R ( ) r , have been taken into account [32]. In other words, they utilize the distribution in static manners, for which β is not a dynamical variable, but just a variable of integration. One of the advantages of the dynamical realization of the superstatistics distribution, as compared with the static realization of the superstatistics distribution, is that we can directly observe the influence of temperature deviations on the physical system, leading to an understanding of the dynamical features of the physical system in nonequilibrium. Another advantage is that we can generate the distribution R r , as well as reweighted distributions, even if the integration in equation (9), which defines x p , R ( ) r , cannot be analytically accomplished (the static method often needs the nontrivial analytical form of R r ). The resulting ability to freely set f guides us to direct interpretations between the simulations and the experiments. In addition, our target is a general density form of equation (8), in contrast to conventional approaches addressing a specific form, On another front, the generation of x p , R ( ) r based on a stochastic EOM, the Langevin equation, was discussed by Beck [17], for which f ( ) b is assumed to exist irrespective of descriptions of the detailed structure. Hahn et al [33] characterized the EOM as stochastic differential equations driven by exchangeable processes. In our method, however, β is a realistic dynamical variable, and the specific mechanism to realize f ( ) b is provided. Moreover, the current deterministic method enables us to monitor the numerical accuracies in integrating the ODE [34].
Note that we shall view the variable β, employed in equation (9), as (not only the inverse temperature but also) a general parameter corresponding to issues (I)-(III) as stated above. In this sense, the term, 'superstatistics distribution' associated with β, will be used. It should also be noted that our aim in the present work is not to dynamically justify superstatistics, but to generate distributions compatible with the superstatistics description as defined above. Such a justification and the relevant discussions related to it are provided in other work, such as the literature mentioned above.
In this paper, we pursue a method for generating a distribution density, equation (8), in a deterministic manner by constructing an ODE. We establish a general framework of the method and consider a mathematical structure of the resulting space of dynamical variables. In section 2 we present a vector field defining the ODE, as shown in equation (15). Then we clarify the structure of the generated probability space, to demonstrate what the realization of the probability density accurately means. In section 3 we define the total phase-space densityfunction, ρ, suitable for our purpose, which includes the term corresponding to (8) (see equation (33)). While ρ contains other required functions that should be specifically defined, its suitable definition will lead us to the attainment of our purpose. This is ensured by theorem 6, which states the fundamental results regarding the realization of equation (8). The realization of R r is established in corollary 7. In section 4 we perform the remaining task, that is, the setting of the required functions in ρ, and discuss the related mathematical conditions. Note that while we can actually perform individual simulations by this procedure, we should focus on situations (I) and (II) due to the specific setting. As discussed in section 5, the numerical simulation of a simple model system demonstrates the validity of the current method, in terms of the generated distribution and the kinetics of the physical and temperature variables. We conclude our study in section 6.

Double density dynamics
Main idea-The target probability density, ρ, of variables ω, including x p , and Q, where Q is a dynamical variable related to β, can be realized by a smooth vector field X that has an . Namely, if we ensure a normal condition for ρ , then for any solution f of an ODE, , generated by a complete field X, and for any P- under an ergodic assumption [35]. Although it is not trivial to construct such a suitable X for a given ρ, the density dynamics (DD) approach [8,36], modeled on the NH method [13,14], gives a solution for this, a kind of inverse problem. Generalizing this idea, we propose a vector field that reflects the two systems: the physical system and the parameter system. The DD for a physical system and the DD for a parameter system will be coupled in a consistent manner.
Specific strategies-To obtain X and the associated probability space we proceed as follows: (1) Construct a vector field X such that X div 0 r = holds, provided that ρ is given; (2) Establish a probability space, which is needed to accurately demonstrate the realization of the density. First, define an ODE via the obtained X along with some technical contexts. Second, assume the ergodicity. Third, calculate the long-time averages of the phase-space functions. Since we will utilize a scaled time average at this point, we should manifest the resulting probability-space structure compatible with this scaling. Accordingly we can say that ρ is realized.

Strategy (1): Construct a vector field
We begin with the definitions of the phase space Ω and then make a general form of the vector field X defined on Ω.
Let the phase space Ω be a domain (an open, connected set) of N  with N n n 2 1 with (·|·) being the inner product in n a ( )  , for a 1, 2 = , and ln . 13 ( ) r Q º -Here h a ( ) and k a ( ) (a 1,2 = ) are arbitrary C 2 -functions on Ω, and they will be suitably chosen according to the problem as seen below. Namely, X is constructed from the given ρ and subsidiary functions, h a ( ) and k a ( ) . If we consider only a = 1 (that is, we only consider the physical system) and put h k 1 . We can confirm the smoothness (C 1 ) of X easily and the validity of the Liouville equation via a straightforward calculation. The flow T t { } generated by X becomes the basis of our method. We call this flow double density dynamics (DDD).

Strategy (2): Establish a probability space
To establish a probability space defined on the 'physical-system-space × parameter-systemspace' and examine its structure, we will emphasize the mathematical issues involved. One can, however, consider a temperature system instead of a general parameter system to easily obtain a physical image in oneʼs mind, although the issues in this subsection are applicable to a general parameter system appearing in situations (I)-(III) as demonstrated in section 1.
To proceed with the specifics, first we rewrite x p x p , , , are the physical variables with n n 1 ( ) º degrees of freedom, and 1 1 ( )   z z º is a control variable for the physicalsystem temperature (here just imagine the NH mechanism for x p , ( ) via a friction variable ζ). On the other hand, Q P x p , , , forms the parameter system yielding a nonequilibrium environment against the physical system, via coordinates Q n m , and a control variable 1  h Î , which again forms an analog of the NH mechanism for Q, ( )  , as seen below. Note that, as suggested, β will be given as a function of Q, that is, the target general 'parameter' β (which is now actually not just a constant parameter, but a dynamical variable) is treated by handling Q. This is due to a technical reason, which will be discussed in section 3.1 in detail. In any case, we treat β in a general manner.
Next, to make an intended ODE, two technical issues are presented. Firstly, we should specifically define the subsidiary functions, and thus we use for a certain strictly-positive, C 2 function T. Although equation (14) is not the unique choice of the function forms, this form may be a natural choice considering the fact that the fundamental variables we consider are x p Q , , ( ) and that equation (14a) is utilized in the original DD, as stated above. Secondly, we multiply X by the scalar field T. In general, scalarfield multiplication, that is, X X k  using a smooth, strictly-positive scalar field :  k W  , is a general concept to re-parameterize the time for each orbit of the flow. Namely, it causes a time rescaling between a solution of the ODE is to obtain a simple form of an EOM, as described below. This is also related to the fact that the DD also use a scalar-field multiplication with the form of T ( ) k w º ¢ = const. Although the resulting rescaling via such a constant field is trivial, the specific effect of the current rescaling on the time averages of the phase-space functions in the DDD should be clarified, and will be investigated below.
Hence we have X : , 15 We therefore get the ODE, which defines the time development of the total system. This is read, e.g. for the x component, as indicating that the derivative of each coordinate with respect to time is given by )has a dimension of energy or temperature (e.g. see equation (18); note that Boltzmannʼs constant is unity). Note that a simple physical view of function T derived from equations (15) and (17) is that the two systems-one the physically based system represented by x p , , ( ) z and the other the parameter system represented by Q, , ), which depends only on the fundamental variables x p Q , , ( )and plays the role of energy or temperature as implied by its physical dimension. The next task is to show that the obtained ODE (17) produces an individual orbit that has a required property: namely, by assuming ergodicity, a time average equals the space average for any phase-space function g. Here, the time average is not the simple one (which is denoted by g ), but a scaled one, due to the scalar-field multiplication adopted above. We shall show g g , 1 9 ( ) = where the lhs of equation (19) is a scaled long-time average of g, defined by g gT T , 2 0 ( ) º and the rhs of equation (19) is the space average of g weighted by ρ and written as We prove equation (19) by proposition 1 and then clarify the resulting probability space compatible with the scaled time average. Before doing so, to accurately state these and subsequent issues we shall use the following definitions.
Definitions and notations. We define projections: : To clarify a variable dependence, we denote T T : , . For any n  Î , l n and n  represent the Lebesgue measure and the Lebesgue measurable sets on n  , respectively; n  represents the Borel measurable sets on n  ; n  Ç W is denoted by n  W and so on. Basically, we consider a measure space P , , using a modified density r ¢ instead of the given density ρ, as described below. Put L P g: Equation (19) is valid as follows: , the following statements hold: (i) for a.e. w Î W there exists a limit, From these facts and the assumption of the completeness of X ¢ , P ¢ becomes an invariant measure with respect )  w Î´W. (i) Thus, according to Birkhoffʼs ergodic theorem, for P ¢-a.e. w Î W there exists a time average of g ¢ , i.e. g gT d t , which is equivalent to gT ( ) w defined in equation (22) (we abbreviate as g gT gT ¢ = º ). Note 'P¢-a.e.' can be replaced by P-a.e. or l N -a.e., since P ¢ , P, and l N become equivalent (i.e. absolutely continuous with each other) because , 0 r r ¢ > , so that we abbreviate them to just 'a.e.' (ii) It follows from the ergodicity that g g dP P ( ) ( ) due to equation (26). (iii) Since the l N -integrability of ρ is nothing but the P-integrability of g 1 º (constant function), the above conditions (i) and (ii) apply to Thus, the division of equation (23) by equation (27) leads to the validity of equation (24) for a.e. ω. , Hence, the accurate description of equation (19) is given in the context of equation (24), which now becomes g g a.e. , 28 In the current method, we thus realize the probability measure P P P ( ) º W and the probability space , under the assumption of the ergodicity and the integrabilities of ρ and r ¢ . This is done through the ability of calculating the expectation value, gdP ò W , for any function g L P 1 ( ) Î via the scaled-time average g . Namely, as is explicitly stated, the probability for any set B N  Î W can be defined and represented as in equation (30) is valid for the l N -a.e. initial point w Î W, and the set of all such points should be written as B W considering the dependence on uncountably many Bs. Thus it Î W . Similar notification applies hereafter). Therefore, using the equivalence of P and l N , the completion of , , Remark 2. We see that proposition 1 holds for any C 2 , strictly-positive :  k W  , instead of T . This means that the probability-space description in a usual ergodic dynamical system, which is generated by a vector field having a smooth invariant density ρ, applies for a dynamical system that is described by a scalar-field-multiplied vector-field along with keeping the density ρ, as long as we consider the scaled-time average. Of course 1 k = reduces to the usual description. Thus we have considered a general situation to obtain a probability-space structure.
p =´´´´´" Î whose rhs implies the marginal distribution. Marginal distributions for other components are given similarly.

Realization of a joint distribution of physical and parameter systems
In section 2.2, the phase-space density, ρ, has been treated in a very general manner, where we have only assumed the following mathematical conditions: smoothness, positivity, and integrability. To perform our purpose, namely, to create the joint distribution density, equation (8), we should set a suitable specific form of ρ along with satisfying the three conditions. This is established in the subsequent section 3.1. The fundamental results, including the creation of the joint distribution density, will be demonstrated in section 3.2.
Regarding the function T, in this section it still retains its abstract form, simply with the condition of the smoothness and positivity, and we will discuss in section 4 its more specific form.

Density form
To create density (8), we set as follows: x p Q Q : , , , Each function constituting ρ and function T have to satisfy the following six conditions:  has a definite signature.
Note that σ is introduced to control the admissible space of β via mapping using variable Q that moves in a whole free space m  , which is convenient in equation (34) appears as the Jacobian deduced from this mapping, and it is unity if one uses mapping σ that is the identity, where Q b = . (C2) Note that the fundamental variables to perform our purpose are the physical variables x p , ( ) and the parameter variable Q, for which their density form should be fixed in view of equation (8). Thus we choose a simple and general expression of the density form regarding the remaining variables , ζ, and η, which seem to be subsidiary but control the fundamental variables in total.
Note that (C4ii) seems to be a natural condition, but will be required only for ensuring equation (55) in the following contexts.
Although (C6i) has already been assumed, we state it here for clarity. We show that ρ defined above satisfies the required properties: , and used the integrability condition (C5) for the finiteness in equation (42c). Therefore we obtain the integrability of ρ and the integral formulae in equations (35a) and (35b).
The measurability of R r follows from the fact that D x p , , (40)) and the fact that this map is R r (equation (41)). The positivity of R r is evident. R r is finite for l n 2 -a.e., due to the integrability, indicated in equation (42c), through Fubiniʼs theorem. Similarly, we see that Z : , -measurable, finite for l m -a.e. β , and strictly positive. So does Z f ¢ , which equals Z (note , which leads to equation (35c). , As found from the above proof, without the integrability condition (C5), both R r and Z are well defined and Borel measurable, and (C5) is equivalent to both the finiteness of the integral in equations (35b) and that in equation (35c).

Fundamental results
The dynamics for flow T t { } ¢ , generated by the vector field X and the scalar function T, and the associated probability space P , , N (ˆ)  W W were formulated through proposition 1, where X was constructed from an abstract form of density function ρ which should satisfy just the normal three conditions. These conditions are satisfied for ρ which is specifically defined by equations (33) and (34) under (a part of) conditions (C1)-(C5), as shown in proposition 4.
In the following text, ρ is the density defined by equations (33) and (34), and T r r ¢ = . In this subsection, we first prove lemma 5 to show that the integrability of r ¢ required in proposition 1 is ensured by (C6). Next, we explicitly show that by using ρ we can realize the joint distribution density (8), that is, we can dynamically realize the superstatistics distribution. This is stated in theorem 6, as a special case of equation (24) ). Integral formulae for B are also shown in lemma 5. Consider the condition: Note, instead of B, we can treat b D : • p , for which similar procedures to the proof of proposition 4 will be used with careful attention to the indefinite signature of B. We have (39)), noting that this map is equal to Also assume (C7ii). Then we have dl n m , becomes l m -integrable, that is, So we can safely use a variable transformation formula to get , , a.e. , .   , the corresponding B satisfies (C7) due to the condition of T, that is, (C6), so that lemma 5 also b x p Q dl dl states , , . From these results and assumption (ii), proposition 1 with g b 1 • p º (this b is the one in equation (47) for a.e. ω. Lemma 5 shows the equality among b 1 • p and equations (43b)-(43d). + In this subsection we shall assume that the conditions (i) and (ii) in theorem 6 hold, so that propositions 1 and 4 are also applicable ( dl N ò r ¢ < ¥ W is concluded) as shown in the proof of theorem 6.
Theorem 6 shows that the joint distribution density, equation (   ere the integrability condition in (C8) can be replaced by a condition such that , , is l n m 2 + -integrable.
Proof. Note the two integrability conditions are equivalent: Equation (51) states that superstatistics distribution R r is realized by calculating the scaled long-time average of the physical quantities. It should be noted that for realizing R r the current scheme needs no explicit form of R r , which is defined by the integration formula (36). It is an advantage of the dynamical realization of superstatistics that R r is realized regardless of whether the integration can be explicitly done or not, as stated in section 1.
Note on ergodicity. So far, we have assumed the ergodicity with respect to T t { } ¢ and P , , Proof of the ergodicity is hard to achieve in general, and we here just note that the following two obstacles [8] of the ergodicity vanish in the current system: fixed points (i.e. For the first issue, we can easily observe from equations (15d) and (15f) and m 0 > that the fixed points do not exist. For the second issue, we have for any x p Q , , , , , Defining a map B : , , we have used an abbreviated form, A TRG Phys r r in equation (55), to represent the scaled long-time average of b : where we have used (C10) and (C4ii) in the last line to ensure the finiteness. Thus we can apply theorem 6 with B B R A º to obtain , is l n m 2 + -integrable, or equivalently, hfZ is integrable, we can apply theorem 6 with putting B B : , , but is an induced measure for m ¢ by a map σ, that is,

Example of function setting
For realizing the joint distribution of the physical and parameter systems in a simulation, we should set the specific forms of the functions: Phys r and Cntr r in equation (33), σ and f in equations (33) and (34), and T in equation (15). These functions should be compatible with (C1)-(C6). First, in section 4.1 we shall consider Phys r , Cntr r , and T. For σ and f, accordingly, their definitions affect only themselves as described below, and they will be discussed in section 4.2.

Fundamental functions
We consider Phys r , Cntr r , and T, while just assuming (C1) and (C4i) for σ and f, respectively.
Here we concentrate on an important case such that x p E x p , , , , ) is the energy of the physical system. Namely, Phys r depends on x and p only through the energy function E. This indicates the situation where the physical state density is determined by the physical-system energy E x p , ( ) and the parameter β, for which they are independently given. We consider that this corresponds to the situations (I) and (II) introduced in section 1. While for situation (III), where any other system is considered as well as the original physical system, the specification in the form of x p E x p , , , , The specific details of the functions are given as follows: , with M being a symmetric, positive-definite square matrix of size n (over  ); O : We see that (C2), (C3), and (C6i) hold. In molecular simulation, potential energy function U is given as the force field and K is the kinetic energy with M such as m m diag , , n 1 ( )  . It follows from (S1) and (S3) that (C6ii) is equivalent to l n m 2 + -integrability of a function , , , ). It should be noted that insightful ideas regarding a general temperature function are seen in [37,38]. Now, applying (S1) and (S2) into equation (33) and substituting the resultant density and temperature-function (S3) into equation (15), we have the following EOM: x (S1-1) : requires the value of U x ( ) in general through the evaluation of F x p Q , , T ( ) (and t )in a general case), the above map ordering ensures that the evaluations of the potential energy U x ( ) and the force F x ( ) is once in this integrator mapping. Maps, t

[ ]
F , at both ends do not need this kind of evaluation but just refer to the values. Note that S 1 and S 2 can be used for an EOM that has a more generalized form than equation (60)); and equation (73) can be extended to apply the scheme into such a generalized form.

Numerical test
To verify the present method, we performed a numerical simulation using the function setting described in section 4 with . This provides the simplest model system, where a onedimensional harmonic oscillator (1HO) is in a heat bath whose temperature is fluctuating obeying the distribution density f f G º (namely, type (I), where the physical system is defined by the 1HO and the parameter system is the temperature system). Note that the 1HO is suitable to validate the method because it gives a typical and important model that describes the physical system behavior around the equilibrium state and because it is simple enough such that the exact distribution is analytically obtained. There also includes a complexity in that the BG distribution is not trivially achieved via the sampling by use of the single NH equations due to the lack of ergodicity [40].
Our main purpose of the simulation is to examine the accuracy of the generated distributions of both the physical and temperature variables. Since this method presents new dynamics including both the physical and temperature variables, we also observed the associated trajectories generated by the simulation. The correlations among the dynamical temperature, the kinetic energy and the potential energy were also investigated. 5.2.2. Results. We show numerical results to confirm the realization of the distribution, equation (33) with (S1)-(S5), via equation (60). Figure 1 shows the distribution densities of the coordinate x, the momentum p, and the control variable ζ for the physically-based system, and the inverse dynamical temperature Q ( ) b s = , the momentum , and the control variable η for the temperature system. The simulated and exact results agree well, as was seen in small standard deviations of the discrepancies from the exact distributions for x, p, ζ, β, , and η, being 2.3 10 3 -, 2.2 10 3 -, 3.6 10 3 -, 4.5 10 3 -, 2.4 10 3 -, and 3.2 10 3 -, respectively. These results indicate that the present method produced a sufficiently accurate distribution. The reweighted result to the BG distributions was also satisfactory, where the standard deviations of the discrepancies from the exact distributions for x and p were 3.1 10 3 and 3.0 10 3 -, respectively.
The behavior of the physical system in the nonequilibrium environment provided by the temperature system is shown in figure 2(a) and (b). The time developments of physical temperature, T K p n p 2 P 2 ( ) º = , are perturbed from the pure 1HO Hamiltonian system, which should have period π for p 2 . Compared with T P , the dynamical temperature, , changes more slowly ( figure 2(a)). Around the time when T D is large in oscillating-like motions, the amplitudes of T P also become larger in an interesting manner. As shown in figures 2(d)-(f), ζ, , and η show oscillation-like behaviors in a nonlinear regime. As expected from the EOM, T T P Ḋ z = -(the 3rd equation in (60)), figure 2(d) shows that ζ mainly obeys T P , which shows faster oscillating motion than that of T D . On the other hand,  mainly changes according to T D (figure 2(e)), as suggested by the relations, Q  = and T Q l exp -. Similar correspondence appears in η (figure 2(f)) and , via 1 2  h = -. Figure 2(b) shows a time development of physical potential energy U x ( ), suggesting a correlated motion with T P . In fact, the correlation with a retarded time is seen by a correlation function (figure 2(c)), defined by G t K pt U x t t t 0 0 0 ( ) ( ( )) ( ( )) å º + . These correlations must be due to the nonlinear feature of ODE (60) and help to surmount nonergodic behaviors seen in conventional methods [41]. Keeping constant in the trajectory of the invariant function (equation (70)) in figure 2(b) shows the success in the numerical integration.

Conclusion
We have developed a new formalism, DDD, to realize the densities of both the physical and parameter systems through the scaled long-time average generated by an ODE. The joint distribution density of the physical system and the parameter system is realized under the ergodic assumption of the flow, and the physical system density R r is realized without its explicit form (that is, the superstatistics distribution is dynamically realized). The functions employed here to specify the vector field are simply an example among many choices. DDD demonstrates its potential for being (i) fruitfully utilized for various problems via the creation of arbitrarily designed temperature functions ((S3)), (ii) exerted for well-planned parameters instead of temperature (equations (15) and (60)), and (iii) generalized for addressing multiple densities hierarchically for, e.g. multiscale problems (equation (12)). These new, general ideas are utilized for efficiently simulating a physical system in nonequilibrium and in various environments, and thus it aims to advance molecular simulation, which is now a standard tool for atomic, molecular, soft-matter, and biological physics. We have focused on a theoretical study, and extensive applications are now in progress.