# Ensemble dependence of Critical Casimir Forces in Films with Dirichlet Boundary Conditions

###### Abstract

In a recent study [Phys. Rev. E 94, 022103 (2016)] it has been shown that, for a fluid film subject to critical adsorption, the resulting critical Casimir force (CCF) may significantly depend on the thermodynamic ensemble. Here, we extend that study by considering fluid films within the so-called ordinary surface universality class. We focus on mean-field theory, within which the OP profile satisfies Dirichlet boundary conditions and produces a nontrivial CCF in the presence of external bulk fields or, respectively, a nonzero total order parameter within the film. Our analytical results are supported by Monte Carlo simulations of the three-dimensional Ising model. We show that, in the canonical ensemble, i.e., when fixing the so-called total mass within the film, the CCF is typically repulsive instead of attractive as in the grand canonical ensemble. Based on the Landau-Ginzburg free energy, we furthermore obtain analytic expressions for the order parameter profiles and analyze the relation between the total mass in the film and the external bulk field.

## I Introduction

Confining a critical fluid by parallel walls gives rise to a critical Casimir force (CCF) acting on the bounding surfaces Fisher and de Gennes (1978); Krech (1994). Here we consider fluids belonging to the Ising bulk universality class (UC), which, accordingly, are described by a one-component order parameter (OP) field . The bulk UC splits up into several surface UCs, describing further universal properties induced by the surfaces Diehl (1986, 1997); Brankov et al. (2000). In a classical fluid, the constituent molecules are generically attracted towards an immersed solid surface. This attraction can be either strong or weak compared with the liquid-liquid interaction. Accordingly, for a one-component fluid the surfaces have a preference either for its liquid phase (in the case of a strong substrate) or the vapor phase (in the case of a weak substrate), whereas for a binary liquid mixture the walls attract that phase which is rich in the species preferred by the surfaces. Near the critical point, this attraction gives rise to the phenomenon of critical adsorption, which, in the limit of infinitely strong adsorption (surface field ), is described by the so-called normal surface UC Fisher and Nakanishi (1981); Liu and Fisher (1989); Flöter and Dietrich (1995). Fluids show also an enhanced molecular order near a solid surface Diehl (1986); Flöter and Dietrich (1995), which is modeled field-theoretically by a so-called surface enhancement parameter . The limit (for finite adsorption strength ) defines the so-called ordinary surface UC, in which the OP effectively satisfies Dirichlet boundary conditions.

While critical fluids are typically strongly adsorbed at container walls Gambassi et al. (2009), by suitable preparation of the surfaces it is nevertheless possible to approach the limit of weak adsorption, corresponding to the ordinary surface UC. In Ref. Nellen et al. (2009), this has been achieved by chemical treatment of the surface, while in Refs. Sprenger et al. (2006); Tröndle et al. (2009, 2010); Gambassi and Dietrich (2011) surface patterning has been used.

The CCF stems from residual finite-size contributions of the free energy of the film. Remarkably, as has been shown in Refs. Gross et al. (2016, 2017), the amplitude and the scaling function of the CCF depend not only on the bulk and the surface UC, but also on the thermodynamic ensemble under consideration. In fact, CCFs are typically studied for fluid films which can exchange particles with their environment—a situation which realizes the grand canonical ensemble. However, global OP conservation, which is applicable for the canonical ensemble, can induce drastic changes of the CCF Gross et al. (2016, 2017). In the present study, building on Ref. Gross et al. (2016) (where critical adsorption has been investigated), we consider Ising-type fluid films within the ordinary surface UC, subject to a global OP constraint. We focus on mean-field theory, within which the effects of fluctuations are neglected and the CCF is a consequence of the presence of a spatially varying OP profile across the film.

In the grand canonical ensemble, a nonzero external bulk field acting in the film does induce a nontrivial OP profile. In the canonical ensemble, instead, a nonzero value of the total integrated OP, henceforth called the mass, is imposed:

(1) |

Here denotes the transverse area of the film, its thickness, and the associated transverse coordinate. We generally assume the film to be homogeneous in the remaining, lateral directions. Henceforth we consider all extensive quantities, such as , as quantities per transverse area , i.e., . We find that the OP constraint in Eq. (1) can change, inter alia, the character of the CCF from attractive in the grand canonical case to repulsive in the canonical case.

