Global Thermodynamics for Heat Conduction Systems

We propose the concept of global temperature for spatially non-uniform heat conduction systems. With this novel quantity, we present an extended framework of thermodynamics for the whole system such that the fundamental relation of thermodynamics holds, which we call ``global thermodynamics'' for heat conduction systems. Associated with this global thermodynamics, we formulate a variational principle for determining thermodynamic properties of the liquid-gas phase coexistence in heat conduction, which corresponds to the natural extension of the Maxwell construction for equilibrium systems. We quantitatively predict that the temperature of the liquid-gas interface deviates from the equilibrium transition temperature. This result indicates that a super-cooled gas stably appears near the interface.


Introduction
The phenomena of liquids and gases close to equilibrium have been extensively studied for a long time. As a macroscopic universal theory for the phenomena, the hydrodynamics is believed to be wellestablished [1], and the connection of the hydrodynamics with the classical and quantum mechanics of atoms have been discussed for over a century [2,3]. Nevertheless, in this paper, we construct a new universal theory for thermodynamic properties in the linear response regime. There are two main messages. First, a new concept, global temperature, is found, with which a novel framework of global thermodynamics is constructed to describe the whole of non-uniform non-equilibrium systems with local equilibrium thermodynamics. This outcome provides a fresh viewpoint for the description of systems out of equilibrium. Second, this formulation provides non-trivial quantitative predictions. As an example, let us consider pure water under a pressure of p ex = 1.013 × 10 5 Pa, where two heat baths of temperature T 1 = 368.0K and T 2 = 378.0K are in contact with the sides of the system. See Fig. 1 for an illustration. Recall that the liquid-gas transition temperature is 373.1 K. Based on global thermodynamics, in this paper, we predict that the interface temperature of the liquid-gas coexistence is 368.3K. This result means that super-cooled gas stably appears in heat conduction, which may be tested in experiments.
In the remaining part of this introduction, we first present a brief summary of development in nonequilibrium statistical mechanics, confirming that the phenomenon described above has never been discussed by established theories. We then provide a review of extended frameworks of thermodynamics so that readers can understand how the global thermodynamics proposed in this paper is different from previous frameworks. At the end of the introduction, we summarize the achievements of this paper.

Non-equilibrium statistical mechanics
In the universal theory of macroscopic irreversible dynamics built by Onsager, a steady-state current is assumed to be linear in the thermodynamic forces [4]. The linear coefficients form a non-negative symmetric matrix, which is referred to as the Onsager matrix. This theory with a variational principle for determining a current and/or a thermodynamic force is called irreversible thermodynamics, which is valid near equilibrium [5]. Onsager theory can be interpreted as a universal theory for the dynamical fluctuations of thermodynamic variables, where the properties of the Onsager matrix are connected to the stability of the equilibrium state and the time-reversal symmetry in microscopic systems. In short, this fluctuation theory is represented by a simple stochastic process with a symmetry property [4,6]. In accordance with this theory, statistical mechanics of trajectories has been studied based on microscopic dynamics, which leads to the expression of probability densities, linear and non-linear response formulas, and macroscopic deterministic dynamics [2,3,7,8,9,10,11].
After these developments, in the last two decades, our understanding of phenomena related to thermodynamics has progressed greatly because of the following two reasons. First, the development of experimental techniques for the measurement and manipulation of biological molecular machines [12,13,14,15,16] naturally leads to the extension of thermodynamics to mesoscopic scales, which is now called stochastic thermodynamics [17,18,19]. Second, the newly discovered relations, such as the fluctuation theorem [20] and the Jarzynski equality [21], have emerged as universal [22,23,24,25,26,27], providing a new starting point for the re-organization of previous theories with greatly simplified derivations of the formulas [28,29,30,31,32]. These two developed directions, stochastic thermodynamics and new universal relations, are related to each other and have given rise to a new connection with large deviation theory [33,34,35,36,37], and information theory [38].
These results are quite useful once we fix the phenomena under study. However, they do not predict a qualitatively new phenomenon that has never been studied to date. Indeed, our prediction on the stabilization of super-cooled gas by heat conduction can be studied based on the modern development of non-equilibrium statistical mechanics [39], while the prediction itself has never been obtained within a framework of non-equilibrium statistical mechanics. This is because the non-equilibrium statistical mechanics essentially comes from the Onsager theory, which is the variational principle on the current and/or the thermodynamic force but not on the thermodynamic variables. For example, the minimum entropy production principle and other related variational principles may characterize the steady state in the linear response regime [40]. These principles may be interpreted as a variational principle in which the statistical ensemble in the linear response regime is obtained from minimization of the entropy production as a function of the probability measure [41]. Nevertheless, this variational principle is not directly related to the determination of thermodynamic variables.

Extended frameworks of thermodynamics
Equilibrium thermodynamics provides a unified description of thermodynamic properties of materials at equilibrium. It also formalizes the second law, which leads to a variational principle for determining equilibrium states [42,43]. The variational principle naturally suggests the law of fluctuation of thermodynamic variables [44], which is formulated as a large deviation theory. From this viewpoint, a framework of statistical mechanics may be constructed in a consistent manner with the fluctuation theory of thermodynamic variables [45].
Therefore, in order to establish a universal theory for thermodynamic properties out of equilibrium, it is natural to consider an extended framework of thermodynamics. A naive attempt is to extend the equilibrium fundamental relation of thermodynamics for the case of simple fluids, where F is the Helmholtz free energy and S is the entropy. Examples of such attempts can be seen in Refs. [46,47,48]. The heart of the problem for the extension is to confirm the two conditions: First, the theory is self-consistent and self-contained; second, new predictions specific to the extension are presented. However, the two conditions are not confirmed in many studies. The derivation of the extended thermodynamics from the microscopic theory, regardless of its importance, does not make sense unless the extended framework satisfies the above two conditions. Here, we do not give a complete review of previous studies on the extended framework of thermodynamics, but study heat conduction systems from the viewpoint of the extended framework of thermodynamics. Heat conduction is described as a spatially extended system where thermodynamic variables slowly vary in space. The local subsystems, which are small but macroscopic, are regarded as a local equilibrium state [1,5,49]. If we followed this standard description, we could not have a strong motivation to seek a thermodynamic framework because all thermodynamic properties can be calculated by the heat conduction equation and the local equilibrium thermodynamics. Nevertheless, we still have two possibilities for considering an extended framework of thermodynamics.
The first attempt is to go beyond the local equilibrium thermodynamics. In this approach, which may be adapted in Refs. [50,51], the equation of state for the local subsystems is modified to contain the influence of the heat flux. Since such a contribution is quite small in the linear response regime, the theory is not useful even if it is correct. Although this approach may still be effective for spatially homogeneous driving systems, in which intensive parameters may describe the balance of extensive variables [51,52,53,54,55], we do not deal with such systems in this paper.
In the second approach, we retain the standard description for a spatially extended system with local equilibrium thermodynamics. On this basis, we then seek a thermodynamic framework for the whole system. As an example, we first consider the extension of the second law, by which a state variable "entropy" is defined, and we then derive the thermodynamic relation which corresponds to an extension of the fundamental relation of thermodynamics. A concrete procedure of the first step was proposed by introducing the concept of excess heat [48,56,57]. This idea was also studied from semi-macroscopic and microscopic theories [58,59,60,61,62,63,64]. However, the fundamental relation of thermodynamics has not been derived from the extended entropy. One reason may be that an interesting phenomenon associated with the extended entropy was not addressed. Nevertheless, one can continue to seek a possible framework without considering certain phenomena. Indeed, for a specific model of heat conduction, the extended entropy was numerically estimated [65], which suggests that the extended entropy is close to the spatial integration of the local equilibrium entropy density. Since the extended entropy can be obtained in experiments, a natural question is "What is temperature ?". If the fundamental relation of thermodynamics is formulated, there should be a temperature satisfying it. We consider this question seriously, putting aside specific phenomena.

