Measuring Nonequilibrium Temperature of Forced Oscillators

The meaning of temperature in nonequilibrium thermodynamics is considered by using a forced harmonic oscillator in a heat bath, where we have two effective temperatures for the position and the momentum, respectively. We invent a concrete model of a thermometer to testify the validity of these different temperatures from the operational point of view. It is found that the measured temperature depends on a specific form of interaction between the system and a thermometer, which means the zeroth law of thermodynamics cannot be immediately extended to nonequilibrium cases.


I. INTRODUCTION
Temperature and entropy are basic concepts of thermodynamics which have clear definitions and meaning in equilibrium but which are not yet fully understood in nonequilibrium situations. In equilibrium thermodynamics, one way to introduce temperature is to define entropy somehow (e.g. through the adiabatic invariant) so that temperature can be introduced as a derivative of entropy with respect to energy; However, since nonequilibrium entropy has never been constructed in a consistent way, we cannot define nonequilibrium temperature in this manner. For example, we don't know whether the entropy is a measurable quantity in contrast to equilibrium cases where entropy difference between two states is measurable by heat produced in quasistatic processes. Furthermore, even if we can measure nonequilibrium entropy, we cannot obtain unique temperature unless we properly set up the thermodynamic state space [1,2,3,4]; i.e. the value of the temperature depends on the choice of variables which we will fix through the differentiation of entropy with respect to energy [1]. Local equilibrium temperature, which we are familiar with, looses its validity for systems where the deviation from equilibrium ensemble is not negligible. Indeed, it is expected that equipartition of energy will no longer be valid, in such a way that different degrees of freedom may have different energy. For instance, some numerical simulations showed that nonequilibrium systems are anisotropic regarding with their kinetic energy [5], which is never explained based on the local equilibrium assumption.
So far, several authors have tried to seek the meaning of temperature beyond the local equilibrium picture utilizing the microscopic expression devised by Rugh [6].
where N is the number of degrees of freedom, m denotes mass of the microscopic particles, and φ is the interparticle potential. Note that we need arbitrary factors a and b to let the dimensions of the two terms (both of the numerator and of the denominator) be the same. We remark that, however, this arbitrariness has no influence on the value of temperature in equilibrium situations. Although Eq. (2) is originally defined in microcanonical ensemble, Jepps et al. [7] generalized this expression for canonical ensemble and presented it in a more general form. Furthermore, they applied it to numerical simulations of nonequilibrium stationary states in the presence of shear flow or heat flow. However, in nonequilibrium systems, we cannot have the unique value of temperature due to the arbitrariness of a and b [8]. This ambiguity seems quite natural since Rugh's expression is essentially based on equilibrium thermodynamic relation Eq. (1). Namely, the problem is carried over from the choice of variables to be fixed: we cannot reach the dynamical expression of nonequilibrium temperature unless the correct thermodynamic state space is set up. There is another way to define temperature which we call operational temperature in this paper. When finite closed systems are in contact, they finally equilibrate to have the same intensive quantity, which we identify with temperature (i.e. the zeroth law of thermodynamics). Hence, it might be possible to measure nonequilibrium temperature by putting equilibrium thermometer in contact to nonequilibrium systems. A gedanken experiment has been proposed in order to clarify the meaning of nonequilibrium temperature from that point of view [9]. In particular, Baranyai has performed numerical simulations on shear-flow or heat-flow systems in contact with thermometers and obtained some explicit values of operational temperature [5]. However, since we don't have a theoretical framework in which the obtained values should be interpreted, those numerical data seem to be left alone. In other words, we cannot theoretically predict the value of operational temperature with which the numerical data should be compared.
Since the meaning of temperature out of equilibrium seems to lack a sound theoretical basis, it is reasonable to pick up the simplest model in order to analyze theoretical aspects more easily. For that purpose, it seems that systems such as shear flow or heat flow are still excessively complicated. In this paper, we adopt a forced harmonic oscillator as a model system. Although it might be regarded as one of the simplest nonequilibrium systems, it is worth noting that time-averaged distribution functions of momentum and position are both Gaussian but characterized by different effective temperatures depending on the forcing frequency [10]. Thus, this system may provide a simpler example for the concept of temperatures than fluid systems under shear where different effective temperatures have been studied so far. We will simulate an experiment analogous to [9], by letting a forced oscillator interact with another nonforced oscillator (in a different heat bath) which acts as a thermometer. Comparison of the respective results may be useful for clarification of the concept of temperature.
The plan of the paper is as follows. First, some statistical properties of a forced harmonic oscillator in a thermal bath are recalled and interpreted in terms of nonequilibrium temperature. In the second section, in order to define temperature in a macroscopic point of view, a forced and an unforced oscillators situated in different thermal baths will be considered and the heat current between them will be calculated. In the final section we discuss the form of the entropy for some different choices of variables, and compare our result with those obtained by Baranyai in the framework of nonequilibrium molecular dynamics of fluid systems of soft spheres in shear flow.