In passing, we recall that, for a critical fluid film within the ordinary surface UC, the critical temperature is shifted from its bulk value to . For Dirichlet boundary conditions and vanishing external fields , the OP profile vanishes above the film critical point, i.e., for temperatures . CCFs for Ising-type systems in the ordinary surface UC (including crossover effects to the normal surface UC) have been previously studied within the grand canonical ensemble in Refs. Krech and Dietrich (1991, 1992a, 1992b); Schmidt and Diehl (2008); Hasenbusch (2011); Diehl and Schmidt (2011); Mohry et al. (2010); Vasilyev et al. (2011); Vasilyev and Dietrich (2013).

In Sec. II, we define the general scaling variables required for the description of the universal critical properties and outline the scaling relations expected for the OP profile. We furthermore introduce the Landau-Ginzburg model which is analyzed in the remaining part of this study. The OP profile resulting from the Landau-Ginzburg model within mean-field theory is determined perturbatively in Sec. III and fully via numerical studies in Sec. IV. The associated relation between the total mass and the external bulk field is analyzed separately in Sec. V. In Sec. VI, the CCF is studied analytically within linearized MFT and numerically within full MFT, focusing on ensemble differences. In Sec. VII the predictions of MFT are compared to Monte Carlo simulations of the three-dimensional Ising model.

## Ii Preliminaries

### ii.1 Scaling behavior

Here, we summarize the general scaling behavior expected for the OP profile and the CCF in a film of thickness . In the following we focus on the so-called ordinary fixed point, at which and, accordingly, the dependence of the scaling functions on drops out. According to standard finite-size scaling arguments, the universal properties of a critical film are expected to be controlled by the following set of scaling variables:

(2a) | ||||

(2b) | ||||

(2c) | ||||

(2d) |

where

(3) |

is the mean mass density of the film, , , and are standard bulk critical exponents, and

(4) |

is the reduced temperature relative to the bulk critical temperature . In the case of a one-component fluid the (dimensionless) external bulk field describes the deviation of the chemical potential from its critical value in the bulk, while for a binary liquid mixture, represents the deviation of the difference in the chemical potentials of the two species A and B from its bulk critical value: . The quantities and (as well as , which we include here for completeness) denote non-universal amplitudes defined in terms of the (bulk) correlation length at zero bulk field and at zero reduced temperature:

(5a) | |||||

(5b) |

The value of is different for , but the amplitude ratio forms the universal number in and in spatial dimensions Pelissetto and Vicari (2002). Except for Sec. IV, we focus on the supercritical regime and therefore in the scaling relations we use solely . The non-universal amplitude is defined in terms of the bulk OP , which, near criticality, behaves as

(6a) | |||||

(6b) |

in the case of a vanishing external field and a vanishing reduced temperature , respectively.

The OP profiles in the grand canonical and the canonical ensemble fulfill the following scaling relations Binder (1983); Diehl (1986); Privman (1990); Brankov et al. (2000); Krech (1994):

(7a) | ||||

(7b) |

where are the corresponding universal scaling functions. In order to simplify the notation, we henceforth drop the superscripts (c) and (gc) on and . The scaling variable in Eq. (2d) is related to the scaling function via

(8) |

The general scaling behavior of the CCF is discussed in Sec. VI.

### ii.2 Model and boundary conditions

We aim at determining the order parameter profile between two parallel plates, located at and subject to the constraint of a specified total mass [see Eq. (1) and recall that here and in the following is considered per area ]. The canonical Landau-Ginzburg (LG) free energy functional for films, in units of per transverse area of the plates, is given by

(9) |

The integral represents the bulk contribution, whereas the terms are surface enhancements giving rise to Robin-type boundary conditions Diehl (1986) on — see Eq. (12) below. Within MFT, the coupling constants and are given by and , where is the reduced temperature [Eq. (4)] and the amplitudes and are defined in Eqs. (5a) and (6a). Within MFT, one has . Equilibrium states minimize Eq. (9), subject to the constraint in Eq. (1). In the grand canonical ensemble the LG functional for films (per and area ) reads

(10) |

which is to be minimized with respect to , taking for the external bulk field (i.e., the chemical potential) a value such that Eq. (1) is obeyed. Minimization of the grand canonical energy functional leads to the Euler-Lagrange equation (ELE)

(11) |

subject to the boundary conditions

(12) |