Summary of results
We first focus on single-phase systems (either liquid or gas) in heat conduction. By assuming local equilibrium thermodynamics, we define the entropy S and the Helmholtz free energy F in heat conduction by the spatial integration of the local density fields corresponding to these variables. The pressure field p is homogeneous in space. Then, the problem is to find a temperatureT satisfying the fundamental relation dF = −SdT − pdV (2) in the linear response regime. In §3, we solve this problem by definingT as the kinetic energy averaged over particles in the system. We call this temperature global temperature. It should be noted that the local temperature T (r), which depends on the position r, satisfies the thermodynamic relations for each point r. In contrast to the local relations, (2) is a global relation applied for the whole system as if the system is at equilibrium. We call such a thermodynamic framework global thermodynamics. This formulation is also interpreted as a mapping of each heat conduction system to an equilibrium system through the novel quantity of the global temperature,T . The formulation for a single-phase system is indeed derived from fluid dynamics with local equilibrium thermodynamics. That is, this formulation is interpreted as a different formulation for describing thermodynamic quantities in the linear response regime. The prediction by using global thermodynamics can also be predicted by hydrodynamics with local equilibrium thermodynamics, in principle. Here, we go one step further. We consider a thermodynamic phenomenon that cannot be described by the standard hydrodynamics with local equilibrium thermodynamics. This is the phenomenon of phase coexistence in heat conduction. In §4, we study this phenomenon and we show how existing theories are not appropriate for determining the thermodynamic properties. The essential point is that the condition for the connection of the two phases is outside of the local equilibrium thermodynamics.
We study the phase coexistence in heat conduction with the framework of global thermodynamics. We first formulate a variational principle for determining the local temperature of the liquid-gas interface. The idea is quite simple. Fixing the global temperature of the whole system, we naturally extend the variational principle for equilibrium systems to that for heat conduction systems. This idea was proposed in our previous paper [66]. Remarkably, by using the solution of the variational equation, in §5, we derive a universal relation among the interface temperature θ, the equilibrium transition temperature T c , the global temperatureT , and the mean temperature of the two heat baths T m , which is We call this relation the temperature relation. From the temperature relation, we find that the temperature of the liquid-gas interface deviates from the equilibrium transition temperature. That is, super-cooled gas stably appears near the interface in heat conduction. This is a qualitatively new phenomenon that has never been considered in previous studies. Since the steady state as determined by the variational principle is expressed as a function of the global temperature, the prediction of measurable quantities for given conditions is indirect. Thus, in §6, we re-express all quantities in the steady state in terms of the two temperatures of the heat baths. We present several formulas of the interface temperature by directly using measurable quantities. Furthermore, we illustrate examples of quantitative results for a van der Waals fluid and pure water.
For the steady state determined by the variational principle, we further develop thermodynamics for heat conduction systems with phase coexistence. First, in §7, we derive the fundamental relation associated with the Gibbs free energy G. At first sight, the result does not seem to contain nonequilibrium extensions. This is because G is not differentiable at the transition point at equilibrium. By performing a careful analysis near equilibrium, we find that the fundamental relation holds in an appropriate equilibrium limit. The heat capacity and the compressibility, which are singular at equilibrium, are also obtained as a regularized form while breaking the additivity. In §8, we derive the fundamental relation associated with the Helmholtz free energy F . Since the free energy F is defined for the coexistence phase in equilibrium cases, its non-equilibrium extension can be written as a perturbation from the equilibrium form. We derive this expression explicitly.
It should be noted that the quantitative prediction is made based on a fundamental assumption of the variational principle. As is often observed in universal theories, one may replace the fundamental assumption by another one. In §9, we formulate the theory starting with another assumption instead of the variational principle. For example, when we assume the fundamental relation of thermodynamics for the whole system by using the global temperature, we can derive the results of the variational principle. As another example, one may focus on how the volume change near the transition temperature at constant pressure exhibits the singularity in the equilibrium limit. Supposing the simplest form of a singularity, we can derive the results of the variational principle. These findings indicate that the theory itself possesses an elegant structure.
In §10, we extend the theory in §3 to cover steady states with an arbitrarily shaped container beyond the linear response regime, while restricting our attention to single-phase systems. Our theory quantitatively predicts a new relation among the global quantities in this setup.

Equilibrium thermodynamics
We consider a macroscopic material at equilibrium. As the simplest example, we focus on a simple fluid whose thermodynamic state is characterized by temperature T , the volume V of a container, and the amount of material N . For a system in contact with a heat bath of temperature T which may be controlled externally, there exists a state variable S(T, V, N ), called entropy, which satisfies the Clausius equality for infinitesimal quasi-static heat d Q from the heat bath. The infinitesimal change of internal energy U (T, V, N ) of the material is determined by which is referred to as the first law of thermodynamics. d W is the infinitesimal quasi-static work required in the infinitesimal change of the volume dV , which is given by The substitution of (4) and (6) into (5) leads to the fundamental relation of thermodynamics: From this expression, we find that it is useful to consider U as a function of (S, V, N ). Indeed, a single function U (S, V, N ) leads to all thermodynamic properties such as equation of state and heat capacity We then define chemical potential µ as Various thermodynamic functions equivalent to U (S, V, N ) can be defined by the Legendre transformation of U (S, V, N ): The fundamental relations associated with these functions are Substituting G(T, p, N ) = N G(T, p, 1) into (16), we find This form together with (16) leads to the Gibbs-Duhem relation The extensivity of thermodynamic quantity leads to the concept of density defined as the quantity per unit volume, similarly to particle density ρ = N/V . For instance, entropy density is defined as By using ρ = ρ(T, p), which is obtained from (8), one may consider the entropy density as a function of (T, p), which is expressed as following the convention in thermodynamics. Similarly to the entropy density, for any extensive quantity A, its density is defined as a(T, p) = a(T, ρ(T, p)).
We here consider free energy density g = G/V . Substituting G = gV into (17) and (18), we obtain sdT + ρdµ = dp, where u = U/V . From these, we have where dT , dρ and dp are not independent, because T , ρ and p are connected by the equation of state (8) such that Furthermore, we define extensive quantities per one particle, aŝ Note thatĝ is equivalent to the chemical potential µ. We then have whereφ is specific volumeφ = 1/ρ.

Setup of heat conduction systems
Throughout this paper, we consider a system of N -particles, which are packed in a rectangular container with lengths of side L x , L y and L z , as shown in Fig. 2(a). We ignore the effect of gravity. L x and L y are fixed throughout this paper. L z is fixed for a constant volume system, or not fixed at constant pressure. We study heat conduction states driven by the temperature difference between two heat baths. As schematically described in Fig. 2(b), a heat bath of temperature T 1 is attached to the left end (x = 0) and another heat bath of temperature T 2 to the right end (x = L x ). Other four boundaries are thermally insulating. We take T 1 ≤ T 2 without loss of generality, and (T 2 − T 1 )/T 2 is assumed to be so small that the system is in the linear response regime. In such an idealized steady state without any convection, the system is regarded as a one-dimensional system. Local states are homogeneous inside any section perpendicular to the x axis, and therefore, local thermodynamic quantities are considered as functions of x.
We introduce a dimensionless parameter ε that indicates the degree of non-equilibrium as where Ξ is the temperature difference and T m is the mean temperature of the heat baths, which are defined by When we focus on the linear response regime around an equilibrium state, we ignore the contribution of O(ε 2 ).

Local equilibrium thermodynamics
For such macroscopic non-equilibrium systems, the hypothesis of local equilibrium thermodynamics works well. That is, local thermodynamic quantities are assumed to satisfy local thermodynamic relations at each space and time [1,5]. Now, suppose that the temperature profile T (x), the density profile ρ(x) and the pressure p(x) are determined by experimental observation. Any local thermodynamic quantities, such as a(x) andâ(x), are expressed as where functions a(T, ρ) andâ(T, p) are determined in thermodynamics, as described in the previous subsection. According to (22), a(x) may be also written as a(x) = a(T (x), p(x)) via ρ(x) = ρ(T (x), p(x)).
Then, local equilibrium thermodynamics means which correspond to the relations (23), (24) and (25), respectively. For the quantities per one particle, we also have which correspond to the local version of the thermodynamic relations (28) and (29).

Global conditions for steady states
For steady state heat conduction, the local pressure p(x) satisfies for any x and x . Especially, for the system at the constant pressure p ex , holds for any x. These equalities for the local pressure may be regarded as global relations because they are not obtained in the local thermodynamics. Furthermore, heat flux J(x) is uniform in x as expressed in the form for any x and x , which may work as another global relation for the local quantities. In order to connect the heat flux with the local thermodynamic quantities, we assume a heat conduction equation in which κ(T, ρ) is heat conductivity as a function of (T, ρ). We assume that there is no temperature gap at the boundaries, i.e., Finally, the conservation of particle number is written as These global relations, (40), (41) (42), (45), together with the equation of state (8), the heat conduction equation (43), and the boundary condition (44), may determine the profiles of the local temperature T (x) and the local density ρ(x) for systems without phase coexistence. As we will see in §4, there is a case where liquid and gas coexist in the container. For this special situation, the above global relations are not sufficient to determine the local states. This problem will be seriously studied in later sections.