II. A MODEL SYSTEM: FORCED HARMONIC OSCILLATOR
We assume that our model system is described by the following Langevin equation; where the mass of the oscillator is taken as unit. The natural frequency of the harmonic oscillator is denoted by Ω, and A sin ωt corresponds to external forcing. The noise term ξ(t) is assumed to be Gaussian white noise which satisfies where β is the inverse temperature of the heat bath. Macroscopic or thermodynamic quantities should be defined by an appropriate averaging; usually ensemble-average or time-average. Throughout this paper, we will take time-averaged quantities as macroscopic variables, since the model system is periodic in time due to sinusoidal forcing.
In our model, the external force gives power input to the oscillator, which may cause different influence on the average energy of momentum and of position. To have a physical idea of this influence, note that we have two kinds of relaxation times, each of which is related to position and momentum, respectively. We write the relaxation time of position as τ x = γΩ −2 and the one of momentum as τ p = γ −1 . When the forcing period τ A = 2π/ω is longer enough than a relaxation time (i.e. τ A ≫ τ x or τ p ), the sinusoidal motion is averaged to yield a distribution function which deviates from the equilibrium one. In contrast, if the forcing period is comparable with (or shorter than) a relaxation time (i.e. τ A ≤ τ x or τ p ), the corresponding motion of position or momentum cannot follow the forcing and the distribution function is indistinguishable from the equilibrium one. For instance, if τ p ≪ τ A ≤ τ x , we may expect that the distribution function of position is not much changed from equilibrium, whereas the one of momentum is modified by the forcing.
To see these circumstances explicitly, we will calculate the potential energy u con and the kinetic energy u kin . Since the dynamics is linear, we can separate the ensemble-averaged motion and the fluctuation. After some textbook-like calculations we get where The corresponding energy dissipation rateẇ iṡ Then we wish to consider the distribution function described by the Kramers equation corresponding to our model The solution independent of the initial condition is where the time axis is shifted suitably. Since this solution corresponds to ensemble distribution at a given instant, in order to calculate time-averaged quantities, the distribution function Eq. (11) itself must be time-averaged, To perform this integral, we expand the integrand in Eq. (12) up to second order in βpωX 0 and βΩ 2 xX 0 (Gaussian approximation), providing that β(ωX 0 ) 2 ≪ 1 and β(Ω 2 X 0 ) 2 ≪ 1. Namely, we assume that gain of the internal energy due to the external forcing is small compared with thermal energy. With this approximation, we get the following distribution function after simple calculations. where We can see there are two kinds of temperature for one oscillator. Hereafter each temperature corresponding to position and momentum is called configurational temperature and kinetic temperature, respectively.