induced by the surface enhancement terms. In what follows, we shall study the limits , for which Dirichlet boundary conditions emerge.

Within MFT, the finite-size scaling variables defined in Eq. (2) turn into

(13) |

in terms of which in Eq. (10) can be expressed as

(14) |

The non-universal amplitude is given by

(15) |

in terms of the amplitudes of the correlation length and the bulk OP [see Eqs. (6a) and (5a)]. The dimensionless form of the ELE, following from Eqs. (11) and (12), reads

(16) |

with the corresponding Dirichlet boundary conditions (obtained in the limits )

(17) |

Equations (16) and (17) are independent of the plate separation and the coupling constant , because these variables can be scaled out such that they appear as prefactors in Eq. (14). In general, via Eq. (15), and can be expressed in terms of experimentally accessible critical amplitudes.

## Iii Perturbative Mean Field Analysis

In order to make analytical progress, we address the nonlinear term of the ELE in Eq. (16) perturbatively by introducing a parameter (eventually to be set to unity):

(18) |

This equation must be solved subject to the Dirichlet boundary conditions in Eq. (17) and under the constraint [Eq. (8)]

(19) |

In a first step, we solve Eq. (18) without this constraint by carrying out perturbation theory in terms of powers of , with the series expansions

(20) |

The boundary conditions from Eq. (17) hold for each term . Concerning the expansion of the mass constraint in Eq. (19), we choose

(21) |

where .

As a side remark, one infers from the structure of the ELE that, if is a solution of Eq. (18) with parameters and , then will be a solution for the parameters and . Thus, the total mass is an odd function of the bulk field , i.e., . This feature is illustrated in Fig. 1 for the full, numerical (non-perturbative) solution of Eq. (18), which must hold also at each perturbative order.

### iii.1 Solution at

At this order, Eq. (18) yields

(22) |

with the solution

(23) |

In contrast to the case of critical adsorption considered in Ref. Gross et al. (2016), the lowest order MFT solution for Dirichlet boundary conditions is well-behaved near the bulk critical point. This is revealed by a series expansion for small , yielding . By using Eq. (13), Eq. (23) can be written in terms of dimensional variables:

(24) |

which will be useful for the analysis presented in Sec. VI.1.

Implementing now the constraint in Eq. (21) selects and fixes, at this order, the value :

(25) |

where the last expression exhibits the asymptotic scaling behavior close to the bulk critical point and for thick films, respectively. Inserting Eq. (25) into Eq. (23) gives the contribution to the constrained order parameter at this order:

(26) |

The asymptotic scaling of this expression,

(27) |

shows that at bulk criticality the lowest order MFT contribution for Dirichlet boundary conditions is a parabolic profile. In turn, away from criticality, the (spatially constant) solution must vanish due to the boundary conditions, which shows that if . Consequently, in Eq. (25) must also vanish away from criticality. Finally, expressing Eq. (26) in terms of dimensional variables, one finds the constrained profile

(28) |

which indeed satisfies the relation .

### iii.2 Solution at

To linear order in , Eq. (18) gives

(29) |

The solution of this differential equation vanishes in the limit . (The full expression is cumbersome and is not shown here.) Implementing the constraint of Eq. (21), one finds the following corresponding specific expression :

(30) |

which exhibits the asymptotic scaling behavior

(31) |

From this the constrained profile for very small and very large can be calculated:

(32) |

At bulk criticality, a polynomial solution obeying the boundary conditions in Eq. (17) is obtained. As it was the case for the contribution , the constrained profile vanishes away from criticality.

The perturbative solution of the ELE to is reported in Appendix A.

## Iv Comparison of perturbative MFT solutions with exact and numerical results

In this section we compare the leading perturbative solution at order with numerical solutions of the full, nonlinear ELE (18). In the case of zero external field, the full solution can be computed analytically (see Sec. IV.1 below). In Sec. IV.2 we consider the unconstrained solution for a given pair of parameters , and compare it with given by Eq. (23). Therefore, in Sec. IV.3 we impose the constraint on the total mass and regard the corresponding solution as a function of the independent parameters . The latter is compared with as given by Eq. (26).

### iv.1 Exact analysis for and location of the film critical point

For an exact expression for the order parameter profile can be obtained in closed form in terms of elliptic functions Gambassi and Dietrich (2006). According to Eq. (16), the associated ELE is

(33) |

