Violation of the virial theorem and generalized equipartition theorem for logarithmic oscillators serving as a thermostat

A logarithmic oscillator has been proposed to serve as a thermostat recently since it has a peculiar property of infinite heat capacity according to the virial theorem. In order to examine its feasibility in numerical simulations, a modified logarithmic potential has been applied in previous studies to eliminate the singularity at the origin. The role played by the modification has been elucidated in the present study. We argue that the virial theorem is practically violated in finite-time simulations of the modified log-oscillator illustrated by a linear dependence of kinetic temperature on energy. Furthermore, as far as a thermalized log-oscillator is concerned, our calculation based on the canonical ensemble average shows that the generalized equipartition theorem is broken if the temperature is higher than a critical temperature. Finally, we show that log-oscillators fail to serve as thermostats for their incapability of maintaining a nonequilibrium steady state even though their energy is appropriately assigned.

Theoretical and numerical investigations of statistical mechanical systems invariably rely upon a suitable modeling of thermostats. Two commonly used models in literature are Langevin thermostat and Nosé-Hoover thermostat (see, e.g., refs 1, 2 and references therein). The former involves in cumbersome calculations of stochastic processes and the latter is a deterministic bath with time-reversible dynamics, which however can give rise to irreversible dissipative behavior. A new model of thermostat, capable of producing controlled dynamics to sample canonical ensemble, is always desirable from both theoretical and numerical viewpoints.
Recently, a logarithmic oscillator (log-oscillator) has been suggested to be a candidate to serve as a Hamiltonian thermostat 3 . In the simplest 1D version, the Hamiltonian of a log-oscillator reads , indicating an infinite heat capacity 3,5 . This implies that its kinetic temperature can maintain a constant value even though the log-oscillator couples with another system, which absorbs energy from the log-oscillator or vice versa. This peculiar property make the log-oscillator an ideal candidate of thermostats.
Although the log-oscillator has the outstanding property, its feasibility as a thermostat has been argued recently [6][7][8][9] . The arguments mainly lie in its long length and time scales involved in practical implementations or numerical simulations. It has particularly been shown that the log-oscillator is incapable of establishing the heat flow in a nonequilibrium system 8 and fails to serve as a thermostat for small atomic clusters 9 .
Actually, it has been known for a long time that Brownian particles with momentum-dependent drift can be described by power-law tail distributions and correspondingly show anomalous diffusion in optical lattices [10][11][12][13][14][15][16][17][18] . In this case, systems are equivalently described by Brownian particles confined inside a logarithmic potential in momentum space. Ergodicity of such systems is accordingly broken 10,11 and the Brownian motion approaches to an infinite covariant density 12 . Distinct phases of the anomous diffusion and nonstationary scale-invariant dynamics have been extensively investigated 11,[13][14][15][16] . The phase-space density can deviate from the Boltzmann-Gibbs statistics at the intermediate energy and exhibit heavy power-law tails 17,18 .
Note that in order to avoid the singularity at the origin, a common way in previous studies is to apply a modified form of the logarithmic potential (see Eq. (3) below) when one performs numerical simulations 3,5,8,9 . Some controversial conclusions on the original log-oscillator (1) given in previous studies were drawn from the simulation results based on the modified log-oscillator. However, less attention has been paid to the effects coming from the modification required in numerical studies. It is therefore necessary to elucidate the role played by the modification of the logarithmic potential and further understand the underlying mechanism for the failure of modified log-thermostat in numerical simulations. In the present study, we study the statistical property of the modified log-oscillator from both microcanonical and canonical aspects. Interestingly, we find that the virial theorem is seemingly violated, illustrated by an unexpected energy dependence of kinetic temperature numerically. Furthermore, as a thermalized log-oscillator is concerned, our analysis shows the generalized equipartition theorem is broken when the temperature is higher than a critical temperature, which is verified by numerical simulations. Finally, we show that modified log-oscillators cannot serve as thermostats to produce a stationary temperature profile inside the FPU-β system even though the energy of log-oscillators is appropriately assigned.