Global thermodynamic quantities
Since thermodynamic quantities in heat conduction are not uniform in the space, they are basically described as local fields. Nevertheless, from the fact that the local states are governed by the global relations as explained in §2.4, we expect that some properties of the heat conduction systems are explained from a global point of view. Toward the characterization of global properties, we here define global thermodynamic quantities for heat conduction systems.
For extensive variables A originally defined for equilibrium states, we define the following global quantities as an extension to those for heat conduction states: where we have used the same notation A as that for equilibrium states. When we interpret A as a state function of heat conduction states, we explicitly write A(T 1 , T 2 , p, N ) which is in contrast to A(T, p, N ) for equilibrium states. For instance, global entropy and global Gibbs free energy are defined as Here, we investigate the reference state dependence of these global thermodynamic quantities. In equilibrium thermodynamics, entropy density and internal energy density are defined up to an additive constant which depends on the choice of their reference state. Thus, the entropy density s(x) and the internal energy density u(x) possess this property. Concretely, letŝ 0 andû 0 be the shift of the additive arbitrary constants of entropy and energy per one particle, respectively, for the change of the reference state. We express this transformation by The shift of other thermodynamic quantities are induced as Note that local thermodynamic relations are invariant under the transformation. Now, we consider the transformation of the global thermodynamic quantities, which are defined by the spatial integral of the local quantities. It is obvious that while for the global Gibbs free energy G, we have Theŝ 0 dependence of G is far from trivial. Here, we introduce the global temperatureT such that the transformation is written as Explicitly,T is given asT which means that global temperatureT corresponds to the kinetic temperature averaged over particles. The transformation (56) suggests the consistency among G, S andT . Indeed, we will show the global thermodynamic relations for these quantities. Similarly, we also define the global chemical potential as Sinceμ = G/N , the global chemical potential corresponds to the Gibbs free energy per one particle.
3 Global thermodynamics for single-phase systems in the linear response regime In this section, we restrict ourselves to a single phase where the heat conduction system is occupied by either liquid or gas at the constant pressure p ex . When the environmental parameters T 1 , T 2 , and p ex are fixed, the steady state is uniquely determined from the equation of state (8), the heat conduction equation (43), and the mass conservation law (45). That is, the values of (T (x), ρ(x), L z ) are determined by Since the local thermodynamic quantities obey equilibrium thermodynamics, the profile (T (x), ρ(x)) leads to any local thermodynamic quantities such as s(x) and µ(x). Then, the global quantities such asT ,μ, and S are calculated for the steady state. Now, we show that the global thermodynamic functions satisfy This means that the global free energy and entropy, which are functions of (T 1 , T 2 , p, N ), are expressed as the equilibrium free energy and the equilibrium entropy when we use the global temperatureT . This result leads to which indicates that the fundamental relation of thermodynamics can be extended to that for heat conduction states in the linear response regime. We call such a framework the global thermodynamics for heat conduction systems. Moreover, any global quantity A defined by (46) is connected to a corresponding equilibrium quantity A(T , p, N ) as 3.1 Proof of (62) and (63) When the environmental parameters (T 1 , T 2 , p ex ) are slightly changed to (T 1 +δT 1 , T 2 +δT 2 , p ex +δp ex ), the solution of the equations (59), (60), and (61) is modified slightly. We express the corresponding change as which leads to the change of any local thermodynamic quantity. For instance, the change of the local Gibbs free energy density is given by The change of the local quantities brings the change of the global quantities as where Here, since (37) holds for the local densities, the variation of each density satisfies a thermodynamic relation such as By substituting the local relations (35) and (73) into (72), we have whereŝ(x) = s(x)/ρ(x) and we have defined Note that n(x)dx is the particle number in [x, x+dx] and that n(x) is estimated as n(x) = N/L x +O(ε). η(x) is the deviation of the local temperature T (x) from the global temperatureT , and is estimated as η(x) = O(ε).
Since any local thermodynamic quantity a(x) is regarded as a(T (x), p ex ) = a(T + η(x), p ex ), we may estimate a(x) = a(T , p ex ) + O(ε). Then, the integrant in the right hand side of (74) is estimated asŝ where δn(x) = O(ε) for N and L x fixed. We find that holds from (57). Obviously, the conservation of the particle number leads to Thus, the integral in (74) is estimated as and then (74) becomes We emphasize that the definition ofT in (57) is essential for this relation.
Here, let us recall that the value of G is determined for a given (T 1 , T 2 , p ex , N ). Since (T , ε) may have the one-to-one correspondence to (T 1 , T 2 ), G is given as a function of (T , p ex , N, ε). The relation (80) implies that G is independent of ε in the linear response regime, and thus given as a function of (T , p ex , N ). Therefore, we conclude (62). Next, we show (63). From the relation (80), the left hand side of (63) is expressed as while the equilibrium fundamental relation results in Combining these two equalities, we obtain (63). From the definition (58) for the global chemical potentialμ and (62), we havẽ We then find whereŝ are entropy per one particle and specific volume of the heat conduction system.

Various global thermodynamic functions
From a local relation f (x) = g(x) − p, we have By using (62), we obtain We thus have the fundamental relation for F in the linear response regime. Next, we consider the global internal energy U and the global enthalpy H. From local thermodynamic relations u(x) = f (x) + T (x)s(x) and h(x) = g(x) + T (x)s(x), we write Remembering that η(x) = O(ε), the integral is estimated as where we have used (77). Thus, we obtain

Clausius equality
Since the thermodynamic relations are extended with keeping the same forms as those in equilibrium thermodynamics, we may define quasi-static heat d Q in an infinitely small quasi-static process (T 1 , T 2 , p ex ) → (T 1 + δT 1 , T 2 + δT 2 , p ex + δp ex ) as This d Q is free from the persistent heat flux in the steady state. For example, if the heat absorbed from the right heat bath is exactly the same as the heat released from the left heat bath at every moment, the net heat vanishes for the heat conduction. Such systems are thought to be adiabatic in the sense of global thermodynamics, and in such a quasi-static adiabatic process in the linear response regime. More generally, we expect that d Q corresponds to the excess heat proposed in [60], which can be measured in experiments. Then, constant-pressure heat capacity is defined as Applying the thermodynamic relation, we obtain Note that C p is defined as the response to the change of the global temperatureT . By changing the temperatures of heat baths,T may change in accordance with absorbing heat corresponding to C p δT .

Correspondence of global thermodynamic quantities to equilibrium quantities
In the previous subsections, we have obtained thermodynamic relations in the linear response regime and found that they are equivalent to the equilibrium ones. The key concept for extending thermodynamics is the global temperatureT and the global chemical potentialμ, by which the heat conduction system of (T ,μ) is mapped to the equilibrium system of (T ,μ). The connection between the two thermodynamic frameworks becomes clearer in the argument below by considering the estimation method of global thermodynamic quantities. In a single phase system, local temperature T (x) is a continuous monotonic function of x. Any extensive quantity A, which is defined as a spatial integral of local density a(x), can be transformed into an integral over temperature: where We expand ψ(T, P ) around T = T m in the form The second term is canceled when the integral (99) is performed. We thus have an estimation For A = V , we have a(x) = 1 for all x. By combining this with (103), we find We then have arrived at a universal estimation of the global quantity A as where A(T m , p) = V a(T m , p). The result (105) directly connects the global quantity in the heat conduction state to the corresponding equilibrium quantity at T m .
We may apply the similar method to the global temperature. The global temperature is also written asT in which By substituting it into (107), we obtainT Thus, in the linear response regime, the global temperatureT in the heat conduction system is equal to the mean temperature T m of the heat baths. The relation (105), together with (110), concludes (66). We also note that (109), which is rewritten as ρ = ρ(T , p) + O(ε 2 ) with ρ = N/V , corresponds to the global version of the equation of state We here remark that the relationT = T m in the linear response regime is a specific feature of systems in a rectangular container. We will show in §10.2 that the global temperature deviates from the mean temperature when the shape of the container is not a rectangle, whereas the relation (66) still holds.

Liquid-gas coexistence in heat conduction
From now on, we study liquid-gas coexistence in heat conduction. In §4.1, we review thermodynamics under equilibrium conditions while paying attention to the description of metastable states. In §4.2, we describe heat conduction states based on local equilibrium thermodynamics. We explain how the temperature profile and the density profile are determined for a given position of the liquid-gas interface. In §4.3, we present global thermodynamic quantities as a function of the interface position. In §4.4, we address the main problem of thermodynamics for the heat conduction systems with the liquid-gas coexistence.

Thermodynamics under equilibrium conditions
We first consider equilibrium thermodynamics. For a given (T, ρ), the pressure p is uniquely determined even when phase coexistence is observed. This equation of state is expressed as so as to emphasize that the pressure is the equilibrium value. See Fig. 4(a). A remarkable fact is that there is a plateau region [ρ G s (T ), ρ L s (T )] for a given temperature T below the critical temperature. From this equation of state, we find that the density ρ exhibits the discontinuous change at T = T c (p) when the temperature is changed with the pressure p fixed, as shown in Fig. 4(a). That is, the whole system is occupied by liquid when T < T c (p), while by gas when T > T c (p). In experiments, metastable states, which also may be called quasi-equilibrium states, are often observed. The most typical approach to the description of such meta-stable states is to employ the van der Waals equation of state: When T is fixed as that below the critical temperature, the equation of state contains the unstable states (∂ρ/∂p) T < 0. See Fig. 4(b). Following the standard procedure of the Maxwell construction, we obtain the equilibrium equation of state from (113). Here, we extract the two stable branches satisfying (∂ρ/∂p) T ≥ 0 from the curve defined by (113), each of which is connected to the equilibrium liquid phase or to the equilibrium gas phase. By solving (113) in ρ for each region, we obtain Hereafter, subscripts L and G are attached to quantities related to liquid and gas, respectively. When ρ G vdW (T, p) > ρ G s (T ), the spatially homogeneous phase of (T, p) corresponds to a meta-stable gas. Similarly, when ρ L vdW (T, p) < ρ L s (T ), the homogeneous phase of (T, p) corresponds to a meta-stable liquid. By plotting ρ as a function of T with p fixed as shown in Fig. 4(b), one finds that the meta-stable gas and the meta-stable liquid are super-cooled gas and super-heated liquid, respectively.
More generally, without using the van der Waals equation of state, we may obtain shows a jump at x = X. The density in the liquid region obeys ρ L (x; X) = ρ L (T (x; X), p), whereas that in the gas region by experimental measurements. See Fig. 4(c). Below we assume (116) and (117), whereas one may interpret them as (114) and (115). We then define densities for the liquid phase and the gas phase as and define quantities per one particle for the liquid phase and the gas phase aŝ