subject to the boundary conditions . Beside the trivial solution , there is a non-vanishing solution for , where

(34) |

denotes the scaled reduced temperature (relative to the bulk critical point) of the film critical point. Here and in the following, when considering the regime , i.e., , we define as

(35) |

One finds that

(36) |

where

(37) |

is the complete elliptic integral of the first kind, is the elliptic modulus, determined implicitly by through , and sn is Jacobi’s elliptic sine (see Refs. Gradshteyn and Ryzhik (2014); Olver et al. (2010) for more details). As shown in Fig. 2, the numerical solution of Eq. (33) perfectly matches the exact solution given in Eq. (36).

### iv.2 Unconstrained profiles

#### iv.2.1 Profiles for

In Fig. 3, a comparison is shown of the unconstrained profiles obtained numerically with the perturbative approach at leading order. The perturbative solution [Eq. (23)] deviates significantly from the numerical solution for large values of the bulk field , with the largest deviations being localized in the middle of the film (i.e., ). Close to the film boundaries at , the inaccuracy of the perturbative solution is mitigated by the fact that the Dirichlet boundary conditions are satisfied for all values of .

#### iv.2.2 Profiles for

The approach outlined above can be followed also for , and in principle the entire phase diagram can be explored. However, the same qualitative behavior encountered for occurs also for . In general, the strongest inaccuracy is observed for (as the phase-separating regime is approached) and for large values of (where nonlinear effects become more dominant due to the term in the ELE).

### iv.3 Constrained profiles at and

Here we consider the constrained profiles obtained by numerically solving the ELE [Eq. (16)] and compare them with the first-order perturbative solution [Eq. (26)]. In Fig. 4, the two cases and [Eq. (34)] are examined, where the latter corresponds to the film critical point. It is interesting to note that for the perturbative profile is not singular, but reduces to a particularly compact form:

(38) |

## V Phase diagrams, equation of state, and scaling

Here we explore the magnetization phase diagram, the equation of state , and, in particular, we compare the film behavior with the one corresponding to the bulk. Exact numerical results are discussed in Sec. V.1, while the validity of the perturbative MFT results is studied in Sec. V.2. In Sec. V.3 we show that the near-critical behavior of the mass can be captured by simple scaling arguments. This scaling behavior can even be applied to the order parameter profiles themselves, as will be discussed in Sec. V.4.

### v.1 Exact numerical results for the mass

While the nonlinear ELE in Eq. (18), subject to Dirichlet boundary conditions, can be solved by standard numerical methods for (i.e., above phase separation in the film) and for sufficiently small (for which nonlinear effects are not too strong), these methods typically become inaccurate outside these regimes, where gradients of the profile can be large. This issue can be addressed by solving the ELE via the so-called symplectic integration method Ruth (1983); Hairer et al. (2010); Gross et al. (2016), which, by construction, yields a spatially constant pressure in equilibrium. Essentially the ELE in Eq. (16) is equivalent to the Hamiltonian “equations of motion”, and the algorithm conserves the Hamiltonian density , which, in turn, allows one to directly extract the film pressure (see, c.f., Eq. (63)). This method has the advantage that it avoids the (inaccurate) numerical computation of . The order parameter profile obtained this way for a pair () of scaling variables can be integrated numerically in order to determine the corresponding mass. The results of this procedure are shown in Fig. 5. The shift of the critical point in the film is clearly visible, as is the symmetry . The super-imposed bulk diagram was obtained by solving Eq. (18), without the gradient term, for .

### v.2 Comparing exact and perturbative results for the mass

#### v.2.1 Mass as function of an external field at

The lowest order perturbative MFT solution for the mass [Eq. (25)] is linear in at bulk , i.e., . In Fig. 6 we compare this result (solid line) with the exact mass computed from the numerical solution of the nonlinear MFT (dots). The lowest order MFT result starts to deviate significantly from the exact result at , whereas the numerical solution gradually approaches the bulk critical behavior with within MFT.

#### v.2.2 Mass as function of with

In the absence of the external magnetic field one can use the exact solution (Eq. (36)) for the study of the mass:

(39) |

The elliptic modulus entering into the exact solution is the positive root of the implicit equation , with defined in Eq. (35). The integration in Eq. (39) can be carried out in closed form by using elementary properties of elliptic functions Gradshteyn and Ryzhik (2014); Olver et al. (2010):

