Shape and interfacial structure of droplets. Exact results and simulations

We consider the fluctuating interface of a droplet pinned on a flat wall. For such a system we compare results obtained within the exact field theory of phase separation in two dimensions and Monte Carlo (MC) simulations for the Ising model. The interface separating coexisting phases splits and hosts drops whose effect is to produce subleading corrections to the order parameter profile and correlation functions. In this paper we establish the first direct measurement of interface structure effects by means of high-performance MC simulations which successfully confirm the field-theoretical results. Simulations are found to corroborate the theoretical predictions for interface structure effects whose analytical expression has recently been obtained. It is thanks to these higher-order corrections that we are able to correctly resettle a longstanding discrepancy between simulations and theory for the order parameter profile. In addition, our results clearly establish the long-ranged decay of interfacial correlations in the direction parallel to the interface and their spatial confinement within the interfacial region also in the presence of a wall from which the interface is entropically repelled.


Introduction
The characterization of the interfacial region separating coexisting phases is a classical problem in statistical physics since the pioneering works of van der Waals [1]. In particular, it has been shown by van der Waals that the interface separating a liquid at coexistence with its vapor phase can be described in terms of a density profile that smoothly interpolates between the bulk values attained far away from the interface [1,2]. More generally, the characterization of fluctuating interfaces is a classical problem in statistical physics, and near to critical points the features of universal behavior also emerge for extended fluctuating objects [3,4]; we refer to [3,[5][6][7][8][9][10][11][12][13] for reviews on interfacial phenomena.
In statistical mechanics, the interface separating phases a and b can be imposed by means of boundary conditions (see e.g. [14] for rigorous treatments). Within a coarsegrained view the interface separating a and b is regarded as a sharp entity whose fluctuations lead to a smooth density profile that interpolates between pure phases. By going beyond the picture of a structureless interface, it is legitimate to expect that interfaces are not sharp but rather they possess some kind of structure. Intuitively, in two dimensions the interface structure can be understood in terms of processes in which the line separating coexisting phases self-intersects or branches by hosting droplets of nonboundary phases c ̸ = a, b, as depicted in figure 1. From simple topological considerations it is clear how branching is possible for the q-state Potts model with q = 3 and q = 4 while for the Ising model the interface must necessarily spilt in three lines hosting a sequence of four phases with the pattern −| + | − |+ [15][16][17][18]. The notion of interface structure is crucial for the understanding of interesting phenomena such as the formation of droplets adsorption along the interface and the spreading responsible for the occurrence of interfacial wetting [6,9,11,12].
The two-dimensional case turns out to be interesting because several results are available from the scaling limit of exact calculations on the square lattice [17]. Such findings proved to be fundamental for establishing the connection between random walks and fluctuating interfaces at equilibrium [4]. Analytical progress on interfacial phenomena in two dimensions has been obtained within the exact field theory of phase separation developed in [18][19][20][21][22][23]. Results for many-body correlation functions have been obtained in [24] and their test in Monte Carlo (MC) simulations has been performed in [25,26]. In particular, interface structure effects have been examined in [25] for the two-dimensional Ising model on the strip. Then, the occurrence of interfacial adsorption and interfacial wetting in q-state Potts models on the strip predicted in [20] has been tested in [27].
From the theoretical perspective, the case of the Ising model on the half-plane has fundamental importance since it provides the first exact solution of a wetting transition in two dimensions [28]. In more recent times it has been shown how it is possible to go beyond the analysis of wetting in the Ising model by means of field-theoretical techniques. In particular, it is has been possible to obtain exact results for order parameter profiles in other universality classes 3 [19] as well as criteria for wedge filling [21,22].
For the half-plane geometry early MC studies on Ising droplets pointed out that finite-size effects are indeed quantitatively relevant [31]. As shown in [31], the droplet profile without finite size corrections-which has derived long time ago [19,28,[32][33][34]exhibits deviations from MC data. This is something that we are going to show with the aid of high-performance MC calculations, meaning that the origin of the discrepancy is not completely due to the early simulations. We will find that by taking into account the interface structure the agreement between theory and simulations is perfectly achieved. The analytic expression of the interface structure correction for the droplet profile has been recently calculated, we refer the interested reader to [35] for further details.
To the best of our knowledge 4 , the simulation studies of [31] provided the first MC results for interfacial correlations but unfortunately the comparison with the theory was not possible because the analytical results for the pair correlation function within the interfacial region was not available at that time. This theoretical gap has been recently filled and only now it is possible to compare the analytical results with the MC ones-which we also perform. In addition to the leading order result, also the interface structure correction has been computed in closed form and compared against simulations in this paper.
The scope of this paper is thus to show how interface structure plays an important role in the study of density correlations in the interfacial region. To do this, we run high-precision MC simulations, improving on the early MC studies. Then, we compare the outcome of our MC results with the analytical predictions calculated within the exact theory of phase separation in [35]. The addition of the recently calculated highorder correction to the density profile is able to yield significant agreement between theory and simulations. Beside the remarkable agreement between theory and numerics that we will present, one of the most significative results of this paper is to demonstrate through the specific example of Ising interfaces how interface structure and the presence of a wall reflects in the long-range character of interfacial correlations.
This paper is organized as follows: first we outline the analytical results which are the object of a separate work. Then, we illustrate the comparison between theory and simulations for the density profile; lastly a section on two-point correlations follows. Conclusions and perspectives are drawn in the last section.

Summary of analytical results
In this section we recapitulate the exact analytic expressions for one-and two-point correlation functions of the order parameter in the presence of a fluctuating droplet. The detailed calculations involved in the derivation of the results presented in this section are provided in [35]. Although the theory encompasses several universality classes within a unified framework, we specialize the analytical results to the Ising model and we test them by means of MC simulations. The system we consider is therefore the two-dimensional Ising model along the coexistence line T < T c and H = 0, where T c is the critical temperature and H is the bulk field. We briefly recall the notations and conventions. In our simulations we consider the Hamiltonian where s i ∈ {±1} for Ising spins and the sum is restricted to nearest neighboring sites of a two-dimensional square lattice. The critical temperature is T c /J = 2/ ln(1 + √ 2) ≈ 2.269 [37] and the bulk correlation length is given by [17], with the reduced lattice coupling K = J/T , and its dual K ⋆ defined by exp(−2K ⋆ ) = tanh K. The spontaneous magnetization, M = ⟨σ⟩ + = −⟨σ⟩ − > 0, is given by Yang's formula, M = ( 1 − (sinh(2K)) −4 ) 1/8 , and in the scaling limit it reduces to M ∼ (T c − T ) 1/8 [38,39]. Green circles indicate the order parameter fields σ 1 and σ 2 appearing in correlation functions considered in this paper. For the Ising model we take b = +1 and a = −1.
Throughout this paper we consider a ferromagnetic coupling J > 0 and without loss of generality we set J = 1 in our simulations.
The geometry we are interested in is the one of the half plane with boundary conditions on the wall x = 0 enforcing the formation of a droplet in the semi-space x > 0. For |y| > R/2 spins on the wall take the value σ i = +1 while for |y| < R/2 the reversed patch of boundary spins is such that σ i = −1, where i is a lattice point on the wall; see figure 2. This protocol defines the portion of wall wetted by a droplet of negatively magnetized phase. The droplet such defined fluctuates in the bulk but its extremities, the pinning points (0, ±R/2), are not allowed to fluctuate. It is implicit that we are dealing with a system in which the volume enclosed by the droplet is not conserved. In general, imposing constraints such as fixing the total mass of the profile leads to interesting effects on free energies [40] and on correlations in interfacial phenomena [41].
In the continuum limit the discrete lattice spin s i is replaced by a field σ(x, y) and the lattice node i is replaced by the coordinates (x, y) in the Euclidean plane. From now on we focus on statistical averages of the order parameter field on the above mentioned geometry, something that we denote ⟨σ(x, y)⟩ B +−+ and ⟨σ(x 1 , y 1 )σ(x 2 , y 2 )⟩ B +−+ , respectively for one-and two-point correlation functions. The order parameter profile is given by [19,35] ⟨σ(x, y) where P 1 (x, y)dx is the probability that the interface intersects the interval (x, x + dx) at ordinate y, its expression is provided on due course. The parameter m, which is the surface tension (in k B T units) of the interface separating the coexisting +/−phases, is related to the subcritical correlation length via ξ b = 1/(2m). The latter is is an exact form of Widom's relation valid for all subcritical temperatures that follows from duality [17,42,43]; see also [44,45]. The coordinates x and y enter in the above expression through the dimensionless combinations 3) The first term on the right hand side of (2.2) is where erf(χ) is the error function [46]. The order parameter profile is a smooth function that interpolates between −1 (for χ = 0) and +1 as χ approaches +∞. The scaling function Υ(χ) already appeared in the literature on Ising droplets (see e.g. [19,28,[32][33][34]) while the correction due to interface structure is a new result, the latter is the term ∝ P 1 [35]. It has been shown in [19] that the scaling function Υ(χ) is not a specificity arising from the Ising model but rather it is universal in the sense that it appears for other universality classes admitting droplet shapes (e.g. the q-state Potts model with q ⩽ 4 [47,48]). Note that for x → +∞ the magnetization profile approaches the value M corresponding to a pure phase obtained from a uniform wall with uniform boundary condition + along it. This result is consistent with the fact that for any finite R interfacial fluctuations are of order R 1/2 and at large x ≫ R 1/2 the droplet does not affect the value of the order parameter. The exception to such a behavior is provided by the wetting transition (see e.g. [19]). The passage probability is with the normalization´∞ 0 dx P 1 (x, y) = 1, and χ is given in (2.4) [19]. Since the second term on the right hand side of (2.2) is proportional to R −1/2 , it follows that it is a finite-size effect. In general, these interface structure effects in the magnetization profile along the x -axis can be organized in the form of a power series whose generic term is of the form R −n/2 U n (η), with n = 1, 2, . . . and U n (η) is a universal scaling function of η = x/λ. Again, the prefactor of U n (η) is not universal. The second term in the right hand side of (2.2) corresponds to n = 1 and the subsequent correction hidden in O(R −1 ) comes from n = 2, and so on. These corrections can be computed systematically within the field theoretical approach, as outlined in [35]. The scope of this paper is to show how the interface structure correction ∝ P 1 (x, y) is essential in order to achieve a quantitative agreement between theory and MC simulations. As we are going to show, in order to make a comparison between theory and simulations it is necessary to take into account the interface structure correction term ∝ R −1/2 . Of course, for sufficiently large R it is -in principle-possible to observe agreement between simulations and theory at the leading order but -de facto-even for the largest value of R we have simulated the effect due to the term ∝ R −1/2 is non negligible.
There is another point to emphasize. Strictly speaking the analytical result (2.2) is valid for x ≫ 1/m. Technically, this restriction comes from the fact that the spin field is treated as a bulk field and this is possible only far away from the wall. In fact, it is known from the scaling limit of the exact solution of the Ising model on square lattice that the wall affects the magnetization in a layer of thickness of order of the bulk correlation length, which is ∝ 1/m [49]. These deviations will be reported in our simulations.
Let us consider now the pair correlation function of the order parameter. In this paper we focus on two particularly symmetric arrangements of spin fields that allow for mathematically simple expressions of the analytic results. We thus focus on the so-called parallel correlation function (∥), denoted ⟨σ(x, y)σ(x, −y)⟩ B +−+ , in which spin fields lie parallel to the wall and moreover are spatially equidistant from the symmetry axis of the system, the x -axis. The second type of two-point function we will examine is the perpendicular correlation function (⊥), ⟨σ(x, 0)σ(x + d, 0)⟩ B +−+ ; in this case both spin fields lie on the x -axis at distance x, and x + d from the wall, respectively. The arrangement of spin fields defining the correlation functions given in (2.6) is illustrated in figure 3.
The aforementioned correlation functions are conveniently written in terms of rescaled variables by defining with η = x/λ and τ = 2y/R. The analytic expression for the parallel correlation function reads where G(η, τ ) is the scaling function at leading order in finite size corrections and C(η, τ ) is a finite-size effect due to the interface structure. The leading-order term for large R can be expressed as a single integral where C(η, τ ) = ∂ η G(η, τ ) is given by and χ = η/ √ 1 − τ 2 . Let us discuss now the perpendicular correlation function. For spin fields widely separated from each other (md ≫ 1) and from the wall (mx ≫ 1) the perpendicular correlation function is The clustering property is easily tested by performing the limit md → ∞, which yields Anticipating some results, we will show that the clustering to the one-point correlation function shown above is confirmed by the numerical simulations provided in the next section.

Magnetization profile
Before presenting the comparison between theory and numerical simulations, we summarize the details of the numerical implementations. We performed MC calculations by using a hybrid scheme which combines the standard Metropolis algorithm (see, e.g. [50]) and the Wolff cluster algorithm [51]. The minimum number of MC steps per site is 10 7 . Parallelization was obtained by independently and simultaneously simulating up to 128 Ising lattices on a parallel computer. An appropriately seeded family of dedicated, very large period, Mersenne Twister random number generators [52], in the MT2203 implementation of the Intel Math Kernel Library, was used in order to simultaneously generate independent sequences of random number to be used for the MC updates of the lattices. More in detail, each independent parallel simulation comprised: i) a warm up phase in which 10 9 numbers drawn from each random number stream in order to warm up the generators followed by a thermalization of the lattice of 3 × 10 6 × N Metropolis steps on random spins, with N the total number of spins in the simulation box. Then production phase steps started. Every production step a realization of the lattice was saved to disk for further processing. Every production step consisted on 10 5 hybrid updates consisting of a Wolff cluster update at a random point in the lattice and 0.5 × N Metropolis spin flips on random spins. The temperature is chosen such that the bulk correlation length ξ b is much smaller than the system size, i.e. ξ b ≪ R. The simulation box is a rectangle with (x, y) ∈ [0, l x ] × [−l y , l y ]. The vertical size is l y = H + R/2 with H/λ ≳ 4 while l x /λ ≳ 7. This design effectively reduces the influence of the boundaries and mimics the field-theoretical setups in the simulations. The first observable we examine is the magnetization profile. According to (2.2), the theoretical prediction for the rescaled profile along the x -axis is given by In figure 4, we plot the rescaled magnetization profile ⟨σ(x, 0)⟩ B +−+ /M as a function of the rescaled horizontal coordinate η = x/λ. For the bulk correlation length ξ b we used the exact expression provided below (2.1). MC data have been obtained for several combinations of the temperature T and system size R. We also emphasize that T and R and the only independent input data and that no adjustable or fit parameters are involved.
In order to compare simulations with field-theoretical results ξ b has to be much larger than the microscopic lattice spacing, this in order to ensure a continuum treatment in terms of fields. Nonetheless, it has been observed a posteriori that the good agreement persists even for correlation lengths of order of the lattice spacing; such a feature has emerged for the q-state Potts model [27] and for the three-dimensional Ising model [53].
We observe in figure 4 how the numerical results fall systematically away from the leading-order profile given by the scaling function Υ(η), a feature that has been reported in the literature [31]. The reason is due to the fact that the interface structure correction ∝ (ξ b /R) 1/2 is actually non-negligible, as anticipated. The quantitative agreement between theory and numerics improves drastically when the subleading correction at order R −1/2 is included. For the sake of completeness we mention that MC data in [31] have been obtained for T ≈ 1.9 and R ≈ 46 while in this paper we consider temperatures closer to the critical point (T = 2, T = 2.1, and T = 2.15) and larger system sizes (ranging from R = 151 to R = 451). The agreement is excellent for any x, provided x is not close to the boundary within a layer of thickness ≈ ξ b . This feature is actually compatible with the applicability domain of the field-theoretical result which, as mentioned, requires mx ≫ 1. Error bars are not indicated but typically their extent does not exceed that of the symbols in figure 4 and subsequent ones.
It is tempting to further test the validity of (3.1) by isolating the interface structure correction. The latter can be extracted from the numerical data by considering the quantity whose theoretical prediction is We recall that m = 1/(2ξ b ) and the exact expression for the bulk correlation length is provided below (2.1). Note that for R → ∞ the term B(x, 0) is actually the passage probability P (x, 0) up to an overall proportionality factor. The comparison between B(x, 0) extracted from the numerical data and the analytic prediction (3.3) is illustrated in figure 5. The overall agreement between theory and numerics is visible although the comparison is not accurate as for the magnetization profile shown in figure 4. The reason is due to the fact that numerical data are inevitably affected by further subleading effects at order R −1 which are amplified by the factor √ mR in (3.2). These corrections at order R −1/2 in (3.3) are responsible for the deviations visible in figure 5.
Another aspect that is often investigated in the context of interfacial phenomena is the liquid adsorption at a wall (see e.g. [7]). The adsorption is quantified in terms of the excess of density with respect to density attained deep in the bulk phase away from the wall. The quantity of interest is the coverage, which in our notations we define as . (3.5) In the following we will examine Γ 0 , which is the coverage computed from the density profile sliced along the x -axis, such a quantity can be directly measured in simulations. The theoretical prediction can be easily worked out from the analytical expression of the magnetization profile. A simple calculation entails we stress that Γ 0 depends on the temperature T via M and ξ b . We further observe that for large system size such that R ≫ ξ b the first term in (3.6) dominates and the

Spin-spin correlation function
Turning to correlation functions of the order parameter field, figure 7 provides the comparison of the theoretical and MC results for the parallel correlation function. The agreement is again noteworthy provided the correction term ∝ C(η, τ ) in (2.7) is taken into account. As expected, the agreement improves by taking system with larger size, as shown in figure 8.
It is instructive to further comment figure 7 in light of the analytic result obtained for separations ξ b ≪ y ≪ R [35], which is For fixed vertical separation y the x -dependence of (4.1) indicates an exponential envelope proportional to the passage probability density. This feature is signaling that interfacial correlations reach the maximal magnitude in correspondence of the most probable point where the interface passes. Secondly, upon varying y with fixed x the long-ranged character of correlations within the interfacial region is captured by the power-law growth proportional to y 1/2 .  An analogous comparison is carried in figure 9 for the perpendicular correlation function. In particular, we consider ⟨σ(x, 0)σ(x + d, 0)⟩ B +−+ as function of the rescaled coordinate η = x/λ for fixed separation d between spin fields. We observe an excellent agreement between theory and MC data for all the values of d we have examined. The clustering argument (2.11) is also successfully tested in figure 9. Equation (2.11), which corresponds to the limit d → +∞, is indicated with the red dot-dashed curve in figure 9. The above limiting curve is progressively approached for wide ranges of η by the MC data upon increasing the distance d, as shown by the curves with finite d ranging from  figure 9. Still focusing on the clustering for d → ∞, the leading order result Υ(η), which is obtained by ignoring the term proportional to P 1 in (2.11), and which is indicated with the dotted line, exhibits deviations from the MC data. Once more, we note that including the subleading correction at order R −1/2 is decisive in order to establish a quantitatively accurate agreement between theory and simulations.

Conclusions
In this paper, we have successfully illustrated the comparison between theoretical and MC results for order parameter profiles and correlations in the presence of a fluctuating droplet. The theory, which is exact in the near critical region, has been specified to the universality class of the two-dimensional Ising model, which we have simulated with high-precision MC simulations. All theoretical predictions are found in very good agreement with MC simulations in the absence of adjustable parameters, confirming thus the predictability power of the exact theory of phase separation. Although the comparison has been performed for the Ising model, there is no fundamental obstruction in testing other universality classes.
The salient feature that emerges from the analysis illustrated in this paper is the importance of interface structure effects. The system sizes and temperatures we considered clearly show how large the system has to be in order to achieve a quantitatively good comparison between theory and simulations. We show that when interface structure effects are added on top of the magnetization profile the old discrepancies between theory and numerics are resolved, as shown by the perfect agreement in figure 4. As discussed while commenting equation (4.1), order parameter correlations are found to be long-ranged along the interfacial direction but at the same time they are spatially confined within the interfacial region. This observation, which for us emerges directly from an exact field-theoretical results in real space, agrees with Wertheim's prediction based on the analysis of integral equation theory for inhomogeneous liquids [54] (see also [1,2,55] and [56] for experiments). Our analysis show that also interface structure yields long-ranged correlations in the presence of entropic repulsion. Although our investigation has been carried out in real space, it will be shown elsewhere how to pass in momentum space by extending the familiar notion of interface structure factor to strongly fluctuating droplets [57].