Heat conduction states
When T 1 < T c (p ex ) < T 2 , we may observe liquid-gas coexistence at the constant pressure p ex , which is in contrast with equilibrium cases. Concretely, the liquid occupies a lower temperature region, while the gas does a higher temperature region, as shown in Fig. 5(a). We assume that both the interface width and the temperature gap at the interface are negligible. See Fig. 5(b). Then, let X denote the position of the interface. A region 0 ≤ x < X is occupied by the liquid and the other region X < x ≤ L x is occupied by the gas. The density profile for a given X is then expressed as See Fig. 5(b). The temperature profile T (x; X) is determined by the heat conduction equation where J st (X) is the steady heat current which is constant in x due to the conservation of energy. The boundary conditions for the equation (121) are Concretely, T (x; X) is obtained as follows. We first express the interface temperature as θ(X): By setting we find that the local heat conductivity is estimated as Then, (121) leads to Solving the second equality of (126) in θ(X), we obtain where ζ and κ eff (ζ) denote the scaled position of the interface and an effective heat conductivity, which are defined by Note that which is obtained by substituting (127) into (126). By using θ(X) given by (127), the temperature profile is written as Substituting this profile into the equation of state (120), we have ρ(x; X) explicitly as exemplified in Fig. 5(c).

Global quantities as a function of X
The system with the phase coexistence is interpreted as a composite system of a liquid region and a gas region. For each subsystem, we may define global thermodynamic quantities introduced in §3.
Quantities defined in the liquid region and in the gas region are subscripted by L and G, respectively. In this section, we express global thermodynamic quantities as a function of X. First, the global temperatures in the liquid region and in the gas region are defined bỹ In each region, the temperature profile and the density profile are smooth. We describe each subsystem as a single phase system discussed in §3. The boundary temperatures of the liquid are T 1 and θ, and those of the gas are θ and T 2 . Thus, applying the trapezoidal rule explained in §3.4 to (131), the global temperature in each region is expressed as where we explicitly write the dependence on the interface position X. See Fig. 5(b). Next, the volume of the whole system V (X) is expressed as and thus the volume of the liquid region V L (X) and the volume of the gas region V G (X) are written as Then, the particle number in the liquid region, N L , is obtained as The particle number in the gas region, N G , is immediately given by By usingT L (X),T G (X), and N L (X), we can express all global thermodynamic quantities in the liquid region and the gas region. Explicitly, the extensive quantities are defined as Noting (66) in §3 with (27), we conclude that the global extensive quantities are written as Finally, we define the global temperature for the whole system by the formula (57), which is rewritten asT

Problem
We have shown that all thermodynamic quantities are determined for a given position X of the liquidgas interface. Since the liquid-gas interface in heat conduction may be at rest at constant pressure, the position X is uniquely determined for (T 1 , T 2 , p ex ). For equilibrium cases, the phase coexistence is observed in a container with V and T fixed. When we start with the van der Waals equation of state, we can determine the interface position of the two phases by the Maxwell construction, which is equivalent to the continuity of the chemical potential at the interface. Recalling this theory, one may impose the continuity of the local chemical potential µ(x) at the interface even in heat conduction [49]. Since µ L (T c (p)) = µ G (T c (p)), it leads to θ = T c (p ex ). That is, under this assumption, the interface temperature is equal to the equilibrium transition temperature.
Here, the important observation is that there is no justification of this condition for heat conduction systems. To our best knowledge, there are no experimental measurements on this issue, no numerical simulations of sufficiently large systems, and no reliable theory for supporting this condition. To be more precise, although the continuity of the chemical potential was reported in numerical simulations [67,68], the system size was too small to draw a definite conclusion. One may also recall from a viewpoint of non-equilibrium statistical mechanics that the local equilibrium distribution may be the leading contribution [37]. However, it should be noted that in the standard approach, fluctuations are described as those around the most probable profile, while the most probable profile itself is undetermined in the current problem. Therefore, the previous studies starting with the local equilibrium distribution do not result in the continuity of the chemical potential without imposing an additional assumption.
Furthermore, as another theoretical approach, one may study the stationary solution of the deterministic equation for the generalized Navier-Stokes equation with the interface thermodynamics [69,70]. We can estimate the discontinuous jump of the chemical potential at the interface as ε(w/L x )T c based on reasonable assumptions [39], where w represents the interface width. This result indicates the validity of the continuity of the chemical potential at the interface, because we may assume w/L x 1. However, as carefully argued in the paper [39], fluctuation effects should be taken into account for the generalized Navier-Stokes equation, so as to quantitatively describe the thermodynamic behavior. Now, based on these facts, we reconsider the continuity of the chemical potential at the liquid-gas interface within the framework of equilibrium thermodynamics. We then find that the continuity is equivalent to the thermodynamic variational principle. If the principle were applied to a local subsystem including the interface, the continuity of the chemical potential would be concluded. However, the constraint of a variational principle is generally applied to a whole system but not to the local subsystem, and therefore no variational principle is expected for the local subsystem near the interface. Thus we cannot justify the continuity of the chemical potential. This suggests an alternative approach to determine the interface position X. Namely, we formulate a variational principle for the whole heat conduction system as a natural extension of the equilibrium variational principle.

Variational principle for determining the liquid-gas interface position
Now, we formulate the variational principle for determining the interface position X and solve the variational problem. In §5.1, we start with the simplest example of the variational principle for equilibrium systems and naturally extend it to that for heat conduction systems. In §5.2 we derive the steady-state value by solving the minimization of the variational function. The solution is then reformulated as the temperature relation in §5.3. By using these results, in §5.4, we express the global quantities in terms of (T , p ex , Ξ).

Equilibrium systems
We start with an example of the variational principle for equilibrium systems. We consider a system at the constant pressure p ex in contact with a heat bath of the temperature T . Since the volume V of the system is not fixed, we determine V for a given (T, p ex ). When the equation of state is assumed, it is determined from the pressure balance equation We attempt to derive this condition from a variational principle. We take a variational function This corresponds to the Helmholtz free energy minimum principle for the composite system of the system and the volume reservoir whose Helmholtz free energy is p ex V + const. We then find that the minimization of G with respect to V for (T, p ex , N ) fixed leads to (141). The volume in the equilibrium state satisfies with which the Gibbs free energy is determined as It should be noted that this principle may be applied to the case that F (T, V, N ) is given by the van der Waals free energy: The variational principle for this case is equivalent to the Maxwell construction.

Extension to heat conduction systems
We study heat conduction systems at the constant pressure p ex . First, we focus on single-phase systems studied in §3. From the equivalence between a heat conduction system and the corresponding equilibrium system, the volume V of this system is determined by the minimization of whereT is the global temperature of the system. Since where L z = V/(L x L y ) depends on V andT should be fixed in the variation of V. Then, the variational principle (143) holds for single-phase systems. Next, we consider a system with a liquid-gas interface. In this case, V is not determined when we use the equation of state for each region, because the interface position X is not determined yet as described in the previous section. We then assume that the minimization of (148) is also valid for such cases, where Ξ = T 2 − T 1 is also fixed in the variation V . The last property is necessary to determine the value V uniquely. That is, we define a variational function and propose that the volume of the system V is determined as with This variational principle was first proposed in Ref. [66].

Remarks
Since the volume V is uniquely determined by the interface position X with (T , p ex , Ξ, N ) fixed, as explicitly shown in (133), the variational principle (150) with (151) is expressed as Similarly, since the particle number N L in the liquid region is uniquely determined by X, we also have another variational principle As the second remark, we provide a physical interpretation of (152). Let us consider a fluctuation of the interface position, X → X + δX. We assume that the motion of the interface is slowest and that the interface motion can be observed in a hypothetical system where (T 1 , T 2 ) is controlled such thatT and Ξ are fixed. See Fig. 6. It is then expected that fluctuations in this hypothetical system may be described by equilibrium statistical mechanics for the system with (T , p ex ). Thus, the most probable position is given by the left equation of (152) and the stability of the interface position is expressed as the right inequality in (152).
Moreover, we interpretT from an operational viewpoint. We virtually attach another heat bath of the temperature T 3 to the rigid top plate as shown in Fig. 7. Using the idealized assumption that the motion of the top plate is sufficiently slower than the other dynamical degrees of freedom, we control T 1 and T 2 with Ξ = T 2 − T 1 fixed such that the heat flux J 3 to the heat bath of T 3 is zero. Here, noting that the total kinetic energy in the bulk is given by 3N RT /2, J 3 is proportional toT − T 3 in the linear response regime. That is, J 3 = 0 is equivalent toT = T 3 , and therefore, the position of the top plate may obey the canonical distribution with temperature T 3 . In this manner, the variational principle withT is expected in this specific setting.

Steady state determined from the variational principle
We solve the variational equation (153). Hereafter, we consider the variation of N L with (T , p ex , Ξ) fixed. We abbreviate A(N L ;T , p ex , Ξ) as A(N L ) or A for the notational simplicity. We express the variational function G as where the free energy of the liquid F L (N L ), the free energy of the gas F G (N L ), and the volume of the system V (N L ) are given as functions of N L with (T , p ex , Ξ) fixed. Here, for any global quantity A in each region, such asT L/G and V L/G , we define We then have From the argument in §3, we obtain Since every global thermodynamic relation obtained in §3 holds in each region, we have Since the variation of the free energy in each region obeys the fundamental relation of thermodynamics for F , we have where By using (160), we rewrite the variation of G in (154) as where we used The global temperatureT L in the liquid region andT G in the gas region satisfy the relatioñ which are derived from the relationsT andT Furthermore, (163) leads to We have used N L + N R = N in the first line. Especially, the formulas (166) and (167), withT , Ξ, and N fixed, bring which simplifies (162) as Here, we confirmμ Let us define a specific temperature T s as By using (163), we rewrite (172) as Thus, the variation (169) is further simplified as This is equivalent to Since the functional form of µ G is different from that of µ L due to the crucial difference between the liquid and the gas, µ L (T s , p ex ) is not identically equal to µ G (T s , p ex ). The equality holds only when