(40) |

From this result one can easily extract the asymptotic behavior of the mass. In particular, we proceed to analyze Eq. (40) for (i) close to film criticality , i.e., for , and for (ii) extreme subcritical temperatures .

(i) : It is convenient to parametrize the deviation from the critical point as

(41) |

where is the film-analogue of the bulk reduced temperature as introduced in Eq. (4). We note that for from Eq. (35) we have , while instead at we recover , as expected. It follows furthermore that . We focus on the regime . Since in this limit (see Eq. (41)), the implicit equation can be substituted by its Taylor expansion around the desired value of . The corresponding small-modulus expansion of the complete elliptic integral is , which implies the following expression for :

(42) |

The solution , corresponding to , is trivially reproduced. For , Eq. (42) gives, to leading order, , which is valid for approaching from below. Higher order corrections can be obtained by iterating this procedure. Inserting this result into Eq. (39) and using the fact that , one obtains the following scaling behavior:

(43) |

with the exponent . Equation (43) is valid in the asymptotic regime where successive corrections , characterized by an exponent , vanish faster than .

(ii) : Since , the roots of accumulate towards . Writing for certain , one has , where is a polynomial in of degree . Hence for large the mass is approximately given by . In order to identify the small parameter in terms of , we note that if the elliptic modulus approaches unity, , so that the implicit equation for exhibits the asymptotic behavior . Accordingly, the mass is , where the last term is negligible because and . This renders the asymptotic result

(44) |

To summarize, we have derived the analytical expression of the mass in the absence of an external field, and its asymptotic behavior close to film criticality () and far from criticality in the two-phase region (). As shown in Fig. 7, the approximate expressions agree well with the analytical result in Eq. (40).

### v.3 Widom scaling for the mass

It is well known Gambassi and Dietrich (2006); Nakanishi and Fisher (1983) and explicitly demonstrated in Sec. IV.1, that in the film geometry the presence of two confining walls induces a shift of the bulk critical point from to . In the present section we discuss in detail the mean-field critical behavior around , resulting from Eq. (18).

It is useful to recall the essential ideas of the static scaling hypothesis, as originally formulated by Widom Widom (1965); Stanley (1971). The film critical point is located at , where [see Eq. (41)] is the reduced temperature of the film relative to . Instead of considering the order parameter profile inside the film, here we are interested in the mass . In the critical region of the film, for a vanishing bulk field one expects the scaling behavior

(45) |

We note that, according to Eq. (41), corresponds to . The critical isotherm follows as

(46) |

The above relations can be considered as a definition of the critical exponents and and of the non-universal amplitudes and . According to the scaling hypothesis, in the near-critical region around the equation of state fulfills a homogeneity relation of the form

(47) |

where are a pair of universal scaling functions and is called the gap exponent Kadanoff et al. (1967); Stanley (1971). Various sections of the phase diagram in the scaling region lead to curves of the type shown in Fig. 8(a). A suitable rescaling of the thermodynamic variables and as prescribed by Eqs. (45) and (46) results in a data collapse onto two single master curves corresponding to the scaling functions (see Fig. 8(b)).

In the previous subsection we have established Eq. (45) in the form of Eq. (43), leading to and to the non-universal amplitude . On the other hand, the results of the complete numerical analysis, shown in Fig. 8, confirm Eq. (46); in fact the critical isotherm in the scaling region near can be approximated well by Eq. (46) with the critical exponent and the non-universal amplitude We thus recover for the gap exponent and obtain an excellent data collapse.

To summarize, our analytical and numerical analysis recovers the expected mean field critical exponents for the film critical point. We remark that the maximum value of the critical profile (in the center of the film) has the same scaling behavior as the total mass, , with and .

This analysis reveals explicitly that, as expected, within MFT the bulk transition in spatial dimension exhibits the same scaling behavior and the same critical exponents as its counterpart in the film which, asymptotically, behaves as an effectively -dimensional system. The inability to capture this actual dimensional crossover is a well-known shortcoming of many analytical approaches, i.e., MFT and beyond Krech (1994); Diehl (1997); Dohm (2009); Diehl et al. (2014), whereas simulations can deal with this issue successfully.

### v.4 Magnetization profiles in the near-critical region: insights from Widom scaling

Here we consider the case and . Since at criticality the profile vanishes, the ELE in Eq. (18) in the vicinity of the film critical point, i.e.,