III. OPERATIONAL TEMPERATURE
The discussion in the previous section deals with the distribution function in the phase space and hence it may be rather microscopic consideration. From the thermodynamic point of view, the problem arises how the different microscopic temperatures of Eqs. (14) and (15) are connected to macroscopic measurements. Since those temperatures differ, the macroscopically measured temperature (i.e. operational temperature) may indicate different values depending on the details of the connection between the system and the thermometer. In this section, we will investigate the problem by devising a concrete model for the temperature measurement.
Here we will examine a situation which bears some similarities with the proposal by Jou and Casas-Vázquez [9]. They defined a prototype of operational temperature in which two systems are in thermal contact: the one is kept in a nonequilibrium steady state by means of heat flux, whereas the other is in equilibrium to act as a thermometer.
In our setting, we consider two coupled oscillators in contact with different heat baths whose temperature can be controlled independently. The one is forced to stay away from equilibrium (the system) while the other is unforced to remain in equilibrium (the thermometer). They are connected through a weak interaction potential V . The schematic picture of our situation is shown in Fig. 1.
The dynamics of such oscillators will be written as where α denotes the mass of unforced oscillator. We let α and ǫ to be small so that the disturbance of the forced oscillator by the thermometer will be weak. The noise terms ξ i (t) are again assumed to be Gaussian white noise but characterized by different temperatures denoted by β −1 i , i.e.
where δ ij is the Kronecker's delta (unity if i = j and zero otherwise). The similar system, in the absence of inertia and the external forcing, was studied in detail by Sekimoto [11]. We will test two kinds of interaction terms; harmonic and bistable potentials.
The heat flux from the forced oscillator to the thermometer is evaluated as When both oscillators are left unforced (A = 0), the heat flux between both systems is proportional to the difference of the temperatures of the corresponding heat baths. When one of the oscillators is forced (A = 0), the unforced oscillator plays the role of thermometer. Our definition of the operational temperature is as follows. We fix the parameters of the forced oscillator (i.e. A, ω, β 1 and Ω), and change the temperature of the heat bath for the thermometer (β 2 ). There should be a certain value of β 2 at which the average heat flux q vanishes. Then we identify the temperature of heat bath β −1 2 with the operational temperature. Throughout the numerical simulations, we set γ = 1.0, A = 1.0, β 1 = 1.0, ǫ = 0.1, and α = 0.1.
The results of the numerical simulations are shown in Figs. 2 and 3 for different parameters of the oscillator, where the values of q are rescaled suitably. Note that the zero-point of heat flux (equilibration point) is different depending on the interaction potential: thermometers indicate different values for the same system. Especially, the one with harmonic potential shows good agreement with the configurational temperature, while the other with bistable potential indicates a value close to the kinetic temperature. This tendency is unchanged for the system with different parameters at which the configurational temperature is higher than the kinetic one (Fig. 3), while the latter is higher than the former in Fig. 2.  In short, the operational temperature of a thermometer with harmonic potential is almost the configurational temperature, and the one with bistable potential is close to the kinetic temperature. It is concluded that various interactions show various temperatures whose values are ranged between kinetic and configurational temperatures. Therefore, it is plausible that the zeroth law of thermodynamics should be formulated only in a very restrictive form for nonequilibrium cases.

IV. CONCLUDING REMARKS
In this concluding section, we will pay some attention to entropy corresponding to the nonequilibrium system, and to a comparison with the results obtained by Baranyai and coworkers by means of the nonequilibrium molecular dynamics.

A. Entropy for nonequilibrium steady states
In this subsection we will see some possible extended Gibbs relations of forced harmonic oscillators. The motivation is to discuss temperature in thermodynamic point of view and to compare it with the operational temperature obtained in the previous section. The distribution function Eq. (13) based on the Gaussian approximation may be considered to give the maximum entropy with the constraint that the second moments of momentum and position are given.
The corresponding Gibbs relation will thus be ds = β con du con + β kin du kin , where s denotes entropy per oscillator. This may be a natural extension of the equilibrium Gibbs relation by considering that each degree of freedom is a thermodynamic system by itself with different temperature. This simple form of the Gibbs relation is due to the Gaussian approximation of the distribution function. If they are not Gaussian, we would need additional independent variables; e.g. higher order moments. Recall, anyway, that systems with two temperature are common in nonequilibrium physics, for example in plasmas (where electrons and ions may have different temperatures), or in metals or semiconductors (where electrons may exhibit a temperature different from that of the lattice). However, Eq. (21) is not the only candidate. We can consider arbitrary linear transformation of u con and u kin such that u = u con + u kin , where α 1 = α 2 . By using these new variables, we can rewrite Eq. (21) as Since u is the total energy of the oscillator, it may be possible to define temperature analogous to Eq. (1); which yields Here we use another notation θ as the (inverse) temperature of the nonequilibrium system, which clearly depends on the choice of α 1 and α 2 . Namely, the thermodynamic temperature θ depends on the choice of the new variable y.
Although the nonuniqueness of temperature has been argued, e.g. in the context of extended irreversible thermodynamics [1], the criteria for the choice of new variables is still unknown. At least there are some necessary conditions for the choice of y: • it must be extensive, and macroscopically observable.
• the entropy must be convex regarding with the new extensive variable y.
In Eq. (21), the convexity is identical with that β con and β kin are nonincreasing functions of u con and u kin , respectively. This is obvious since β con = (2u con ) −1 and β kin = (2u kin ) −1 . Of course there are other choices satisfying convexity. Say we set α 1 = 1 and α 2 = 0, then where the thermodynamic temperature coincides with the kinetic temperature. Another possibility is to adopt the entropy production rate σ where Eq. (9) is recalled. Since dσ = 2γβdu kin , from Eq. (21) we get where the thermodynamic temperature becomes the configurational temperature. Using Eqs. (14) and (15), we can further rewrite Eq. (29) as where we can see the second-order contribution of the flux σ to the entropy. Note that the above expression with energy and entropy production rate are analogous to that in extended irreversible thermodynamics [1], where the usual thermodynamic variables and the fluxes are taken as independent variables. In this case we have β con as the thermodynamic temperature. In addition, we remark that the second terms of Eqs. (27) and (29) vanish at resonance where ω = Ω. It is identical with Also the entropy production rate σ is maximum at resonance. As we have seen so far, thermodynamic temperature defined through extended Gibbs relation depends on the choice of the new variable. We have expected that the operational temperature would be the criterion for choosing the new variable, which was one of the motivations of our study. However, in the previous section, we have seen that different thermometers read different temperatures, which means impossibility of defining the unique nonequilibrium temperature even in this simple model system consisting of one degrees of freedom. The absence of a unique operational temperature can be a serious problem for the construction of thermodynamics: at least the formulation of the zeroth law is not immediate and, if possible, would be a very restrictive form in contrast to equilibrium cases. Sasa and Tasaki have already pointed out this kind of operational restriction which results from the anisotropy of pressure in a macroscopic heat conducting system [12].