Gas Liquid
Liquid-gas Thus we conclude that the variational principle (153) results in the unique steady value N L formulated by (178). Next, we consider the second derivative of G. By using (175) with (172), the second derivative of G is obtained as At the steady state satisfying (177), the above entropy difference is connected to the latent heatq per one particleq Therefore, we estimate the second derivative as Since Ξ > 0 andq(p ex ) > 0, we conclude that in the linear response regime.

Careful analysis of ε → 0
The solution (178) of the variational principle includes an undetermined term of O(ε). Nevertheless, (178) provides additional information to equilibrium behavior. To clarify the situation, we consider the limit ε → 0 with keeping the phase coexistence. This may be formalized by fixing a parameter r defined as whereas (178) indicates As 0 < N L < N , r satisfies In Fig. 8, we draw straight lines connecting from (T , Ξ) = (T c + rΞ, Ξ) to (T c , 0) for several values of r in the parameter region where the liquid-gas coexistence is observed. The convergence of these lines at (T c , 0) indicates that the equilibrium state (T c , 0) behaves as a singular point. By specifying the value of r, we identify the corresponding line which undoes the degeneration at equilibrium. N L is determined by considering the equilibrium limit ε → 0, while it is not uniquely determined at equilibrium.

Temperature relation
We obtained the steady state value of N L as (178) by solving the variational principle (153). Substituting this solution into (173) with (177), we have the global temperature for each region as We sum up the two relations in (186) and substitute (132) intoT L/G . We then obtaiñ which we call the temperature relation. This non-trivial relation, which results from the variational principle, is a simple condition that the steady states satisfy. Generally,T is not equal to T m , because the particle density of liquid is larger than that of gas, ρ L > ρ G . Thus, the relation (187) leads to θ = T c (p ex ), which indicates that the liquid-gas interface may not be observed at the position whose local temperature is equal to T c (p ex ). This means that metastable states stably appear near the liquidgas interface as schematically shown in Fig. 9. Now, suppose that (187) holds. Summing up the two equalities in (163), we obtain Here, it should be noted that we do not use the variational principle for the derivation of (163). By substituting the temperature relation (187) into this result, we derive (178) without the variational principle. In this sense, we may say that the temperature relation (187) is equivalent to the variational principle.

Global quantities as functions of (T , p ex , Ξ)
We determine the steady state values of several quantities. Let V L/G s (p ex ) denote the equilibrium saturated volume defined as with super-cooled gas Fig. 9 Super-cooled gas stabilized by the steady heat current.
The volume of the liquid (gas) region is given by Since N L is determined with an error of O(ε) in (178), (191) is the best estimate of V L/G in ε → 0. Substituting (178) into (191), we obtain which corresponds to the solution of the variational principle (150). Then, the steady position X of the liquid-gas interface, is determined as which is the solution of the variational principle (152). Then, the interface temperature is obtained from the elimination of T m from (187) and (127) as where ζ is given in (194). The deviation of θ from T c and the relation between X and θ will be demonstrated in §6 for explicit examples. In this section, we express any quantity in the steady state as a function of (T 1 , T 2 , p ex ). In §6.1, we derive (T , ζ) as a function of (T 1 , T 2 , p ex ). Since we already have the relationsT lead to the expression Particularly, we consider the interface temperature θ and the jump of the chemical potential ∆µ at the interface as a function of (T 1 , T 2 , p ex ) in §6.2 and §6.3, respectively. In §6.4, we discuss how the super-cooled gas appears near the interface quantitatively. The last subsection §6.5 is devoted to the demonstration of examples. In what follows, we use the notations for simplicity.
6.2 Interface temperature as a function of (T 1 , T 2 , p ex ) Substituting (203) into (201), we obtain another expression ofT : Then, from the temperature relation (187), we find that the interface temperature is written as This formula clarifies that the interface temperature θ(ζ) generally deviates from T c (p ex ) when the temperature gradient is imposed. Since the particle density satisfies ρ L > ρ G , i.e., u > 1, and the heat conductivity is expected to be κ L = κ G , the interface temperature θ is deviated from the equilibrium transition temperature T c (p ex ) in the order of ε, which is not negligible.
Here, we show the simplest expression of the interface temperature θ: whereφ = V /N . The derivation of this formula is shown below. We first recall (126) and (127). From these, Combining these two relations, we have another form of T m as with an error of O(ε 2 ). Second, the result of the variational principle (178) is written as As N L = ζV ρ L (T L , p ex ) and N G = (1−ζ)V ρ G (T G , p ex ) with ρ L (T L , p ex ) = ρ L c +O(ε) and ρ G (T G , P ) = ρ G c + O(ε), the relation (212) is further expressed as with an error of O(ε 2 ). By subtracting (211) from (213) and using the temperature relation (187), we have two forms of the interface temperature: By substituting (214) and (215) into the first term and the second term of the right-hand side of the trivial identity we obtain (208). 6.3 ∆µ as a function of (T 1 , T 2 , p ex ) Next, we discuss the chemical potential jump at the interface. Note that which gives the definition of the transition temperature T c . Thus, θ = T c means the imbalance of the chemical potential at the interface, which is quantitatively expressed as Since the entropy difference is connected to the latent heatq per one particle aŝ the jump of the chemical potential is proportional toq. Now, substituting (207) into (219), the amount of the jump is estimated as Another expression of ∆µ is also obtained by substituting (208) into (219). The result is

Super-cooled gas in the liquid-gas coexistence
When T 1 < T c (p ex ) < T 2 , we find the position X c at which the local temperature satisfies T (X c ; T 1 , T 2 , p ex ) = T c (p ex ). Suppose θ < T c (p ex ). Then, the position X of the interface satisfies We observe the liquid in the region 0 < x < X and the gas in X < x < L x . Then the local temperature of the gas in X < x < X c is less than T c (p ex ). This means that a super-cooled gas is observed in X < x < X c , which is not stable in equilibrium.   We plot (204) and (207) in Fig. 10 for three sets of (u, v). It is clearly seen that the interface temperature θ deviates from the equilibrium transition temperature T c (p ex ). The deviation of the global temperatureT from T m shows the same figure with Fig. 10(b) due to the temperature relation (187). The jump of the chemical potential at the interface also exhibits the same dependence on ζ = X/L x since it is proportional to θ − T c (p ex ) as shown in (221). In the presented three sets of (u, v), a super-cooled gas region appears near the right side of the interface, as shown schematically in Fig. 9.
According to Fig. 10(a), the stable position of the interface is shifted continuously as we increase the temperature of both heat baths. This feature is different from the numerical solution of the variational equation reported in Ref. [66], where the interface discontinuously appears. A more accurate numerical calculation was necessary for the cases with high ρ L c /ρ G c .

Examples
In this subsection, we illustrate some examples of the liquid-gas coexistence in heat conduction quantitatively. First, we consider the van der Waals equation of state of a = 0.365 Pa· m 6 /mol 2 and b = 4.28 × 10 −5 m 3 /mol for CO 2 at a high pressure p ex = 4 × 10 6 Pa [71]. The gas constant is R = 8.31 J/K·mol. Substituting these parameters into (224), we obtain the transition temperature as T c = 262.7 K, and the mol density at T c as ρ L c = 1.38 × 10 4 mol/m 3 and ρ G c = 2.70 × 10 3 mol/m 3 , which results in u = 5.1. From the database [72], we set the heat conductivity of the liquid as κ L = 0.1 W/m· K and that of the gas as κ G = 0.02 W/m· K, which yields v = 5. The result is close to the result (u, v) = (5.5) which was already shown in Fig. 10. The temperature profile becomes linear both in the liquid and gas regions in this example due to the constant heat conductivity. Then the volume fraction of the super-cooled gas is simply obtained as where V scG is the volume of the super-cooled gas, i.e., V scG = (X −X c )/L x . The interface temperature θ is given by (207) and T 2 results from (203). In Figs. 11, we show the interface temperature θ and the volume fraction V scG /V as a function of the position of the interface when Ξ = 10 K. We notice that the interface temperature θ may deviate about 3.3 K from T c , and that the volume fraction of the super-cooled gas may exceed 1/3. Second, more importantly, we do not need the equation of state, but have only to know the value of u, v, T c (p ex ), when we obtain the phase diagram for a given material. For example, consider pure water at 1.013 × 10 5 Pa, where T c (p ex ) = 373.1 K. From the database [72], we have u = 1604 and v = 27.06. From these, we can predict the interface temperature θ and the volume fraction V scG /V as a function of the position of the interface as shown in Fig. 12. As an example, when T 1 = 368.0 and T 2 = 378.0, we obtain θ = 368.3 and ζ = 0.4086.

Global thermodynamics for systems with a liquid-gas interface
Hereafter, all the quantities are defined in the steady state. When a global quantity A is determined for (T , p, N, Ξ), we express this relation as Recall that A(T , p, N, Ξ) = A(T , p, N ) + O(ε 2 ) for single-phase systems, where A(T , P, N ) is the equilibrium function, as discussed in §3. For the systems with a liquid-gas interface, Ξ dependence appears even in the linear response regime. In this section, we extend the global thermodynamics introduced in §3 so as to describe the liquid-gas coexistence determined in the previous sections. In §7.1, we derive formulas as preliminaries for later sections. In §7.2, we formally derive a fundamental relation of thermodynamics for systems with a liquid-gas interface. However, since the derivative of G is not defined at (T c (p), p), the formal expression is not properly defined. Then, in §7.3, we perform a careful analysis of the limit ε → 0. In §7.4, we show the form of the entropy and the volume in the appropriate limit. In §7.5, we present formulas for constant pressure heat capacity and compressibility.

