First order phase transitions and the thermodynamic limit

We consider simple mean field continuum models for first order liquid–liquid demixing and solid–liquid phase transitions and show how the Maxwell construction at phase coexistence emerges on going from finite-size closed systems to the thermodynamic limit. The theories considered are the Cahn–Hilliard model of phase separation, which is also a model for the liquid-gas transition, and the phase field crystal model of the solid–liquid transition. Our results show that states comprising the Maxwell line depend strongly on the mean density with spatially localized structures playing a key role in the approach to the thermodynamic limit.


Introduction
In this paper we revisit the topic of equilibrium first order phase transitions and elaborate on the origin of the famous Maxwell or 'equal areas' construction that applies in the thermodynamic limit (TL), i.e. for infinite systems. In particular, we examine in detail how this construction emerges as the system size becomes larger and larger, thereby gaining additional insight into this construction. According to Ehrenfest's classification, a first order phase transition is characterized by the appearance of a discontinuity in a first derivative of the free energy with respect to some thermodynamic variable, e.g. for the solid-liquid phase transition, the derivative with respect to pressure becomes discontinuous, implying a jump in the density. For a second order transition all first derivatives are continuous and a discontinuity occurs in second derivatives. More modern classifications define first order transitions as transitions that involve latent heat or as transitions where an order parameter changes discontinuously [1,2].
At a first order phase transition two different phases (e.g. a gas and a liquid with distinct densities) can coexist and the characteristics of the coexisting states can be calculated by performing the Maxwell construction on the relevant thermodynamic quantity plotted as a function of the order parameter. This can for example be the pressure as a function of the volume or an appropriate free energy specified by the constraints on the system, such as conservation of particle number at fixed volume, or similar. Strictly speaking, the Maxwell construction is only valid in the TL. Likewise, the related discontinuities only arise in this limit (i.e. at infinite system size and particle number) and at equilibrium, after all transients have decayed. The equilibrium states and their transitions are then represented in the form of state diagrams showing various thermodynamic quantities as a function of a relevant control parameter and phase diagrams that show the location of the various stable phases in two-or higher-dimensional parameter space. The Cahn-Hilliard equation can be written in the form of conserved gradient dynamics for a scalar order parameter field f [16,32] where Q(f) is a positive mobility function (not relevant for steady states), and the underlying free energy functional is [17] where V is the domain volume and dx is a volume element in V. Here, the field f(x, t) corresponds to a local concentration, i.e. a scaled linear combination of local particle number densities. The first term in (2) captures the energetic cost of interfaces (κ0) and the local free energy density is the double-well potential , 3 2 4 obtained on making a Taylor expansion of the true potential about the critical point [33]. Here b>0 while the parameterã can change sign. We denote the temperature at the critical point as T c and write˜( ) = a a T T c , a>0. The critical concentration is f c =0. Note that, for uniform concentration systems, the free energy is just expression (3) multiplied by the volume.
Approximating the mobility Q by a constant Q c , proportional to the diffusion coefficient, one obtains from equation (1) the standard form of the Cahn-Hilliard equation [16] [ The term in the square brackets represents (in general, in nonequilibrium) a nonuniform chemical potential that is made up of an interfacial contribution proportional to the Laplacian Δf and the local term ∂ f f (f).
To study steady states, i.e. time-independent uniform or nonuniform concentration profiles, we set ∂ t f=0 in equation (4) and obtain after two integrations where the integration constant μ represents a Lagrange multiplier to enforce the constraint that the average concentration ( ) takes a specified value. This constraint reflects the fact that the total number of particles of each of the two species in the binary mixture is separately conserved. This conservation stems from the form of the dynamics in equation (1). In the following f 0 is used as the relevant control parameter for obtaining stationary solutions. The integration constant after the first integration is set to zero as appropriate for systems with no flow across the boundaries. See, e.g. [34][35][36] for situations where this condition is not fulfilled. We note in passing that in situations where the total concentration is not controlled, the parameter μ becomes a relevant control parameter representing an external field or imposed chemical potential. However, this changes the properties of the associated bifurcation diagram [26]. Strictly speaking, μ is actually a scaled chemical potential difference, since the local concentration f is a scaled difference in the local densities, although below we refer to it simply as the chemical potential.
Note that the Cahn-Hilliard model can also be thought of as a simple model for gas-liquid phase separation; in this case the order parameter f represents a scaled density change from its critical value.

PFC model
The conserved Swift-Hohenberg equation with cubic nonlinearity, also known as the PFC model [19,37], provides the simplest phenomenological microscopic mean-field continuum description of the dynamics of the transition between a liquid state and a crystalline state. This local (i.e. partial differential) equation may be derived from a truncated gradient expansion in the DDFT description of an undercooled system undergoing crystallization [19,[38][39][40]. The governing equation also takes the form of conserved gradient dynamics for a scalar order parameter field f, as in equation (1), this time with the underlying free energy functional that has higher order spatial derivatives than the Cahn-Hilliard free energy functional (2). The above two forms of the free energy are related by partial integrations with appropriate boundary conditions. Note that the coefficient of | | f  2 is now negative and so favors gradients in the order parameter, while the system is regularized by the strictly positive higher order term that limits the steepness of the variations in f. The parameter q represents the dominant wave number and r is the undercooling. This model has a critical point at r c =0, f c =0, and so is expected to fail before this point is reached [38]: in experiments the freezing transition is first order. The origin of this failure lies in the approximations made in deriving the PFC, but for r=r c the PFC model provides a qualitatively correct description of the freezing transition region. The resulting kinetic equation is of sixth order and we refer to this conserved Swift-Hohenberg equation as the PFC model. Here Q c again represents a constant mobility. Other sign conventions as well as other nonlinearities are discussed, e.g. in [19,20,41,42]. As before, steady states of (7) are studied by setting ∂ t f=0. After two integrations we obtain where the integration constant μ again represents a Lagrange multiplier for particle number conservation, i.e. the chemical potential. Thus, the steady states of the PFC equation correspond to steady states of a nonconserved Swift-Hohenberg equation. However, the chemical potential μ plays an important role through the properties of the associated bifurcation diagram [20,26] unless μ is set to zero as appropriate for studies of the nonconserved system [8,43].

Numerical approach
For both of the models studied here, the Cahn-Hilliard and the PFC equations, we are only interested in the various possible steady states and their bifurcations. Branches of steady state solutions are determined using pseudo-arclength path continuation techniques [26,44,45] employing the packages AUTO07 [46,47] for the Cahn-Hilliard equation in 1D and PDE2PATH [48,49] for the Cahn-Hilliard equation in 2D. The latter is used for the PFC model in both 1D and 2D. In the 2D case, we choose numerical domains  num , on which we apply the continuation method, that are fractions of the physical domain  shown in the figures. This fraction is determined by the symmetry of the state considered. This procedure lowers the computational cost, but the results of the accompanying numerical linear stability analysis have to be interpreted with care as only perturbations fulfilling the assumed symmetries are admitted. As a result bifurcations that break the symmetry of the state are not detected. In particular, in section 3 we employ the Cahn-Hilliard equation to study liquid-liquid phase decomposition in 1D and 2D domains with periodic boundary conditions. Consequently the critical domain size is defined as L c =2π/q + with q + being the upper limiting wave number of the band of wave numbers with positive growth rate in the linear regime. For all calculations, equation (5) is rescaled such that q + =1, i.e. L c =2π. In section 3.3 we consider states on a square domain with periodic boundary conditions. For large domains we use an adaptive mesh to guarantee the convergence of the results and improve computational performance. Since the drop and stripe states under consideration fulfill reflection symmetry in both directions, the calculation is done on one quarter of the physical domain imposing no-flux boundary conditions. Furthermore, in the case of stripe states we use the translation invariance of these states in one direction and directly calculate the associated branches in the 1D system. However, this procedure fails to capture the linear stability of the stripe state with respect to transverse perturbations. Instead, we use the branches of modulated stripes that connect translation-invariant stripe and drop states obtained in 2D calculations to deduce the full linear stability properties of the translation-invariant stripes.
In section 4 we employ the PFC equation to consider crystallization in 1D and 2D. Here, the critical domain size is defined as L c =2π/q with q representing the dominant wave number, where the maximal growth rate first crosses zero, i.e. at the instability onset. The model is scaled so that q=1 and therefore once again L c =2π. In section 4.3 we analyze localized and periodic states with hexagonal symmetry (see, e.g. figure 12 below). This allows us to use a numerical domain  num defined as a 1/12 angular section (angles from 0 to π/6) of  with no-flux boundary conditions. This implies that we impose a π/3 rotational symmetry on the states considered as well as reflection symmetry with respect to the median lines. At small amplitude the periodic hexagonal states can be thought of as resulting from the superposition of three harmonic modes, all with the same wave number q=1 but orientations rotated by 2π/3. In this case it is more convenient to characterize the domain size by the length of the hexagon side. The critical side length for pattern formation is . Since the Cahn-Hilliard and PFC models are both continuity equations and we only consider steady states with periodic or no-flux boundary conditions, we can use the integrated equations (5) and (8), respectively. There are then two possibilities for studying the dependence of the various steady states on the mean concentration f 0 [26]. Direct continuation in the control parameter f 0 is incorporated through an integral condition. This additional equation forces the use of an additional free parameter (the Lagrange multiplier μ) which is adapted as part of the continuation procedure. Alternatively, the parameter μ is employed directly as a control parameter without the need to include any further condition on the continuation procedure. However, an integral has to be evaluated at each μ in order to determine the corresponding solution measure f 0 . The two approaches lead to the same set of steady states, although their arrangement into solutions branches depends on which parameter is used as the control parameter. This reflects the different stability properties obtained through the different procedures (see conclusion of [20]). Here, only the first approach results in the correct stability properties corresponding to the original conserved dynamics described by equations (4) and (7), respectively. In contrast, the second approach allows for perturbations which alter the mean concentration at fixed imposed chemical potential (e.g. via an external reservoir), a situation that does not correspond to conserved dynamics.

Liquid-liquid phase decomposition
3.1. Phase behavior First, we review the phase behavior of a binary mixture as described by the free energy (2). We minimize F[f] under the constraint of fixed total number of particles, i.e. at fixed mean concentration f 0 . This is equivalent to requiring the grand potential to be minimal. We consider passing to the TL, i.e. the limit of taking the system volume V=L d to infinity, where d is the dimension of the system, whilst proportionally increasing the particle number, so that f 0 remains constant as the limit is taken. In the TL the contribution to the free energy due to any interfaces can be neglected, since their contribution scales as L d−1 , and so the condition for a minimum becomes where ω h is the grand potential density. Moreover, the entropy density at fixed These results are illustrated in figures 1 and 2 using thick solid blue lines. For temperatures T below the critical temperature T c , equation (10) can have three solutions and depending on T and f 0 the system may phase-separate into regions with concentrations f=f + and f=f − . These coexisting concentration values are called the binodals and are given by a Maxwell (or double-tangent) construction on f. This construction results directly from the minimization of [ ] f F at fixed f 0 and volume V and implies that the chemical potentials μ=∂ f f and pressures p=−f+μf in the two phases f + and f − are equal, i.e. that the conditions for thermodynamic equilibrium hold. For the present symmetric potential (3), this results in the binodals . These results are illustrated in figures 1 and 2, using thick black lines. Note that the homogeneous state also exists in the concentration range between the binodals where it is not the thermodynamically (or globally) stable state. When crossing the binodal, the uniform state initially remains metastable, i.e. it represents a local minimum of the free energy, before becoming unstable on crossing the locus of ¶ = ff f 0. The corresponding f values are called the spinodals. For the present f they are . These metastable and unstable homogeneous states are illustrated in figures 1 and 2 using blue dashed and dotted lines, respectively, and the spinodal points are marked by filled circles.
Owing to the overall concentration constraint, coexisting states can only exist for f − f 0 f + . So, for any given f 0 , phase coexistence is only possible below . This condition defines the phase transition temperature T b (f) and implies (i) that in the 'critical' case (f 0 =0), T s , T b and T c all coincide and moreover that the order parameter ( ) (first order phase transition). The resulting representation of the phase transition in the phase plane spanned by f 0 and T and the respective continuous and discontinuous dependence of δf on T are shown in figure 3. Note that the jump in order parameter is not accompanied by a jump in ∂f/∂T (see figure 2(b)). This indicates that liquid-liquid phase decomposition is an example where the classical Ehrenfest definition of first order phase transition fails, while that based on a discontinuity in the order parameter holds. The corresponding entropy is 2 . This is as expected, as in demixing the system reduces its internal energy, paying a trade-off in entropy. As the system is not isolated, this process is accompanied by a flux of heat through the boundary.  (

Bifurcations in finite domains-phase separation in 1D
In the previous section we described in a compact manner the well-known thermodynamic behavior of a simple binary mixture close to the critical point for liquid-liquid phase decomposition. These results are valid in the TL, i.e. in the limit of diverging system size and particle number where the mean concentration is the relevant control parameter, in addition to the temperature. Implicitly it is also assumed that fluctuations have eliminated all transients, including metastable states. However, real systems have a finite domain size and finite observation times, and in such systems finite-size effects and transients may become important. Finite systems are often discussed in terms of bifurcation diagrams instead of phase diagrams, following dynamical systems theory and the theory of pattern formation [4,5]. Here we discuss the behavior of the Cahn-Hilliard equation for a finite domain size combining linear stability and bifurcation analyses, and focus on how the Maxwell construction emerges as the domain size is increased. For simplicity, we present first the case of a 1D domain before moving on to 2D.
In section 3.1 it is mentioned that a homogeneous state f=f h loses stability when ∂ ff f=0. The linear dynamic behavior may be studied by introducing the ansatz i.e. the growth or decay rate of a harmonic perturbation as function of its wave number q = |q|. Here h and κ are always positive. In contrast, can be positive or negative depending on the curvature of the local free energy density. Note that λ(q) is always real since equation (1) has gradient form. The system first becomes unstable when = 0 , corresponding to a long-wave instability. Above the instability threshold, there exists a band of unstable wave numbers 0<q<q + with the largest λ at = + q q 2 max . The zero crossing at q=q + corresponds to a steady state bifurcation from the homogeneous state f=f h .
Fixing the domain size to L selects the wave number p = q n L 2 L n , where n=1, 2, ..., and the zero crossing of (13) determines the threshold value of f L n (at fixed T) or of T L n (at fixed f). Here, we focus on the first option and obtain the thresholds . Panel (a) shows the liquid-liquid demixing phase diagram in the plane spanned by the mean concentration f 0 and the temperature T. The binodals and spinodals are shown using thick black and thin dashed red lines, respectively. The one-and twophase regions lie above and below the binodal, respectively. The vertical solid, dashed and dotted lines indicate the range of thermodynamically stable, metastable and unstable homogeneous states for a specific choice of f h =f 0 . The limiting temperatures for this value of f h are indicated by red symbols. Panels (b) and (c) show the dependence of the order parameter δf on T in the case of a discontinuous transition (f 0 =0.4, first order) and a continuous transition (f 0 =0, second order), respectively, in the thermodynamic limit. In (b) the states on the binodal that cannot be realized for the chosen f 0 are indicated by a black dotted line.
all converge to f s . Numerical continuation allows us to determine the emerging branches of heterogeneous steady states. The most relevant is that with n=1 whose chemical potential and mean free energy are shown for several values of L in figure 4, where the TL is also included.
We now discuss how the character of the bifurcation curves (with f 0 as the primary control parameter) changes when the secondary control parameter consisting of the domain size L is increased: at the smallest L that allows a heterogeneous state to develop (L min =2πℓ where ℓ ) a pair of supercritical symmetry-breaking pitchfork bifurcations appears in a codimension-2 bifurcation at f 0 =0. With increasing L these primary bifurcations move apart, generating in the f 0 range between them an ever longer branch of heterogeneous states (see e.g. the L/ℓ=7 case in figure 4) that represent the lowest energy state at these f 0 (see figure 4(b)). At ℓ p = = L L 10 ss both pitchfork bifurcations become subcritical 8 , i.e. the branch of heterogeneous states emerges towards the linearly stable homogeneous states f=f 0 . Close to the bifurcation, these heterogeneous states are then linearly unstable and of higher energy than the homogeneous state at identical f 0 , as is clearly visible for curves with L/ℓ>10 in figure 4. These unstable states correspond to threshold or nucleation solutions that have to be overcome in order to jump between the linearly stable homogeneous state and the linearly stable heterogeneous states that exist beyond the saddle-node bifurcation at f=±f sn where the branch of heterogeneous states turns around and acquires linear stability. Shortly after turning, the heterogeneous state becomes the global free energy minimum. Typical examples of concentration profiles can be found in the literature; see e.g. figure 3 of [52] and figure 2 of [53].
On further increasing L, the primary bifurcation points move further away from f=0, ultimately converging on the spinodals f s as  ¥ L . At the same time the saddle-node bifurcations move outward and converge on the binodals f b . The branch of unstable heterogeneous states becomes longer and ultimately approaches the states represented by a dashed line that correspond to the metastable states in the TL. These states represent a branch of critical nuclei for phase separation. Finally, the branch of stable heterogeneous states between the saddle-nodes at ±f sn becomes increasingly straight and horizontal and converges to the Maxwell line for  ¥ L . We emphasize, however, that even in the TL there is no bifurcation between the homogeneous and the heterogeneous state at the binodal and that the branch of unstable nuclei is an intrinsic part of the overall picture.
The overall manner in which the Maxwell construction emerges from the bifurcation scenario in the limit of ever larger domain size is not influenced by the additional branches emerging at f L n with n>1, since these always correspond to states of higher energy than the n=1 branch and never connect to it. This is because for large L the states on the n=1 branch in a periodic domain consist of a single region where f(x)≈f b together with a second region where f(x)≈−f b , with two interfaces between them. For n>1 the stationary periodic states consist of a larger number of single phase regions with a correspondingly larger number of interfaces. Note that although the n>1 states never appear as global free energy minima, they may appear as transients in nonequilibrium time simulations. However, this 'simplest possible' picture does not hold in higher dimensions as discussed below.
In figure 5(a) we display the value of the order parameter for various states with f 0 =0.4 as the temperature is varied. This is the same as figure 3(b), but now with the heterogeneous states for various values of L also included. We see that as  ¥ L these tend to the curves displayed in figure 3(b). However, we also have a branch of unstable states connecting the two stable branches.
Similarly, one can construct a finite-size equivalent of figure 3(a)-the phase diagram in the TL. First, we need to define what corresponds to the finite-size equivalents of the spinodal and binodal lines. The equivalent of the spinodal corresponds to the linear instability threshold of a homogeneous state of a finite system and is given by equation (14) with n=1. However, for the binodal lines the question is more involved as these correspond to the two coexisting states in the TL which is unaffected by interfaces. Since interfaces are an intrinsic part of (inhomogeneous) states in finite-size systems, the finite-size equivalent of the binodal needs to be defined as a property of these states. We use the concentration values where the inhomogeneous state becomes the global free energy minimum. Figure 5(b) gives the result for two values of the system size. A notable feature is that for finite systems there exists a finite concentration range where the transition is 'second order', together with the corresponding 'tricritical points' where the transition becomes 'first order'. This concentration range increases with decreasing domain size.

Bifurcations in finite domains-demixing in 2D
We now discuss the emergence of the TL for a demixing system in 2D. In contrast to the 1D system discussed in section 3.2, we must now include additional periodic structures in the discussion. As a result, the corresponding bifurcation diagram is richer and the transition to the TL that takes place as the domain size  ¥ L exhibits more interesting features. The inhomogeneous steady states that can arise vary considerably and depend on the domain size, but the most typical ones consist of a circular drop of one phase surrounded by the other or, when f 0 is closer to zero, of slabs (stripes) connected via the periodic boundary conditions. Figure 6 shows the chemical potential μ and rescaled mean free energy density f F L a 2 ref 2 for cluster (drop or hole) states as well as for stripe states as a function of the mean concentration f 0 . Square domains of size L×L are employed and results for two different values of L are given. The stripe states (green lines in figure 6) correspond to the 1D states discussed in section 3.2 that are now extended into the second dimension in a translation-invariant manner. As a result, they have identical bifurcation curves and identical behavior as  ¥ L as the 1D states. However, in 2D new perturbation modes are available and therefore their linear stability properties differ from those in 1D. In particular, the stripe states may now be unstable to transverse perturbations (equivalent to the Plateau-Rayleigh instability of liquid ridges on a solid substrate [54,55]). Here, this implies that the stripes no longer stabilize at the saddle-node bifurcation at f sn stripe where the subcritical part of the branch turns around, as in the 1D case. Instead they become linearly stable at a subsequent secondary pitchfork bifurcation where an unstable branch of transversally modulated stripes (magenta dashed lines in figure 6) emerges subcritically from the branch of translation-invariant stripes.
Additionally, on the square domain considered here, one finds cluster states, consisting of round clusters of phase 2 in phase 1 ('drops', generally for f 0 <0) and similar clusters of phase 1 in phase 2 ('holes', generally for f 0 >0). These states form a single curve in each bifurcation diagram (red lines in figure 6) and extend between the same primary bifurcation points as the stripe states. As for the stripes, at small L the branch emerges supercritically (not shown) while at larger L the primary bifurcation is subcritical (see, e.g. the case of L=4L c in figures 6(a) and (b)). The branch stabilizes at a subsequent saddle-node bifurcation (at f  sn 1 cluster ) at which it turns around, before losing stability again at a secondary pitchfork bifurcation where the branch of modulated stripes terminates. This shows that the branch of modulated stripes is related to the exchange of stabilities between the stripe and the cluster branch, namely, it corresponds to critical saddle states that have to be overcome in order to transition between these two linearly stable states. On further increasing L, the cluster state branch ceases to be monotonic between the two saddle-node bifurcations at f  sn 1 cluster and undergoes a hysteresis bifurcation, at With increasing domain size L the saddle-node bifurcations at f  sn 1 cluster on the cluster branch and the corresponding bifurcations at f  sn stripe on the stripe branch gradually move outwards and approach the branch of homogeneous states at the Maxwell point. The saddle-node bifurcations at f  sn 2 cluster on the cluster branch likewise move outwards towards larger | | f 0 . As  ¥ L they approach , a prediction that follows from the sharp interface limit 9 . During this process the two plateaus in μ become longer and their difference in chemical potential, m m D º 2 sn 2 cluster , decreases as L −1 . The convergence of f  sn stripe , f  sn 1 cluster and f  sn 2 cluster as well as the way all μ sn approach zero with increasing L is illustrated in the log-log plots of figure 7. There we see that m sn 1 cluster decreases as L −2/3 .
Note that the stripe state (parallel slabs of the two coexisting phases) is the state of lowest free energy state in a region centered on f 0 =0, while for larger | | f 0 the cluster state has the lowest free energy (see figures 6(b) and (d)).
The sharp interface limit indicates that in the limit  ¥ L the stripe state [cluster state] represents the minimum energy state provided | .
The values of f f -¥ sn sn also decrease as power laws as  ¥ L , apparently in the same manner as μ sn , as revealed in figure 7(b). In this section we have seen that close to the primary bifurcation the emerging curves of 2D inhomogeneous states behave much like in the 1D case. However, in the region around f 0 =0 the behavior differs strongly as more states are accessible to the system in 2D, resulting in a different convergence pathway to the Maxwell line: several μ plateaus develop that converge to the single horizontal line of the Maxwell construction only in the TL. However, different regions on this line continue to represent different states.

Phase behavior
As in section 3.1 for liquid-liquid demixing, we first review the crystallization phase behavior described by the free energy (6). Minimizing F[f] at fixed mean density f 0 gives where μ represents the chemical potential that determines the mean value of the order parameter f 0 , here defined as , where V is the volume of the domain. The liquid state is homogeneous and the corresponding spatially uniform solution of equation (15) figure 8). In the r range directly below the critical point, the crystallization phase transition is predicted by the PFC model to be of second order. The transition is of first order only below the 1D tricritical point situated at ( tri tri [20]. Since freezing is in reality a first order transition [59], the PFC model is not correct for r>r tri . This is a consequence of the approximations made in deriving the model [38]. At (r tri , f tri ) the binodals emerge (black solid lines in figure 8) that limit the coexistence region and specify the densities f b li and f b cr of the liquid and crystalline states that coexist at a given undercooling r, i.e. the states with the same chemical potential and pressure (i.e. the same grand potential density). The binodals can either be calculated for particular spatial periods or wavelengths of the crystalline structure or for an infinite domain. In the latter case, the wavelength of the crystalline state is that which minimizes the free energy. The binodals in figure 8 are calculated at fixed structural length corresponding to the critical wave number L c =2π/q. Note that in the following we always use q=1, i.e. L c =2π. The results of the two approaches cannot in general be L (see legend) and describe how the three overlapping μ plateaus corresponding to drop, hole and stripe states approach one another to form a single Maxwell line in the limit. 10 In the sharp interface limit, the interface energy of a stripe state in a L×L domain is ∼2L, independently of f 0 . For the circular cluster state the interface energy is ∼2πR where the radius R depends on f 0 . Averaging the density for a fully decomposed state with f=±f b one , implying that the interface energy of the cluster state is lower than that of the stripe state when distinguished on scales such as that used in figure 8. As for the Cahn-Hilliard equation, the binodal lines are determined using numerical path continuation as done for the 1D case in [20]: first, we identify the coexisting values of the densities f b li and f b cr at some particular undercooling r. We then continue these values as a function of r. Figure 9 shows the chemical potential μ and the mean free energy density F/L of the 1D liquid and the spacefilling crystalline states at fixed r as a function of f 0 using, respectively, blue and black lines 11 . Equation (15) also possesses solutions corresponding to finite size portions of crystal that coexist with a liquid background. These crystallite 'localized states' are present in the coexistence region (and slightly outside) and are discussed in greater detail in [20] (see their figures 3 and 5 for typical solution profiles). These states exhibit the phenomenon of slanted snaking [13,61,62] but take the form of standard homoclinic snaking when the order parameter f 0 is plotted as a function of the chemical potential, a property that is expected to carry over to other models with a conserved order parameter. See the conclusion of [20,26] for further discussion. Note that figure 9 highlights the states representing global energy minima using thick line segments, an approach also taken for other systems showing multiple steady states, as done, e.g. in [14] for various buckling states.
In 2D the phase diagram of the PFC system is much richer (see figure 10 of [20]). In addition to the uniform or liquid state, the PFC model now exhibits stripe-like states, and two distinct periodic states with hexagonal coordination, one of which consists of density maxima on a hexagonal lattice (referred to as bumps, i.e. the crystal state) and the other with similarly arranged density minima (referred to as holes). As a result, the system  may exhibit phase coexistence between the uniform state and the bump state, between the bump state and stripes, between stripes and holes and between holes and the uniform state. Here, we only consider the first scenario arising from thermodynamic coexistence between the liquid (homogeneous) state and the crystalline (bump) state, and the spatially localized structures associated with this coexistence. Figure 8(b) shows the corresponding part of the phase diagram, with the spinodal line indicating the onset of linear stability of the liquid state and the binodals specifying the coexisting liquid and hexagonal crystal states at a given r. As in 1D, the coexistence region is associated with the presence of 2D localized crystallites.
In the next two sections we show that the localized states present in both 1D and 2D are intimately related to the emergence of the Maxwell construction in the TL for the liquid to crystal phase transition. The results represent a third way in which the Maxwell construction is approached as  ¥ L . In both cases, 1D and 2D, we use the values of the chemical potential at coexistence as reference values.

Bifurcations in finite domains: crystallization in 1D
The properties of the liquid (homogeneous) and crystalline (periodic) bump states and of the localized structures associated with their coexistence are summarized in figure 9 in terms of a bifurcation diagram for a relatively small domain of L=16L c . The figure shows (a) the chemical potential μ and (b) the relative mean free energy density (F-F 0 )/L where F 0 =F[f 0 ] as a function of f 0 , for r=−0.9. Thin solid and dashed lines are employed for linearly stable and unstable states, respectively. The figure also indicates, using thick solid lines, the thermodynamically stable state for each value of the mean concentration f 0 , i.e. the global free energy minimum. These states are identified in figure 9(b). For alternative representations using the norm or grand potential as order parameters and typical solution profiles, see [20]. Figure 9 shows that the periodic (or domain-filling crystalline) state (black line) bifurcates in the direction of decreasing stability of the uniform state f=f 0 (blue line), i.e. supercritically. This crystalline state loses stability at small amplitude to a pair of spatially localized structures (red and green lines) resembling crystalline states of finite length embedded in the background uniform liquid state. These states are distinguished by their behavior at x=0, i.e. at the center of the pattern, with the states in red having a maximum at x=0 and those in blue having a minimum at x=0, and both exhibit snaking. This behavior in turn implies the presence of interconnecting branches (resembling rungs on a ladder) of unstable asymmetric structures that are computed in [20] but not shown here. The figure indicates the linear stability properties of each state shown, with thin solid lines indicating linearly stable states (or local free energy minima) and thin dashed lines indicating unstable states. Note that the localized states coexist with the stable uniform state, a possibility that only arises because of mass conservation [63], and that they represent the global free energy minimum over a large part of their range of existence ( figure 9(b)). The inset of figure 9(b) indicates that the transition between successive global minima occurs via swallow-tail-like structures. Figure 9(a) shows that the localized states corresponding to global free energy minima are clustered within a band of width Δμ≈0.031 08 within the snaking region, until superseded by the periodic state at large f 0 . Figure 10(a) displays only these global free energy minima for various domain sizes L. Connecting adjacent branches of such minima generates a zigzag curve, and figure 10(a) shows such zigzag curves for eight different domain sizes ranging from L=6L c to 512L c . We see that the slopes of the slanted segments of the curves are all the same and independent of the domain size L. However, the length of these segments decreases with increasing L leading to a corresponding increase in the number N(L) of slanted segments. Figure 10(b) shows that the resulting stability interval Δμ follows, with excellent accuracy, the power law Δμ∼aL b with exponent b∼−1. This exponent is consistent with the observation that the slope of the slanted portions of the curve is independent of L and moreover that the number of oscillations in each snaking curve is proportional to L. The latter is in turn a consequence of the fact that each oscillation in the curve results in the addition of the same wavelength on either side of the localized state as the structure grows. We conclude that in the limit  ¥ L the stability interval shrinks to zero and conjecture that in the limit the resulting curve is nowhere differentiable. We also confirm that the location of the limiting chemical potential corresponds precisely to the Maxwell construction for these parameter values. Thus, from this point of view the Maxwell construction involves a continuum of localized structures, all coexisting at the same value of the chemical potential μ.

Bifurcations in finite domains: crystallization in 2D
We next consider the case of two-dimensional domains focusing on the transition between a liquid state and a hexagonal crystalline state. For this purpose we use hexagonal domains with specific side lengths as physical domains . As explained in section 2.3, numerical continuation is then applied on a domain  num that corresponds to one twelfth of  with appropriate boundary conditions reflecting the symmetries of the states considered. Employing other domain shapes and boundary conditions corresponding to different symmetries can affect aspects of the results. For instance, a 'wrong' domain shape can rule out the existence of the periodic hexagonal state and bifurcations may become imperfect. However, if the domain is sufficiently large one always finds the snaking branches of hexagonal patches discussed below as long as these are small compared to the domain. The only part of the bifurcation structure that is affected by the domain corresponds to the transition from the crystalline patch state to a periodic crystalline state, and reflects the interaction of the patches with the boundary. For hexagonal domains of different side lengths L h , figure 11 displays the corresponding bifurcation diagrams of steady homogeneous, periodic and localized states for r=−0.9. Examples of localized states with hexagonal coordination are presented in figure 12. At first sight, the diagrams in figure 11, which are in terms of the chemical potential and the free energy density, have a similar appearance to those for the 1D case displayed in figure 9. However, they differ in several significant aspects: following the stable homogeneous state in the direction of increasing density, we see that this time it becomes unstable in a subcritical pitchfork bifurcation where a solution branch of defect-free domainfilling hexagonal patterns emerges (black dashed curve). A branch of nearly rotationally invariant localized target-like solutions (dark green dashed curve) emerges from this branch at very small amplitude, even though the domain does not have this symmetry (see figure 12(a)). The target-like structure grows in radius along the branch, although only part of it is shown in figure 11. Although nominally axisymmetric the target structures do in fact reflect the symmetry of the domain which is responsible for the presence of a very slight D 6 deformation of this state.
The branch of localized hexagonal patches (red curve) bifurcates in a tertiary steady-state bifurcation from the target pattern branch. The bifurcation is a D 6 -equivariant pitchfork bifurcation and so breaks the (nominal) axisymmetry of the target state. However, the above-mentioned slight D 6 deformation of the target structures makes the bifurcation slightly imperfect as can be ascertained from continuation runs with a smaller stepsize.
The states presented in figure 12 (with the exception of the first profile) all lie on this latter emerging branch and show that as one follows the branch the hexagonal patch gradually increases in size. This process occurs via the addition of new bumps at preferred locations along each side and is followed by the addition of further bumps on either side until each row is complete. The first two profiles in the second row show some of the intermediate states on the branch segment extending from a 4 bumps-on-a-side patch to a 5 bumps-on-a-side patch. Hexagonal patches in the nonconserved Swift-Hohenberg equation grow in the same manner [64]. The appearance of each new cell is associated with a fold in the snaking branch, i.e. as the patch increases in size the density of the branches in the bifurcation diagram increases. and show behavior that is, in principle, similar to that in figure 10(b). However, the resulting structure is more rugged and irregular, and the width Δμ has to be carefully defined as the individual inclined straight segments are not all centered about the coexistence value. We therefore determine the maximum and minimum μ value for each inclined segment and average them to obtainm max andm min , thereby obtaining¯m m m D º max min . Figure 14(b) shows the resulting widths of the region of thermodynamically stable patches Δμ as a function of   1 . Therefore, for sufficiently large domains one expects to find the power law m D~-L h 2 . However, in our calculations for the hexagonal patches we are only able to go to a maximum side length of L 32 c h . Despite this limitation the observed dependence of Δμ on L is well described by the full expression derive above, as illustrated in figure 14(b).
To strengthen our argument that the different scaling behaviors can be attributed to the fully 2D nature of the hexagonal patches, we consider stripe-like arrangements of hexagonally ordered bumps. See [64,65] for an analysis of similar states in the case of the standard nonconserved Swift-Hohenberg equation. We find that these exhibit scaling behavior that is identical to that of truly 1D structures. Figures 15 and 16 show the corresponding (f 0 , μ) bifurcation diagrams, the associated zigzag curves and the scaling of Δμ.  figure 15(d)).
These results show that quasi-1D localized structures behave in a manner that is very similar to strictly 1D structures; the approach to the Maxwell construction in the limit  ¥ L x is therefore also similar although, naturally, the coexistence value of the chemical potential that is approached corresponds to the 2D Maxwell line for hexagonal crystal-liquid coexistence and not the 1D Maxwell line for stripe-liquid coexistence. In particular, we expect that in this case the resulting Maxwell curve consists of identical linear segments of vanishing length   and is likewise nowhere differentiable. However, the case of fully 2D hexagonal structures is different ( figure 13), since there the zigzag segments have different lengths, although all appear once again to have the same slope. This distribution of lengths reflects the different effects on the chemical potential of adding bumps in the middle of a face (if the face has an even number of bumps) or two bumps on either side of the middle (if it has an odd number of bumps) and then of adding further bumps at the remaining sites along each face. These differences from the quasi-1D case imply that the zigzag structure does not repeat as one traverses the zigzag curve from one end to the other, a fact that is responsible for the departure of the exponent b from its 1D value b=−1.
Note, finally, that the inset in figure 16 reveals that the zigzag curves for L 64 x c and L 128 x c are both slightly off the exact coexistence value of μ indicated by the red horizontal line. We believe that this results from the resolution of the spatial grid that cannot be sufficiently refined for full convergence. This effect is also visible in the curve for = L L 16 h c h in figure 13. Despite this, Δμ still follows a clear power law ( figure 14(a)).

Conclusions
In this paper we have revisited the Maxwell construction that predicts the location of a first order phase transition in thermodynamic systems at finite temperature. Our work sheds new light on the process whereby a finite size system approaches the TL. We have considered two basic mean field models, a phase separation model described by the Cahn-Hilliard equation and the PFC model or conserved Swift-Hohenberg equation that describes the process of crystallization from a melt. The former case is simpler since the two phases involved are both uniform. As a result any departure from the TL arises from the presence of interfaces whose contribution to the free energy vanishes in this limit. Despite this well-known property our analysis sheds new light on the TL: for example, in 2D it shows that the Maxwell construction (which we recover) involves different states (stripe and cluster states) depending on the mean concentration f 0 (see figure 6). In particular, in finite size domains, however large, the minimum energy state does not have the same spatial structure for all values of f 0 .
The second model, the PFC model, reveals additional complexity. In 1D, finite systems are characterized by a large number of spatially localized states lying on a pair of intertwined branches that straddle the Maxwell point between the homogeneous (liquid) phase and the periodic (crystalline) phase. As in nonconserved systems these states gain and lose linear stability at successive folds along these snaking branches, and in a f 0 (μ) plot these folds are located at particular values μ ± of the chemical potential (m m m < < -+ Maxwell ), i.e. the folds are aligned [20]. This behavior is independent of the domain length L once L is large enough, provided the localized structures remain localized, and do not fill the whole domain ( figure 9). However, what does depend on L is the number of available states which increases in proportion to L generating more and more oscillations across the Maxwell point. In addition, the interval in μ, Δμ, within which these states correspond to the global free energy minimum decreases to zero, to good accuracy, as L −1 . This is because the minimum energy state jumps discontinuously to longer and longer localized structures as f 0 increases and each of these states remains a global minimum for the same μ interval Δμ. In the limit  ¥ L one therefore finds that the Maxwell line consists of a dense set of transitions between different and successively longer localized states. Thus in this case, too, the Maxwell line corresponds to a succession of distinct states (figure 10).
In 2D the PFC model is even more interesting since the cluster states now typically consist of hexagonal patches. In a f 0 (μ) plot these patches snake in the same way as in nonconserved systems [20], implying that they grow by adding a bump in the middle of each interface followed by the addition of bumps on either side of these first bumps until a new and larger patch has been created. Because of this gradual increase in area the growth process does not repeat in a periodic manner, with the number of back and forth oscillations across the Maxwell line increasing with f 0 . Here too the interval Δμ within which each patch corresponds to minimum free energy is well defined and it too decreases to zero as the domain area grows. For hexagonal domains we find that m D~-L h 2 , where L h is the domain scale ( figure 14). As in the 1D case in the limit  ¥ L we recover the Maxwell construction, with the Maxwell line consisting of a dense set of transitions between distinct hexagonal patches ( figure 13). We have argued that the different scaling exponent, b∼−2, is a consequence of the 2D nature of the patches since a similar analysis of hexagonal stripes, i.e. localized regions of hexagons with only one extended direction, exhibit the same Δμ scaling as the strictly 1D problem.
Both of the mean field models we have considered can be derived as local (gradient expansion) approximations of density functional theory (DFT) [19,38,66,67]. Formally, DFT is derived by averaging over all states of the system, i.e. including all thermal fluctuation effects. However, since in practice one must make an approximation for the free energy functional, the mean field models do not include the full contribution of thermal fluctuations. However, since at least some fluctuation contributions are already included and averaged over, one should not add additional fluctuating terms. Instead, to obtain a quantitatively more accurate theory, one should employ better approximations for the free energy functional. This would result in quantitatively different results, but the qualitative behavior regarding the emergence of the Maxwell construction would likely remain unchanged. With these mean field models complete bifurcation diagrams can be determined including the unstable and metastable states that computer simulations are normally unable to capture. To see how the full spectrum of fluctuations amends the picture one would need to pair our approach with large scale computer simulations, as done in the context of partially wetting droplets on solid substrates in [68].
Finally, we discuss another limitation of our study. We have considered finite size systems with ideal boundaries, i.e. periodic or Neumann boundaries (for an extensive discussion of their relation see [69]). The latter do not influence the phase behavior at the boundaries, but only break the translation invariance of an infinite system. An interesting question for future work concerns the role played by rigid boundaries in the approach to the TL. In the case of mixtures these may represent preferential adsorption of one component or a changed interaction between components at a boundary [70]. Such boundaries may result in a very rich surface phase behavior (see e.g. figure 6 in [71] for a phase diagram, figures 3-14 in [53] for bifurcation diagrams and 1D concentration profiles, and [58] for 2D results). Incorporating such boundaries into the present study will most likely add a new level of complexity to some aspects of the bifurcation diagrams, allowing one to investigate the interplay between different bulk (liquid-gas, demixing, crystallization) and surface (wetting, pre-wetting, surface freezing, pre-melting) phase transitions in finite systems and the corresponding transition towards the TL.