Storage of Energy in Constrained Non-Equilibrium Systems

We study a quantity T defined as the energy U, stored in non-equilibrium steady states (NESS) over its value in equilibrium U0, ΔU=U−U0 divided by the heat flow JU going out of the system. A recent study suggests that T is minimized in steady states (Phys.Rev.E.99, 042118 (2019)). We evaluate this hypothesis using an ideal gas system with three methods of energy delivery: from a uniformly distributed energy source, from an external heat flow through the surface, and from an external matter flow. By introducing internal constraints into the system, we determine T with and without constraints and find that T is the smallest for unconstrained NESS. We find that the form of the internal energy in the studied NESS follows U=U0∗f(JU). In this context, we discuss natural variables for NESS, define the embedded energy (an analog of Helmholtz free energy for NESS), and provide its interpretation.


Introduction
The basis of equilibrium thermodynamics relies on the existence of the equilibrium state. The equilibrium state can be characterized by a set of appropriate parameters and some kind of energy-based function of these parameters and internal constraints. The constraints allow comparing this function in the state of equilibrium with states of constrained equilibrium [1]. For a monoatomic system, the internal energy U(S, V, N) is a function of three parameters of state, namely entropy S, volume V, and the number of particles N, which fully characterize all thermodynamics changes that can occur in the system. For an unconstrained isolated system, S(U, V, N) is maximized at constant U, V, N with respect to all states obtained by internal constraints.
A prerequisite for any system to become non-equilibrium is a continuous energy flow. This macroscopic flow of energy leads to an increase of the system energy up to the point when the energy flow into the system matches exactly the flow out of the system. At this point, the non-equilibrium steady state is reached. Two parameters characterize the NESS: the flow J U and the internal energy U. We show that U = U 0 * f (J U ), where U 0 is the energy at equilibrium. We make three case studies: (i) a system internally heated between two parallel walls of the same temperature; (ii) a heat flow between two parallel plates of different temperature; and (iii) a Poisseulle flow between two parallel plates.
Non-equilibrium states are ubiquitous in nature and truly equilibrium states are exceptions. However, despite many decades of study, we have not reached the same status of understanding of non-equilibrium states as we have for equilibrium ones. There is no systematic approach for dealing with NESS. Attempts to create such approaches include: minimum/maximum entropy production principle [2], steady state thermodynamics [3], and driven lattice gas systems [4]. The heat flowing into the system is recognized as a source of entropy increase in [5][6][7]. In information theoretical techniques and extended thermodynamics, the heat flow appears as a natural thermodynamic variable in non-equilibrium steady states. The entropy of some ideal systems such as ideal gases, photons, phonons, and ideal harmonic chains, among others, in the presence of a heat flow is studied in [8][9][10][11][12]. However, the energy that has to be stored in NESS has not been recognized as potentially a function of state, from which in principle we could derive all properties of NESS [13].
In this paper, we attempt to address the latter issues. In a recent paper, a quantity is shown to be minimized in steady states for three different systems [13]. This quantity has the dimension of time. In [13], T is shown to coincide with the characteristic time scale of the system energy dissipation immediately after the shutdown of external energy flow. The minimization is demonstrated through introducing a constraint into the system and showing that T for the unconstrained system is always less than in the constrained system. In this paper, we analyze energy storage and T in Systems (i)-(iii) defined above (these systems are different from the ones in [13]) and we arrive at the same conclusions as in [13]. Moreover, we introduce the embedded energy, which is an analog of the Helmholtz free energy for NESS, and provide its interpretation. We point out that, in this paper, we extensively use the temperature profile to obtain the stored energy U − U 0 . The local temperature is defined from the ideal gas law. It would be interesting, however, to consider using effective temperature in non-equilibrium systems and to study their role in energy storage [14,15].

Models and Results
We consider an ideal gas driven out-of-equilibrium by three different ways of energy delivery that are common in physical realizations. In Case (i) the energy is delivered through a homogeneous energy source, in Case (ii) by an external heat flow, and in Case (iii) by an external matter flow.
In steady states, the local energy does not change in time. Therefore, from local energy conservation, we have where σ E ( r) is the local energy source at the position r. Here, we assume the Fourier's law for the local heat flux, where k is the heat conductivity and T( r) is the local temperature. We further assume that in the NESS the ideal gas law is fulfilled locally and that the pressure (and hence the energy density) is constant. From these assumptions, we obtain the following relation between the energy density and the temperature profile, where V is the volume of the system and 0 and T 0 are the energy density and temperature at equilibrium, respectively. In this paper, we denote the corresponding equilibrium value of a variable with a subscript 0. We show the derivation of Equation (4) in Appendix A. From , we define the stored energy as Without performing work (as is the case for our systems), all out-going energy flowĖ out is in the form of heat, which we denote as Φ out , where S is the area through which the heat flows out andn is the unit normal vector. In the steady state, the total energy flow into the system equals the total energy flow out of the systemĖ ss in =Ė ss out = Φ ss out . In the following, we denote the out-going energy flow in the steady state by J U ≡ Φ ss out . By introducing geometrical constraints, the system is partitioned into two subsystems. These constraints do not change the local expressions for J( r) and T( r). In addition, for each subsystem, definitions of the stored energy U i , i = 1, 2 and the out-going heat flow J U i , i = 1, 2 remain the same. On the other hand, the subsystem energy density depends on the constrain in general. However, in all three cases, the number of particles in each subsystem is kept proportional to the volume of the system, i.e., N i /V i = N/V = n 0 . As a result, the expression of i has the same form as Equation (4) (see Appendix A). For the constrained system, we define the stored energy as and the total out-going heat flow as J tot ≡ J U 1 + J U 2 . For every case studied in this paper, we compare the ratio T 1|2 for the constrained system, with the ratio T (see Equation (1)) for the unconstrained system.