Preliminaries
We first study global extensive quantities It is obvious that they possess the additivity where A L and A G are defined as global thermodynamic quantities in the liquid region and the gas region, respectively. By usingT L andT G , we write Now, we consider an infinitely small change (T , p, N, Ξ) → (T + δT , p + δp, N + δN, Ξ + δΞ), which may be realized by controlling the experimental condition (T 1 , T 2 , p, N ). For any function φ(T , P, N, Ξ) of steady states, δφ is defined as We then have We sum up the above two relations (233) and (234) and substitute (166) and (167) into them. Then, the sum becomes The fourth line is transformed intô By summarizing the terms proportional to δ(N L /N ) in (235) and (236), the coefficient of N δ(N L /N ) becomesâ Here, we note that the sum of the two relations forT L andT G in (163) results iñ where the transformation from the first line to the second line is found by notingT L = (T 1 +θ)/2+O(ε 2 ) andT G = (T 2 + θ)/2 + O(ε 2 ) from (110). Now, using the temperature relation (187), we find that the right-hand side of (238) is T c (p). We thus obtain Here, the last term should be expressed as a linear combination of δT , δp, δN and δΞ. Indeed. from (178) which gives N L as a function of (T , p, Ξ), we derive We leave the formula (239) and (240) as they are for later convenience.

Formal derivation of fundamental relation
We set in the formulas in the previous subsection. Then, we recall the thermodynamic relations whereφ is the specific volume. We also define whereq is the latent heat defined in (180). By using them, we rewrite (239) as Here,μ is the global chemical potential defined by (58). It should be noted that whereμ L/G is the global chemical potential in each region such thatμ L = µ L (T L , p) + O(ε 2 ) and µ G = µ G (T G , p) + O(ε 2 ). Now, noting the relations and (240), we obtain δG = −(S + O(ε 2 ))δT + (V + O(ε 2 ))δp + (μ + O(ε 2 ))δN − (Ψ + O(ε)) δΞ By formally considering the infinitely small change δ, we obtain which corresponds to the fundamental relation of thermodynamics.
At first sight, one may be afraid that the relation (248) does not provide a non-equilibrium extension because the relation involves the error of O(ε). However, it includes non-trivial information. The point is that some equilibrium thermodynamic quantities are singular at T = T c (p). For example, S and V are not uniquely defined at (T, p) = (T c (p), p) in equilibrium. Nevertheless, (248) provides a definition of S at T = T c (p) as the limit ε → 0 of the derivative of G with respect toT while fixing (p, Ξ, N ). We study (248) more carefully by going back to (247).

Careful analysis of ε → 0
We consider the limit ε → 0 with keeping the phase coexistence. The parameter r defined in (183) identifies the line terminating to the equilibrium state as exemplified in Fig. 8, and undoes the degeneration of the equilibrium state. Since (183) is rewritten as Then, let us define the equilibrium limit of entropy and volume by S eq (r) ≡ lim ε→0+ S(r, p, N, Ξ), where we explicitly write only the r dependence for S eq and V eq . Substituting (249) into (247), we obtain where we have used (184). This relation leads to Since we find from (249) that the left hand side is equal to we obtain lim ε→0 ∂G ∂T p,N,Ξ = −S eq (r).

Explicit forms of S eq and V eq
By using (184), we express any extensive quantity A as This leads to where we have used the formula (180) for the latent heatq. Similarly, we obtain where we have used the Clausius-Clapeyron relation Note that V eq (r) in (265) is consistent with V (T , p, Ξ) in (192) at the limit ε → 0.

Heat capacity and compressibility
Let H(T , p, N, Ξ) be the enthalpy for the heat conduction system with the phase coexistence. Since H is given by the spatial integral of the local enthalpy, it satisfies This leads to in an infinitely small quasi-static process at constant pressure. We interpret this δH as the quasi-static heat d Q in this process. We thus define heat capacity at constant pressure as where the specific heat for the liquid (gas) is defined aŝ By using (178) and (269), we rewrite (270) as with the latent heatq =ĥ G c −ĥ L c . Note that the third term mainly contributes to C p in (272) because it diverges as Ξ → 0. This is consistent with the singularity at T c (p) for equilibrium cases, where C p is given by the derivative of a discontinuous function (enthalpy) at T = T c (p). The first and the second terms correspond to the heat capacity of the liquid region and the gas region, respectively. Thus, the third term may be interpreted as the interface contribution. The expression (272) also clarifies the violation of additivity for C p .
Similarly to the heat capacity, we can also study the additivity of other response functions. As one example, we show the singular nature of the compressibility from the viewpoint of the violation of the additivity. The compressibility in heat conduction is defined as By taking A = V in the formula (239), we obtain Similarly to C p , the main contribution in Ξ → 0 is the non-additive term, which diverges in the equilibrium limit.
super-cooled gas

Fundamental relation
Since F and G are given as the spatial integral of the local free energies, the uniformity of the pressure leads to Substituting this into (247), we have  (T , p, N ). Therefore, this term does not exhibit the singularity in the limit ε → 0. Indeed, the second line is estimated as where we have used ∂p Thus, the last term in (276) remains to be negligible as O(ε 2 ) term, and we conclude It should be noted that there are the errors of O(ε 2 ) in contrast to (248).

Free energy for the coexistence phase
In the equilibrium coexistence phase, as shown in Fig. 14, the saturated volumes, V L s (T ) and V G s (T ), and the saturated pressure p s (T ) are defined for a given temperature T . Even for the coexistence phase in heat conduction, we define the saturated volumes, V L s (T , Ξ) and V G s (T , Ξ), at which the liquid and the gas start to coexist. Note that the pressure is not kept constant in V L s (T , Ξ) < V < V G s (T , Ξ), differently from the equilibrium case.
We first express V L s (T , Ξ), V G s (T , Ξ) and p(T , V, Ξ) in terms of V L s (T ), V G s (T ) and p s (T ). Since the steady state satisfies (178), we have Solving this in p, we obtain For equilibrium cases, the relations N L + N G = N and By using these expressions, we rewrite (281) as where V m (T ) = (V L s (T ) + V G s (T ))/2. Figure 14 shows a p-V curve described by (283), in which p(T , V, Ξ) is linearly decreasing with V in the range V L s (T , Ξ) < V < V G s (T , Ξ). Configurations in the phase coexistence are exemplified in Fig. 15. Then, the saturated volumes under heat conduction are expressed as where V = V L (T, p) and V = V G (T, p) are the equilibrium equation of state for the liquid and the gas, respectively. Next, we study the free energy F (T , V, Ξ). We first consider the free energy difference Since the saturated state is in a single phase, we have We then transform where we have used the Taylor expansion in Ξ to obtain the second line, and substituted the equilibrium relation to obtain the last line. Second, from (279), we have          Performing the integration with (283), we obtain