(48) |

can be approximated by the linearized equation

(49) |

because the cubic term is smaller that the linear terms. Equation (49) with Dirichlet boundary conditions is solved by

(50) |

However, the amplitude cannot be fixed by Eq. (49), because the cubic term has been neglected. Nonetheless, we can determine by considering a suitable limit of the exact solution. For we can use the reduced temperature from Eq. (41). We recall that within this limit the elliptic modulus is and that for vanishing the Jacobi elliptic function reduces to a standard sine function. Therefore in the limit Eq. (36) produces exactly

(51) |

Spatial integration yields the mass for . The behavior of the scaled bulk field is less obvious. In the previous section we noted a scaling behavior for the maximum value of the magnetization profiles, namely . The same behavior extends, with remarkably good agreement with the numerical results of Fig. 3, also to . We find that Eq. (51) follows an analogous scaling, i.e.,

(52) |

Combinining Eq. (51) and Eq. (52), in the scaling region around the film critical point we have

(53) |

with a scaling function . Since , one has

(54) |

where , so that Eq. (47) is recovered, for which we identify as . Thus, the scaling functions computed for the mass equation of state in Fig. 8 capture well the spatially integrated order parameter profiles in the near-critical region.

## Vi Critical Casimir Force

In this section we study the critical Casimir force (CCF) in the grand canonical and the canonical ensemble, for a film subject to Dirichlet boundary conditions. We therefore briefly recall the general definitions and protocols for computing the CCF, as set out in Sec. III of Ref. Gross et al. (2016).

In general, the equilibrium CCF provides the derivative of the residual free energy, or, in terms of the stress tensor, quantifies the change of the free energy of the film upon shifting the position of the boundaries. In the first case, one has

(55) |

where we have decomposed the free energy of the film according to

(56) |

in terms of the bulk pressure , the surface free energy , and the residual free energy (all per transverse area and ). Generally, the bulk term scales and the surface term , while the residual terms vanish exponentially for (see, e.g., Refs. Gross et al. (2016); Privman (1990)).

In the second case, the CCF is the difference between the film pressure and the pressure of the surrounding bulk medium in which the film is immersed:

(57) |

The bulk pressure is naturally defined as

(58) |

where the limit is performed by keeping fixed the relevant thermodynamic control parameters (i.e., the chemical potential for the grand canonical ensemble, and the mass density for the canonical ensemble) ^{1}^{1}1In the canonical ensemble the actual thermodynamic control parameter is the total mass . However, as discussed in detail in Ref. Gross et al. (2016), for the purpose of determining the finite-size limit in Eq. (58), instead the mass density should be kept fixed..

The CCF (per transverse area and ) in the grand canonical and the canonical ensemble takes the following scaling form Krech (1994); Brankov et al. (2000); Gross et al. (2016):

(59a) | ||||

(59b) |

where and are scaling functions, which will be determined below for and within MFT, i.e., for spatial dimensions.

Instead of using Eq. (58), the film pressure can equivalently be obtained from the stress tensor :

(60) |

where is computed from the order parameter profile minimizing (Eqs. (9) and (10)). Note that here we have assumed the boundaries of the film to be normal to the -direction. Analogously, the bulk pressure can be obtained from the corresponding bulk order parameter at equilibrium, . Therefore Eq. (57) allows one to compute without explicitly evaluating derivatives of free energy functionals. In the grand canonical ensemble, the definitions of in Eqs. (55) and (57) yield equivalent results, whereas differences may appear due to additional surface contributions in the canonical ensemble Gross et al. (2016).

A core result of Ref. Gross et al. (2016) is that the stress tensor in the canonical ensemble can be computed using a grand canonical stress tensor in which the chemical potential takes the value , satisfying the mass constraint in Eq. (1),

(61) |

in terms of the solution of the ELE. By construction, this yields equal film pressures in the two ensembles, . In the grand canonical ensemble, the mean field stress tensor corresponding to the free energy functional in Eq. (10) is Krech (1994)

(62) |

giving rise to the film pressure

(63) |

where the dimensionless variables from Eq. (13) have been re-introduced; is given by Eq. (15). In turn, the bulk pressure in the grand canonical ensemble,

(64) |

is obtained by solving the bulk equation of state (i.e., the ELE without gradient terms),

(65) |

in order to find the spatially constant solution