Energy Source
In Case (i), we consider a three-dimensional ideal gas placed between two diathermal walls of area A (A → ∞). The walls are kept at temperature T 0 and are fixed at x = ±L. The energy source is distributed homogeneously over the system with σ E ( r) = λ. As internal constraints, we choose a diathermal wall and fix it at x 1 ∈ (−L, L). This wall separates the system into two subsystems 1 and 2 with volumes V 1 = A(L + x 1 ) and V 2 = A(L − x 1 ), respectively. A scheme of the system is shown in Figure 1. Consider first the unconstrained system. As the coordinates y and z do not influence the temperature profile, it is sufficient to consider x−dependence. The temperature profile T(x) is obtained by solving Equation (2), which now has the form −k∂ 2 x T = λ. Using dimensionless variables λ = λL 2 /kT 0 , T(x) = T(x)/T 0 and normalizing x to x = x/L, we obtain Using Equation (4), we find the energy density to be As stated above, the out-going heat flow equals the in-coming energy flow,Ė in = 2LAλ = J U . Combining with Equation (10), we find In the presence of the diathermal wall, the boundary conditions at the constraint are Solving for the subsystem temperature profile with corresponding boundary conditions, we find that T i (x) is not changed by the constraint, i.e., T 1 (x) = T 2 (x) = T(x), in their respective domains. Therefore, we obtain the energy densities as As the total energy source does not change, the out-going heat flow is not changed either, J U 1 + J U 2 = J U = 2LAλ. Together with Equation (12) and (13), we have Now, we compare Equations (11) and (14). The relation reduces to Dividing both sides of Equation (15) by gives where a = Arctanh( We multiply this by both sides of Equation (16) and rearrange the terms to obtain Since (a − x 1 ) 2 ≥ 0, we have verified for this model

Heat Flow
In Case (ii), an ideal gas is in contact with two walls at different temperatures T 1 ≥ T 0 . The walls are of a large area A = H × Z (with height H and width Z) and are placed at x = 0 and x = L (see Figure 2). The steady state is driven by a constant heat flow through the system with no bulk energy supply, i.e., σ E ( r) = 0. For the constraints, we choose an adiabatic wall. We consider two situations. In the first, the wall extends from left to right in a zigzag manner. Its position is given by w 1 (x) = h/2 − hH(x − L/2) where H(x) is the Heaviside function. We refer to this constraint as vertical (Figure 2b). In the second, the wall is a straight line with a slope k, i.e., w 2 (x) = k(x − L/2). We refer to this constraint as linear (Figure 2c). Both constraints are fixed at x = L/2, so that the subsystems are symmetric in shape. Furthermore, they are chosen to ensure a non-zero heat flow in each subsystem. In other words, each subsystem is always in contact with both boundaries. For the unconstrained case, the temperature profile only depends on x. Solving Equation (2) (which is now ∂ 2 x T(x) = 0) with boundary conditions T(x = 0) = T 1 and T(x = L) = T 0 , we have The energy density is then Since we choose T 1 ≥ T 0 , the heat flow passes through the system from left to right. The unit normal vector of the left (right) boundary isn = (−1, 0) (n = (1, 0)). The local heat flux is J( r) = −k(∂ x T( r), ∂ y T( r)). Hence, the heat flow going through the system can be calculated using either of the following expressions, which gives, For the constrained system, the temperature profile depends also on the y-coordinate. It satisfies the equation ∇ 2 T(x, y) = 0 with the following boundary conditions, At x = 0 and L, the system is in contact with the plates. This is represented by the Dirichlet boundary conditions. In addition, since the constraint is an adiabatic wall, we have Neumann boundary conditions at w i (x), i = 1, 2. Finally, at the boundaries far away from the constraint, we expect the effect of the constraint to diminish. In other words, at y = ±H/2, we assume that the heat fluxes are parallel to the x-axis. To ensure this, we need to set H/2 w i (L) and H/2 w i (0) for i = 1, 2. The temperature profiles are obtained numerically using the finite element method. In this method, the system is separated into small domains called mesh and the function is approximated using polynomials [16]. Examples of the contour plot of temperature profiles are shown in Figure 3. After obtaining temperature profiles, the stored energy density is calculated using Equation (4). The total heat flow is obtained using either the left or right boundary according to Equation (21), where y * = w 1 (x = 0) (w 2 (x = 0)) for the vertical (linear) constraint.

Matter Flow
In Case (iii), the ideal gas is flowing between two parallel walls located at y = ±h (see Figure 5). The flow is assumed to be laminar and the fluid incompressible. It is driven by a constant pressure gradient along the x-axis, ∂ x P(x) = −P. Such a flow is known as the Poiseuille flow [17]. Both walls are kept at temperature T 0 . An adiabatic slip wall is introduced as the constraint into the system. It is placed at y = y 1 with 0 ≤ y 1 ≤ 1.
In the steady state, the velocity profile and the temperature profile can be obtained from the Navier-Stokes equation. We note that, due to the presence of the external pressure gradient and since the mass density of the incompressible fluid is homogeneous ρ = ρ 0 , the energy density is not constant throughout the system. We first obtain the velocity profile v = (v(y), 0) from  Figure 5. Schemes of (a) unconstrained and (b) constrained Poiseuille flow. The system is bounded by two plates with a fixed temperature T 0 and area A that are placed at y = ±h. A constant pressure gradient is applied across the system. In (b), the system is divided by an adiabatic slip wall placed at y = y 1 .
where µ is the viscosity. Given the non-slip conditions at the boundaries v(±h) = 0, we find Secondly, from the momentum equation, we obtain the dissipation density function φ = µ(∂ y v) 2 = P 2 y 2 /µ. The dissipation density function governs the rate at which the mechanical energy of the flow is converted to heat. The out-going heat flow J U is given by, where A is the area of the plates and V = A × 2h is the volume of the system. Moreover, we assume that the heat transfer obeys the Fourier's law. From (Equation (2)), together with boundary conditions T ±h = T 0 , we obtain the temperature profile, Finally, we assume that the internal energy locally obeys the ideal gas law given by Equation (A4). The total energy of the system consists of the kinetic energy and the internal energy, where the number density n 0 = ρ 0 /m, with m the mass of a single atom or molecule. Combining Equations (27), (30) and (31), we obtain For the constrained system, additional boundary conditions are dv 1 /dy | y=y 1 = 0, dv 2 /dy | y=y 1 = 0, dT 1 /dy | y=y 1 = 0 and dT 2 /dy | y=y 1 = 0. Following the same method, we find the velocity profiles as and the temperature profiles as From these equations, we obtain Comparing T and T 1|2 , the relation reduces to Analysis shows that T ≤ T 1|2 .

Energy Density as Function of Heat Flow
It is interesting to note that, for all the above studied models, the steady state energy density is a product of the equilibrium energy density 0 and a dimensionless function of the heat flow J U . For the ideal gas system with a homogeneous energy supply, where J U = 2LAλ, can be written as (compare Equation (10) Next, for the heat flow model, with J U = (T 1 − T 0 )Ak/L (compare Equation (22)), the steady state energy density can be expressed as Lastly, for the matter flow model, with J U = 2AP 2 h 3 /3µ (compare Equation (27)), the steady state internal energy density can be expressed as (compare Equation (31)), Thus, in all studied steady states, we find U = U 0 * f (J U L/(AkT 0 )).

Conclusions
We use the ideal gas model with three different energy delivery methods to test the hypothesis that ∆U/J U is minimized in steady states. The results in all models confirm that ∆U/J U ≤ ∆(U 1 + Further, in all studied steady states, we find U = U 0 * f (J U L/(AkT 0 )) and therefore J U is a parameter of NESS. By making a Legendre transform of U with respect to J U , we get an analog of the Helmholtz free energy for NESS, especially since J U /T 0 is the entropy flow leaving the system through the wall at temperature T 0 . We introduce a quantity U − (dU/dJ U )J U = U * , which we call the embedded energy, since it is the stored energy minus the outflow of energy in the characteristic time τ = dU/dJ U . Thus, U * represents the part of the energy that must stay in the system for all times to keep the outflow of energy, while τ J U is the energy that constantly flows through the system in time τ.

Conflicts of Interest:
The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

Appendix A
Here, we provide the calculation of the energy density from the temperature profile T( r) with r = (x, y, z).
The system has a fixed number of particles N and a fixed volume V. It obeys the ideal gas law, where P is the pressure, n = N/V is the particle number density, k B is the Boltzmann constant, and is the (internal) energy density. For steady states in Cases (i) and (ii), we assume that the pressure (and hence the energy density) is homogeneous across the system, that is, The energy density can be obtained by observing that where n 0 is the number density at equilibrium and we used n 0 V = V d 3 rn( r) = N. We denote the equilibrium variables with subscript 0. The energy density is thus When introducing the constraints, the subsystems are separated in such a way that N i /V i = n 0 . For each subsystem, we have V i d 3 rn i ( r) = N i and, thus, Therefore, the expression for i is of the same form as , . (A8)