K X 6 b h W e C P E E I W m U j g g v 4 Y = " > A A A C Y 3 i c h V G 7 S g N B F D 1 Z 3 / E V o 4 U g Q j A o V u G u j Y 8 q K I i l M U Y F F d l d J 3 F w X + x u A h r 8 A W 0 V C y s F E f E z b P w B i / y A I J Y R b C y 8 2 S y I i n q H m T l z 5 p 4 7 Z 2 Z 0 1 5 R + Q F S L K S 2 t b e 0 d n V 3 x 7 p 7 e v v 7 E Q H L N d 8 q e I Q q G Y z r e h q 7 5 w p S 2 K A Q y M M W G 6 w n N 0 k 2 x r u 8 v N P b X K 8 L z p W O v B g e u 2 L a 0 k i 2 L 0 t A C p n K L O 4 k 0 Z S i M 1 E + g R i C N K J a d x A 2 2 s A s H B s q w I G A j Y G x C g 8 9 t E y o I L n P b q D L n M Z L h v s A R 4 q w t c 5 b g D I 3 Z f R 5 L v N q M W J v X j Z p + q D b 4 F J O 7 x 8 o U x u m R b q l O D 3 R H z / T + a 6 1 q W K P h 5 Y B n v a k V 7 k 7 / 8 X D + 7 V + V x X O A v U / V n 5 4 D F D E T e p X s 3 Q 2 Z x i 2 M p r 5 y e F 7 P z 6 2 M V y f o i l 7 Y / y X V 6 J 5 v Y F d e j e u c W L l A n D 9 A / f 7 c P 0 F h K j O b U X O U z s 5 H P 9 G J E Y x h k p 9 7 G l k s Y R k F P l b g B K c 4 i z 0 p P U p S G W q m K r F I M 4 g v o Y x + A A E Y i Z Y = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " 9 4 b q j 4 K X 6 b h W e C P E E I W m U j g g v 4 Y = " > A A A C Y 3 i c h V G 7 S g N B F D 1 Z 3 / E V o 4 U g Q j A o V u G u j Y 8 q K I i l M U Y F F d l d J 3 F w X + x u A h r 8 A W 0 V C y s F E f E z b P w B i / y A I J Y R b C y 8 2 S y I i n q H m T l z 5 p 4 7 Z 2 Z 0 1 5 R + Q F S L K S 2 t b e 0 d n V 3 x 7 p 7 e v v 7 E Q H L N d 8 q e I Q q G Y z r e h q 7 5 w p S 2 K A Q y M M W G 6 w n N 0 k 2 x r u 8 v N P b X K 8 L z p W O v B g e u 2 L a 0 k i 2 L 0 t A C p n K L O 4 k 0 Z S i M 1 E + g R i C N K J a d x A 2 2 s A s H B s q w I G A j Y G x C g 8 9 t E y o I L n P b q D L n M Z L h v s A R 4 q w t c 5 b g D I 3 Z f R 5 L v N q M W J v X j Z p + q D b 4 F J O 7 x 8 o U x u m R b q l O D 3 R H z / T + a 6 1 q W K P h 5 Y B n v a k V 7 k 7 / 8 X D + 7 V + V x X O A v U / V n 5 4 D F D E T e p X s 3 Q 2 Z x i 2 M p r 5 y e F 7 P z 6 2 M V y f o i l 7 Y / y X V 6 J 5 v Y F d e j e u c W L l A n D 9 A / f 7 c P 0 F h K j O b U X O U z s 5 H P 9 G J E Y x h k p 9 7 G l k s Y R k F P l b g B K c 4 i z 0 p P U p S G W q m K r F I M 4 g v o Y x + A A E Y i Z Y = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " 9 4 b q j 4 K X 6 b h W e C P E E I W m U j g g v 4 Y = " > A A A C Y 3 i c h V G 7 S g N B F D 1 Z 3 / E V o 4 U g Q j A o V u G u j Y 8 q K I i l M U Y F F d l d J 3 F w X + x u A h r 8 A W 0 V C y s F E f E z b P w B i / y A I J Y R b C y 8 2 S y I i n q H m T l z 5 p 4 7 Z 2 Z 0 1 5 R + Q F S L K S 2 t b e 0 d n V 3 x 7 p 7 e v v 7 E Q H L N d 8 q e I Q q G Y z r e h q 7 5 w p S 2 K A Q y M M W G 6 w n N 0 k 2 x r u 8 v N P b X K 8 L z p W O v B g e u 2 L a 0 k i 2 L 0 t A C p n K L O 4 k 0 Z S i M 1 E + g R i C N K J a d x A 2 2 s A s H B s q w I G A j Y G x C g 8 9 t E y o I L n P b q D L n M Z L h v s A R 4 q w t c 5 b g D I 3 Z f R 5 L v N q M W J v X j Z p + q D b 4 F J O 7 x 8 o U x u m R b q l O D 3 R H z / T + a 6 1 q W K P h 5 Y B n v a k V 7 k 7 / 8 X D + 7 V + V x X O A v U / V n 5 4 D F D E T e p X s 3 Q 2 Z x i 2 M p r 5 y e F 7 P z 6 2 M V y f o i l 7 Y / y X V 6 J 5 v Y F d e j e u c W L l A n D 9 A / f 7 c P 0 F h K j O b U X O U z s 5 H P 9 G J E Y x h k p 9 7 G l k s Y R k F P l b g B K c 4 i z 0 p P U p S G W q m K r F I M 4 g v o Y x + A A E Y i Z Y = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " 9 4 b q j 4 K X 6 b h W e C P E E I W m U j g g v 4 Y = " > A A A C Y 3 i c h V G 7 S g N B F D 1 Z 3 / E V o 4 U g Q j A o V u G u j Y 8 q K I i l M U Y F F d l d J 3 F w X + x u A h r 8 A W 0 V C y s F E f E z b P w B i / y
Here, by using (282), we can confirm Combining it with the Clausius-Clapeyron relation we derive where Ψ was introduced in (243). Summarizing these, we schematically draw the graph F in Fig. 16. For equilibrium cases, the free energy of the liquid-gas coexistence phase is expressed as a common tangent at V L s (T ) and V G s (T ) as shown in Fig. 16(a). For heat conduction cases, the free energy of the coexistence phase is expressed as a common tangential quadratic curve at V L s (T , Ξ) and V G s (T , Ξ) as shown in Fig 16(b), and therefore F keeps the convexity on V .
Since the convexity of F is concluded, the Gibbs free energy G is expressed as the Legendre transform of F : The graph F in Fig. 16(b) yields the graph G in Fig. 17(b), while it contains an error of O(ε 2 ) as well as F formulated explicitly in (290). The liquid-gas coexistence states are expressed as a common tangential quadratic curve to the graph G at p s (T + Ξ/2) and p s (T − Ξ/2). We have studied the dependence of the steady states on V for (T , Ξ) fixed. However, when we fix values of T 1 and T 2 , the change of V yields the change ofT . Here we propose a protocol to control T 1 and T 2 , by whichT is fixed when changing V . In other words, we derive T 1 and T 2 as functions of (T , V, N, Ξ).
We also give the interface temperature θ and its deviation from T c as functions of (T , V, N, Ξ).
In (203), T m is expressed in terms of T c (p), Ξ, u, v, and ζ, where the value of u and v are determined for each p. For the coexistence phase in heat conduction, the pressure is determined as (283) and can be replaced by p s (T ) for the argument of u and v. We thus represent u = u(p s (T )) and v = v(p s (T )), and simply write u(T ) and v(T ). Then, we rewrite (203) as with Furthermore, from and These expressions of u, v and ζ lead to with which (296) is rewritten as Here, let us recall that T c (p) =T − Ξr defined in (183) and that as used in (283). Substituting this formula into (303) and noting we finally arrive at Since T m = (T 1 + T 2 )/2 and Ξ = T 2 − T 1 , we obtain These formulas inform how to operate T 1 and T 2 in order to fix the global temperatureT when changing V and Ξ.
Substituting the temperature relation (187) into (303), we obtain the interface temperature as Similarly, (306) and (187) lead to the deviation from the transition temperature, 9 Other assumptions leading to the results of the variational principle To this point, it has been assumed that the steady states are determined by the variational principle proposed in §5. In this section, we introduce other assumptions from which the result of the variational principle is obtained. First, in §9.1, we formulate the variational principle for the constant volume systems and show that its solution is equivalent to the previous one (178). In §9.2, we start with the fundamental relation of thermodynamics without assuming the variational principle. We then derive the result of the variational principle. In §9.3, we notice a singularity relation which is a simple assumption for the singularity in the limit ε → 0. We confirm that this is equivalent to the result of the variational principle. Finally, in §9.4 we argue some scaling behavior of the system which also leads to the result of the variational principle.

Variational principle for constant volume systems
In this subsection, we study constant volume systems, where the steady state is characterized by (T , V, Ξ, N ). As shown in Fig. 13, this steady state is equivalent to the steady state of the system at the constant pressure p ex when the value of p ex is chosen as the pressure p(T , V, N, Ξ). We assume the mechanical balance everywhere. That is, for any x and x in [0, L x ], where p(x) = p(T (x), ρ(x)). This condition determines V L and V G = V −V L when N L and N G = N − N L are determined. The problem is then to determine N L . After briefly reviewing the corresponding variational principle for equilibrium cases, we propose the variational principle for heat conduction systems and determine the value of N L .

Equilibrium systems
We define the variational function as in which F L/G (T, V, N ) is the Helmholtz free energy of the liquid (gas). V L and V G are determined by and V L + V G = V . Then, the variational principle for determining N L is which is equivalent to the equality of the chemical potential.

Heat conduction systems
Let us consider a variational principle for the heat conduction systems characterized by (T , V, N, Ξ). We assume that the steady state with a liquid-gas interface is determined by the variational function where we have used a similar formula as (228) with (229) and (230) to obtain the second line. As we formulated in §4, thermodynamic quantities are expressed as a function of N L such thatT L/G = T L/G (N L ;T , V, N, Ξ). The variation of F is then written as Remembering that the relation (160) holds in each region, we have Substituting δN G = −δN L , δV L + δV G = 0 and (168) into (317), we rewrite (316) as where T s is given in (172). Then, the variational principle determines N L as that satisfying This is equivalent to (176) so that (320) leads to (178). Therefore, the variational principle for the constant volume system is equivalent to that for the system at constant pressure.

Gas Liquid
Liquid-gas

Thermodynamic relation
In this subsection, we do not assume any variational principles, but assume the fundamental relation of thermodynamics for the liquid-gas coexistence phase in the linear response regime.
First of all, it should be noted that (235), (236), (237) and (238) hold regardless of the variational principle. Then, we have instead of (239). Letting A = G andâ =μ, and using G = F + pV , we obtain where Since N δ(N L /N ) = O(ε 0 ), the assumption (321) leads to This means the temperature relation (187). As shown in §5.3, (187) leads to (178) so that the thermodynamic relation is equivalent to the variational principle.

Singularity relation
We here consider a protocol to shift the liquid-gas interface from X = L x to X = 0 with fixing Ξ and p. This protocol is obtained by varying (T 1 , T 2 ) as In this protocol, (T , N L ) is changed as For sufficiently small Ξ, it is natural to assume that This differential equation is written as which we call a singularity relation. Solving the differential equation (328) with the boundary conditions andT we obtainT which is the same form as the result of the variational principle (178). The relation (329) indicates which corresponds to an expression of the singularity associated with the first-order transition for equilibrium cases. This singularity is consistent with the discontinuous change of N L at T = T c (p). It should be noted that (333) is connected to the singularity of constant pressure heat capacity C p and compressibility α T as shown in §7.5.