Results and Discussions
Isolated log-oscillator: Violation of the virial theorem. For the sake of avoiding the singularity of the logarithmic potential at x = 0 particularly when one performs numerical simulations, as mentioned above, an usual way is replacing the logrithmic potential of Eq. (1) with the following modified potential 3,5,8,9 : One can see that the presence of the parameter ε eliminates the singularity and make numerical simulations feasible. However, some peculiar effects come from the presence of ε as shown below. The equation of motion is solved by using the velocity-Verlet algorithm and a time step of 0.0001. This algorithm is symplectic and energy conservation is numerically satisfied in our simulations. In what follows, m, k B and b are set equal to unity unless otherwise stated. The results of numerical simulation shown in this section implying that the small parameter ε lead to violation of the virial theorem for the log-oscillator.
By defining kinetic temperature ≡ T k k p m 1 B 2 , the virial theorem for the modified system (3) yields where 〈·〉 refers to a long-time average with respect to a dynamical trajectory.
For an ideal log-oscillator, kinetic temperature T k should be energy independent according to the virial theorem. As a small ε is introduced in a way like Eq. (3), energy dependence of T k is notable only when energy is very low and comes to vanish as energy increases. Surprisingly, we observe an anomalous energy dependence of T k which seemingly deviates from the prediction from the virial theorem when energy of the system is large. As shown in Fig. 1, we plot the dependence of kinetic temperature T k on the total energy E of the system. One can see that there is a plateau region for a certain range of energy, indicating that the kinetic temperature of log-oscillator is equal to τ and independent of the energy of system. In this region, the log-oscillator nearly has an infinite heat capacity. However, T k increases linearly with E as E > E c , where the critical energy E c is independent of the magnitude of ε. Our numerical studies show that the plots coincide for different value of ε, even for the case that ε is as small as 10 −8 . For the increasing part, linear fitting gives that T k = 2(E − E c ), which indicates E c = 〈V log 〉. The fitting relation T k = 2(E − E c ) simply comes from the energy conservation of the Hamiltonian system. Unlike a "normal" thermodynamical system, for which the potential energy quickly obtains the increment of energy while the kinetic energy is kept constant in average, the log-oscillator has a slow relaxation of energy. This explains why the kinetic temperature is proportional to the energy since the potential energy in average gains little the energy increment while kinetic energy gains almost all of them. Note that when energy E is small, as shown in the inset of Fig. 1, T k deviates from the value of τ while the virial theorem, i.e., T k = P k , is valid in this case as expected. The deviation is more considerable as ε increases. The presence of ε leads to the correction in the density of states in the low energy regime (see, e.g., Fig. 2 in ref. 3). To reduce the deviation, it has been argued that the total energy E should be large 3,6 . However, our results indicates that total energy should not be "too large", i.e., E should be less than E c . Otherwise, the modification of logarithmic potential seemingly indicates a violation of the virial theorem, which is required for the log-oscillator to serve as a thermostat.
We also study the effect of average time on the behavior of T k , as shown in Fig. 2. One can see that the plateau region is broaden with the increase of average time. However, the broaden process is extremely slow. Note that the critical energy E c shifts roughly a constant amount through increasing the simulation time by a factor of 10, which indicates an exponential dependence of the required simulation time on energy. As the integration goes to infinity, it is unrealistic through numerics to conclude that whether the broaden process would stop at some critical energy (similar as the critical temperature given in the section below) or keep ongoing, which requires a solid analytical study. However, since the model described by Eq. (3) is proposed mainly for serving as a thermostat for numerical simulations, it is safely to argue that virial theorem is practically "violated" for the modified log-oscillator as E > E c , which results in its failure to serve as a thermostat.
In order to understand the underlying mechanism further, we study the relaxation process of T k by setting the log-oscillator's total energy inside the plateau region (E = 5) and the linear increase region (E = 25), respectively. Meanwhile, we choose two typical trajectories for which the log-oscillator is located either near the origin (trajectory 1) or far from origin (trajectory 2) initially. For the case E = 5, T k relaxes quickly to the predicted value in terms of the virial theorem ( Fig. 3(a)). However, for the case E = 25, the relaxation process is extremely slow shown by our long-time simulations up to 10 12 time steps (Fig. 3(b)). Its relaxation to the equilibrium state is unavailable through practical simulations. The slow relaxation lies in the exponential growing of oscillation time with the energy of the log-oscillator. Therefore, in a much shorter timescale attainable in numerical simulations, the log-oscillator essentially moves with a constant velocity, which explains the linear behavior with energy shown in Figs 1 and 2.
As to make sure the violation does not originate from numerical errors, we also check our results by two higher-order symplectic algorithms, namely, Omelyan-Mryglod-Folk algorithm 19 and Laskar-Robutel algorithm 20 . The two algorithms have four-order and ten-order symplectic integrators, respectively. They recover the results by the velocity-Verlet algorithm with the total-energy drift as small as 10 −10 in a trillion-timestep run, namely, the accuracy is increased by three to four orders of magnitude in comparison with the Verlet algorithm.
Here H is the Hamiltonian of a system, x i,j stand for canonically conjugate (generalized) coordinates or (generalized) momentums. In this section, 〈·〉 denotes the canonical ensemble average at equilibrium temperature T. A special case of GET is the well-known equipartition theorem for a system of quadratic terms only. GET ranks among the basic results of statistical mechanics and has been extensively studied until today 22,23 . Note that although GET is similar with the virial theorem (see, e.g., Eqs (2) or (4) where β = 1/k B T, Z p and Z x are reduced participation functions with respect to p and x, respectively. In the section above we show that the virial theorem is practically violated for the log-oscillator when we replace the logarithmic potential by Eq. (3). In this section we will analytically show that GET is also broken for the modified log-oscillator and verify the result by numerical simulations.
Calculations shows that there exists a critical temperature T = τ for the integration above, for which one have P μ = k B T is valid only for T < τ (see details in the section of methods). P μ deviates from the linear behavior and approaches to an upper bound k B τ for T ≥ τ. This means that GET with respect to coordinate x is broken.
In order to verify the above analytical results, we perform numerical simulations by thermalizing the log-oscillator with a Langevin thermostat at temperature T. The equation of motion is given by log where γ is the friction constant and η denotes a Gaussian white noise with zero mean and variance of 2γk B T. Equation (8) is integrated by using a stochastic velocity-Verlet algorithm. Firstly, we study the dynamical behavior of K μ and P μ for some typical equilibrium temperatures, as shown in Fig. 4, where a canonical ensemble average is taken for any given time. For T = 0.1, after a quick relaxation process, K μ and P μ saturates at the value predicted by GET, i.e., K μ = P μ = k B T. This is also the case for T = 0.5 with a slightly slower relaxation process. For T ≥ τ, while K μ keep satisfying GET, one cannot see a sign for P μ converging to GET prediction. One can see the transition clearly from Fig. 5, where we plot K μ and P μ as a function of thermostat temperature. Unlike K μ , P μ approaches to the upper bound k B τ as T increases, which is consistent with the analysis above. Note that the transition of P μ near T = τ is not sharp but smoothly curvy. The deviation from the theoretical prediction may come from finite number of canonical sampling and finite relaxation time in our simulations. Our simulations show that P μ will gradually close to the analytical value by increasing the relaxation time.
Failure of log-thermostats to maintain a nonequilibrium steady state. In this section, we apply the log-thermostat to see if it can drive the system into a nonequilibrium state. It has been shown that Hamiltonian thermostats, including the log-thermostat, fail to promote heat flow 8 . Although there is a linear energy dependence for T k when E > E c , it would be still interesting to learn if the log-oscillator can serve as a thermostat when its energy is appropriately set inside the plateau region (Fig. 1). It seems promising since T k is practically energy independent and the log-oscillator nearly has infinite heat capacity in this case.
As the usual setup for heat conduction, the system is sandwiched between two log-thermostats at different temperatures through a weak coupling. The whole Hamiltonian is given by The system of interest is a one-dimensional FPU-β chain of oscillators   where γ ± are the strengths of the couplings between the log-thermostats and the FPU-β system. In our simulations, the energy of the log-oscillators is initially assigned in the low-energy zone of their own plateau region, namely, E − = 1 and E + = 20. As a comparison, we also examine the case that the energy of respective log-oscillators are initially to the edge of the plateau. The initial energy of each FPU oscillator is randomly prepared around an average energy equal to 8. After a relaxation process, the system is coupled with the log-thermostats weakly. The evolution of the temperature profile is shown in Fig. 6. One can see that the log-thermostat can drive the system into a nonequilibrium state only for a short time, indicated by a linear temperature gradient. As the average time increases, the log-oscillators cannot maintain the temperature gradient. Our numerical results shows that the FPU system, except the boundary oscillators, reach an energy equilibration after a relaxation time, i.e., heat flow cannot be promoted inside the system. Note that the kinetic temperature of the thermostats is kept at T ± = τ ± for the initial setup that E − = 1 and E + = 20, as is expected. However, if the energy is set near the plateau region (E − = 40 and E + = 80), the log-oscillators can easily gain energy from the system and jump out of the plateau region. In this case, T ± are quite different with τ ± . Thus τ ± cannot be regarded as measures of the thermostat temperature. Inspired by the Nosé-Hoover chain thermostats 2 , we also applied thermostats of a chain of log-oscillators at different temperatures as to see if nonequilibrium steady states can be established. However, this again gives negative results even though the "chaoticity" of the thermostats is increased in this way. Our results indicate that log-thermostats fail to maintain a nonequilibrium steady state (stationary temperature gradient) even though the energy of log-oscillators are assigned in the plateau region.
In summary, a particular property of the log-oscillator is its infinite heat capacity in terms of the virial theorem, which make it a candidate to serve as a Hamiltonian thermostat. With this study, we demonstrated that kinetic temperature of the modified log-oscillator has a linear dependence on system energy if energy is larger than a critical energy, which indicates a possible violation of the virial theorem. The critical energy increases with the average time in a slow way, with which its infinite-time behavior cannot be concluded by numerics. Therefore, it would be safely to argue that the virial theorem is practically violated in finite-time numerical simulations. As a counterpart of the virial theorem, GET was further investigated for the thermalized log-oscillator with respect to the canonical distribution. Based on the canonical average, we found that there exists a critical temperature, above which GET is violated, namely, P μ approaches to an upper bound instead of increasing linearly with the increasing of temperature. Our numerical simulations are consistent with the analysis. Finally, we show that log-thermostats fail to maintain a nonequilibrium steady state although we assign the energy of log-oscillators inside the plateau region. Note that it has been shown that equipartition is not satisfied for a Brownian particle under both harmonic confinement in positional space and logrithmic confinement in momentum space 17,18 . The violation comes from the correction of the the Boltzmann-Gibbs density distribution for such a system. For T > τ, the distribution is unnormalizable 10 . In this sense, one would not expect GET to hold in the first place due to the absence of a normalizable distribution. However, although their model is slightly different with ours, it is still interesting to see that our analysis based on the canonical ensemble average also shows a critical behavior at temperature T = τ, which is consistent with the analysis 18 . Further studies of the deviation from Boltzmann-Gibbs statistics in different situations are thus interesting. The last equality in Eq. (13) is obtained by using the property of Gamma function, i.e., Γ(y + 1) = yΓ(y) for y equal to any positive real number. Therefore, GET with respect to coordinate x, i.e. P μ = k B T is valid for T < τ. For T ≥ τ, the participation function is non-normalizable and an analytical form for P μ is so far unavailable. However, P μ has an upper bound k B τ by noting that 0 < Δ < 1, which indicates that GET is invalid with respect to x provided T ≥ τ.