B. Comparison with the results of nonequilibrium molecular dynamics
In this paper we have examined kinetic, configurational, and operational temperatures in a forced harmonic oscillator. As was mentioned in the introduction, similar situations have been examined by Baranyai [5] based on the same motivation. While he studied systems consisting of soft spheres under shear flow or in the presence of heat current using techniques of nonequilibrium molecular dynamics, the system analyzed here is much simpler than that. It must be noted that Baranyai used the Nosé-Hoover type dynamics which removes energy from the system as dissipation, whereas we adopt the Langevin equation to represent the effect of heat baths. Despite these differences, we think that it is still worth comparing these results since we are looking for general thermodynamic concepts which should be largely independent of the microscopic details of the system.
The definitions of temperature in the works of Baranyai is based on the Rugh's microscopic expression of Eq. (2). As we have mentioned in the introduction, the expression itself cannot define temperature uniquely in nonequilibrium states due to arbitrary factors a and b. Instead of the unique temperature, again we have two kinds of temperature; configurational and kinetic temperatures which can be defined without the arbitrariness.
Theses expressions yield β con = (2u con ) −1 and β kin = (2u kin ) −1 which coincide with our results obtained by Gaussian approximation. Baranyai calculated these expression for temperatures in the mentioned fluid system of soft spheres at various values of the shear rate, and found that the configurational temperature is higher than the kinetic temperature, whereas the situation is opposite in systems with charge current. Furthermore, these temperatures turned out to be anisotropic: i.e. they take different values for different spatial directions. In our situation, the relation between the configurational and the kinetic temperatures depends on the ratio Ω/ω. In general, as Baranyai has discussed, they will depend on the characteristics of the system and the external forcing responsible for the nonequilibrium situation. In addition, Baranyai has proposed an operational temperature by devising the concrete model which emulates a physical thermometer in contact with the fluid. The thermometer consists of the same particles as the fluid's, but do not feel the effect of shear flow nor the thermostatting: they interact only with the fluid particles. This thermometer seems to read definite values of temperature regardless of the mass and the number of thermometer particles. However, the dependency of interaction potential between the fluid and the thermometer is not discussed. Taking our result into consideration, the operational temperature will depend on the interaction between the system and thermometers. Indeed, Hoover et al. have discussed the ideal gas thermometer which reads the kinetic temperature [13], while Baranyai's thermometer reads the value which is closer to the configurational temperature. However, since Hoover et al. ignored the anisotropy of the kinetic energy, it is not apparent what value the ideal gas thermometer reads when it is actually in contact with a nonequilibrium system. More numerical simulations and real experiments on the operational temperature are needed for the clarification of nonequilibrium temperature and the underlying nonequilibrium thermodynamics.
In summary, we have found that the operational temperature depends on the details of interaction between the system and a thermometer, which is never seen in equilibrium situations. The fact may shove a strong restriction on the extension of the zeroth law for nonequilibrium systems. Of course, as our analysis is confined to one-dimensional systems, there may arise another problem for two-or three-dimensional systems regarding with the relation between the anisotropy of temperatures and the operational temperature. In addition, the Gaussian approximation of the distribution functions make especially easy to define temperatures, while additional conceptual problems would appear if the distribution function deviates far from the Gaussian. For instance, there have been some maximum-entropy analyses of nonequilibrium radiation, where nonequilibrium temperature or quasitemperature per mode have been defined in terms of the nonequilibrium populations of the different modes, in the context of a generalized Planck statistics instead of a classical Boltzmann statistics [14]. The present situation has the advantages of the higher simplicity and of the possibility of devising numerical simulations concerning operational temperature which have not been done before.