Scaling relation
We characterize the coexistence phase by (p ex , N L , Ξ, N ). As examples, we writẽ From the homogeneity in the direction perpendicular to x,T L andT G are invariant for (p ex , N L , Ξ, N ) → (p ex , λN L , Ξ, λN ). We thus writeT where By notingT G −T L = Ξ/2, we expressT L andT G as for small Ξ ≥ 0. Now, we attempt to extend the gas region while keeping the thermodynamic state in the liquid region. See Fig. 19. First, we increase Ξ as Ξ → λΞ withT L fixed, where we take λ ε −1 such that the system is still in the linear response regime. We estimate the λ dependence of quantities in this asymptotic regime. First, the heat flux is proportional to λ as a common value to the liquid and the gas region. SinceT L is fixed, it is plausible that the temperature difference in the liquid region remains to be O(λ 0 ), and then this leads that the horizontal length of the liquid region is proportional to 1/λ. On the other hand, in the gas region, the temperature difference is proportional to λ and the horizontal length is L x in the leading order estimate. Thus, the total volume saturates at the finite value with which the whole system is occupied by the gas with the temperature less than T c (p ex ) + λΞ. From these, we find that the volume of the liquid region is proportional to 1/λ. Although its proportional coefficient is not determined in the asymptotic region, we assume the scaling relation that the volume of the liquid is V L /λ as it is consistent with the case λ = 1. Since p ex andT L is fixed in this operation, the density of the liquid is also fixed. Thus, the particle number of the liquid becomes N L /λ. Next, we increase N as N → λN , Then, we have N L /λ → N L and V L /λ → V L . That is, for a series of processes (p ex , N L , Ξ, N ) → (p ex , N L , λΞ, λN ), we haveT L (p ex , N L , λΞ, λN ) =T L (p ex , N L , Ξ, N ). Fig. 19 Scaling to extend the gas region with keeping the state in the liquid region and the horizontal length Lx. In order to fixT L , the temperature of the left heat bath is changed from T 1 to T 1 .
Since this condition cannot be concluded from the local equilibrium thermodynamics, we impose (342) as a requirement for the steady state.
10 Generalization for single phase systems In previous sections, global thermodynamics has been developed for systems in a rectangular container as shown in Fig. 2. In this section, we consider an arbitrarily shaped container D in contact with two heat baths of T 1 and T 2 . Here, we assume that the response of a local quantity is not singular for the change of the system parameters (T 1 , T 2 , V ) → (T 1 + δT 1 , T 2 + δT 2 , V + δV ), where V is the volume of the container. We also assume that local steady states inside the container are well described by the local equilibrium thermodynamics of local quantities T (r), u(r), s(r), and so on. In §10.1, we define global thermodynamic quantities as they are consistent with the definitions given in §3. We construct global thermodynamics in the linear response regime in §10.2. Then, we extend the thermodynamic formulas beyond the linear response regime in §10. 3. In what follows, we restrict to single-phase systems without liquid-gas interfaces.

Global thermodynamic quantities in an arbitrarily shaped container
We define all global thermodynamic quantities similarly to the case of rectangular containers. Explicitly, the global temperature and the global chemical potential are defined as andμ More generally, global extensive quantities A are defined by the spatial integration of local thermodynamic quantities a(r) per unit volume orâ(r) per one particle as It should be noted that all global quantities transform consistently for the change of the reference state on entropy and internal energy as discussed in §2.5.
Due to the mechanical balance, the local pressure p(r) is homogeneous. That is, We here note that the pressure p(r) does not satisfy the local equilibrium equation of state beyond the linear response regime. Concretely, the contribution out of the local equilibrium pressure was explicitly calculated in the kinetic regime by analyzing the Boltzmann equation [73], and also the long-range correlation yields the additional pressure depending on the system size [74,75,76]. Thus, we express where p LE (r) is the pressure determined from the local equilibrium equation of state such that .
This may deviate from the global pressure as

Equivalence of non-equilibrium global quantities with equilibrium quantities
We show that the results in §3 hold in the container with an arbitrary shape in the linear response regime. As we did in §3.1, we introduce a function η(r) = O(ε) as η(r) ≡ T (r) −T .
This indicates that all global thermodynamic quantities are equivalent to those in equilibrium by adopting the global temperatureT in (352), and that global thermodynamics for heat conduction systems are generally mapped to equilibrium thermodynamics, regardless of the shape of the container. It should be noted that we cannot apply the trapezoidal rule argued in §3.4 to such general configurations. Thus, in contrast to the case of rectangular containers, we havẽ even in the linear response regime.

Fundamental relation beyond the linear response regime
Since global thermodynamics relies on local equilibrium thermodynamics, its validity depends on the extent how the local states are close to the local equilibrium. For instance, the local equilibrium thermodynamic relations may be valid up to O(ε) and a non-equilibrium driving may bring a slight deviation of O(ε 2 ) from the local equilibrium. We here expect that the deviation from the local equilibrium is small even beyond the linear response regime, and attempt to extend the global thermodynamics to the second order of ε. Indeed, in §10.3.1 and §10.3.2, we prove where δν is either δT , δV or δN . Ψ is the conjugate variable to Ξ, which is identified as This is a quantity of O(ε) with constant pressure specific heatĉ p . α is a geometrical factor independent ofT and p as examined in §10.3.3. That is, we can extend the thermodynamic framework with an error of O(ε 3 ) beyond the linear response regime.
where the first term in the right hand side of the first line turns out to be zero from the definition of T . Thus, (385) is rewritten as which corresponds to (365).

Concluding Remarks
In this long paper, we have developed a thermodynamic framework for heat conduction states, which we call global thermodynamics. The key concept of the framework is the global temperature defined by (57). We describe spatially inhomogeneous systems by a set of global thermodynamic variables. Although the method is rather different from a standard continuum description, global thermodynamics is built on local equilibrium thermodynamics and the transportation equation. Thus, it is not contradictory with established theories. Nevertheless, our framework provides a new prediction on the thermodynamic properties of a liquid-gas interface in heat conduction. As we explained in §4.4, the rule of connecting two phases at the interface is out of the standard hydrodynamics with local equilibrium thermodynamics. Even for this case, global thermodynamics can determine the way how to connect the two phases by the variational principle that is a natural extension of the Maxwell construction for equilibrium cases. Although global thermodynamics is self-consistent, self-contained, and natural, it does not ensure the validity of the framework. Similarly to other universal theories such as thermodynamics, hydrodynamics, and statistical mechanics, global thermodynamics involves a fundamental assumption on the determination of the steady state. In other words, if our quantitative prediction on the phase coexistence is denied by experiments, global thermodynamics could not be valid. In this case, we will attempt to understand the experimental results carefully so that we have a correct description of phenomena out of equilibrium. This is of course quite stimulating, and we sincerely wait for experiments. Numerical simulations of microscopic dynamics such as molecular dynamics also provide useful materials for further studies. To the present, no deviation of the interface temperature from the transition temperature was observed [67,68]. We conjecture that this is due to the insufficiency of the separation of scales, because the deviation should be observed in systems with enough separation of scales, as discussed in Ref. [39].
Putting aside the ultimate validity of global thermodynamics, we emphasize that the framework itself is quite fruitful. This paper presents only the minimum backbone, but exhibits non-trivial simple relations among several quantities. These somewhat miraculous statements may suggest a deeper structure behind global thermodynamics. Moreover, importantly, the story does not end. We did not argue the Clausius equality for the cases with a liquid-gas interface, while we confirmed it for single-phase systems in §3.3. The equality connects "heat" to the entropy change, and the "heat" should be defined as the renormalized one, called "excess heat", after subtracting the persistent heat as we reviewed in Introduction. It is quite natural to study the liquid-gas coexistence by considering an extended form of the Clausius relation. Since we find that this provides another source of further development of the framework, we have decided to argue it in a separate paper.
From a theoretical viewpoint in a broad context, we expect that we may derive global thermodynamics based on a more microscopic description. Since we know the basic framework of non-equilibrium statistical mechanics, e.g. the linear response theory, we obtain a linear response formula of the temperature of the interface. It is quite formal, and it is not evident that the result satisfies the temperature relation (187). Related to this problem, recently, we have calculated the temperature of the interface for a stochastic order parameter equation describing the phase coexistence in heat conduction [39]. The result shows that the interface temperature deviates from the equilibrium transition temperature. It should be noted that global thermodynamics can be formulated for systems that exhibits an order-disorder transition in heat conduction. By finding a connection between the two theories, we will understand a microscopic mechanism of the interesting macroscopic phenomenon.
Our theory may be naturally generalized to that for heat conduction systems where a liquid flows in and a gas flows out. When a steady state is observed, the system exhibits perpetual evaporation at the interface. The latent heat is persistently generated at the interface and as a result additional heat flux occurs from the interface to the liquid region and gas region. This phenomenon was experimentally studied for pure water [77], and a macroscopic temperature gap as much as about 8K was found at the liquid-gas interface. As a related phenomenon, the temperature gradient in the inverse direction to the imposed temperature difference was observed for the liquid-gas coexistence of a heat conducting Helium far from equilibrium [78]. These phenomena may be closely related to our prediction on the phase coexistence in heat conduction. We will attempt to explain such phenomena in global thermodynamics. Furthermore, we also expect that the description of global quantities may be useful for more complicated systems. Even for phenomena far from equilibrium such as motility-induced phase separation [79], there may be an appropriate set of global variables on the basis of local variables whose space-time evolution may be described by another evolution rule. We will accumulate such examples.
We hope that this paper is a good starting point of many studies.