Four-point interfacial correlation functions in two dimensions. Exact results from field theory and numerical simulations

We derive exact analytic results for several four-point correlation functions for statistical models exhibiting phase separation in two-dimensions. Our theoretical results are then specialized to the Ising model on the two-dimensional strip and found to be in excellent agreement with high-precision Monte Carlo simulations.

Within the language of statistical physics, the characterization of thermal fluctuations exhibited by extended objects is inevitably formulated in terms of correlation functions of certain quantities. This is also the case for liquid-vapor interfaces. In a seminal paper by Wertheim [14] it was shown that density fluctuations in the presence of phase separation become long-ranged along the interface separating coexisting phases and that these correlations are confined within the interfacial region [2,15]. This situation happens to be in sharp contrast with the exponential decay of correlations which characterizes pure phases 1 [18,19]. The above features exhibited by correlations within interfacial region separating coexisting phases have been at the center of numerous studies based on coarse-grained descriptions in terms of effective interfacial models [20][21][22][23][24] followed from the so-called capillary wave model [25]. We refer to [26] for a historical account and to [22][23][24] (and references therein) for the state of the art on effective interfacial models.
Going beyond effective models, and focusing on exactly solvable models, for long time the study of interfacial behavior has been limited to the planar Ising model on the lattice [27]. The possibility of obtaining analytic results from the scaling limit of exact solutions on the lattice proved to be essential for the formulation of heuristic interpretations based on the analogy between fluctuating interfaces in two dimensions and random walks [28] as well as formulations based on Solid-On-Solid models [29,30] in which spin interfaces are identified with fluctuating Onsager-Temperley strings [31] and statistical sums are implemented by path integrals [32,33].
The two-dimensional case turns out to be very interesting because the exact analytic form of interfacial correlations can be obtained for several models of statistical mechanics [34]. In particular, it is possible to characterize in a mathematically precise form the anisotropic character of interfacial correlations, their long-range character, as well as their confinement, providing thus a quantitative analysis of the scenario outlined by Wertheim within the context of non uniform fluids [14,15].
In this paper, we show how to find closed-form expressions for certain four-point correlation functions of the order parameter field in the presence of phase separation. The systems considered in this paper are generic models of two-dimensional statistical mechanics at phase coexistence close to a second-order phase transition point. The analytic expressions derived in this paper follow from the exact theory of phase separation in two dimensions 2 [34,37,38] and wetting phenomena [39][40][41][42] and, more specifically, to recently obtained results for n-point correlation functions in the presence of an interface [43]. In particular, it has been shown in [43] how the first-principle field-theoretic calculation for n-point correlation functions admits an exact probabilistic interpretation in which the interface fluctuates as a Brownian bridge. The results obtained in this paper are thus self-contained and derived by following the probabilistic interpretation, which is exact at both leading and first subleading corrections in finite-size corrections.
This paper is structured as follows. In Sec. 2, we recall the probabilistic interpretation and we introduce the specific four-point correlation functions considered throughout this paper. We then derive the corresponding analytic expression for each order parameter correlator and discuss some analytic properties they satisfy. In Sec. 3, we present the comparison between the analytic expressions and the results obtained from Monte Carlo simulations for the Ising model. Conclusive remarks are summarized in Sec. 4. Appendix A collects some mathematical details involved in the calculations presented in Sec. 2.

Interfacial correlations
We consider a two-dimensional statistical system at phase coexistence close to a second order phase transition point. To be definite, the system we consider is a ferromagnetic spin model defined on the two-dimensional strip (x, y) ∈ R × (−R/2, R/2) with finite width R in the y direction. Boundary conditions with spins fixed in state a for x < 0 and b for x > 0 are used to enforce the emergence of an interface separating coexisting phases; see Fig. 1.
The quantity which we are going to examine is the four-point correlation function of the order parameter field σ(x, y) where · · · ab stands for statistical averages in the strip geometry with ab boundary conditions corresponding to Fig. 1. It has been shown in Ref. [43] how general results for n-point interfacial correlation functions with arbitrary n can be derived from the underlying field theory associated to the scaling model. In this paper, we specify the above mentioned results to n = 4 with the four points arranged in particularly symmetric configurations. This choice will allow us to write the exact results of Ref. [43] in an explicit form which is particularly suited for both analytical considerations as well as for comparison with numerical simulations. The configurations we will examine throughout this paper are those depicted in Fig. 1, and the corresponding four-point correlation functions are those defined in (2.2)-(2.4). The correlation function G (y 1 , y 2 , y 3 , y 4 ) = σ(0, y 1 )σ(0, y 2 )σ(0, y 3 )σ(0, y 4 ) ab (2.2) corresponds to order parameter fields aligned along the straight line which connects the pinning points, as illustrated in Fig. 1a. The correlation function corresponds to order parameter fields arranged in a rhomboidal pattern with diagonals of length 2x and 2y parallel to the coordinate axes, as indicated in Fig. 1b. Lastly, the correlation function corresponds to order parameter fields arranged in a rectangular pattern with edges of length 2x and 2y parallel to the coordinate axes, as depicted in Fig. 1c. Throughout this manuscript we will refer to G (y 1 , y 2 , y 3 , y 4 ), G (x, y) and G (x, y) as the parallel, rhomboidal and rectangular correlation functions, respectively. We are now in the position to introduce the theoretical framework employed to calculate the correlation functions we are interested in. To this end, we recall those results obtained in [43] which are essential for the discussion in the present paper. It has been shown in [43] that exact results for n-point correlation function -such as (2.1) when n = 4 -can be represented within a probabilistic interpretation in which the correlation function is reconstructed by averaging over interfacial configurations weighted with the probability density of a Brownian bridge. This probabilistic interpretation is valid for those universality classes and boundary conditions in which phase separation in two dimensions occurs through a single interface [37,38,44].
More specifically, the probabilistic interpretation allows us to express the correlation function (2.1) as follows (2.5) Some comments are in order. The subscript j in σ j denotes the j-th spin field entering in (2.5). In the more general setting, the spin field σ j can have more than one component and j can be used to label them. For instance, this situation happens in the quantum field theory associated to the scaling q-state Potts model [37]. Then, P 4 du 1 du 2 du 3 du 4 is the probability for the interface to cross all the intervals (x j , x j + dx j ) at y = y j with j = 1, 2, 3 and 4. As shown in [43], the joint passage probability P 4 is the one which characterizes a Brownian bridge in which the random walker starts its motion in x = 0 at time y = −R/2 and comes back to x = 0 at time y = R/2 (the specific form of P n is provided in (2.12)). The profile gives the magnetization at point (x, y) when the interface is conditioned to cross the point (u, y). In (2.6), θ(x) denotes Heaviside unit step function, i.e., θ(x) = 1 if x > 0 and θ(x) = 0 if x < 0, and δ stands for Dirac delta function. The first two terms in the right hand side of (2.6) correspond to a step profile in which the two coexisting phases are sharply separated by the interface. According to this picture, the interface sharply separates two regions in which the order parameter exhibits the vacuum expectation values σ i a and σ i b . The subsequent term proportional to the Dirac delta takes into account the first correction beyond the picture in which the interface is regarded as a structureless entity. The interface structure term originates a subdominant correction which contributes at order R −1/2 and whose overall factor A (σ i ) ab|i is known for integrable field theories [37]; we refer to [45] for recent numerical simulations. We also mention that subsequent corrections can be systematized in an exact expansion in powers of the small parameter ξ/R and that corrections at order R − /2 with = 2, 3, . . . can be obtained within the field-theoretical method outlined in [43]. Since in this paper we are interested in the leading-order form of interfacial correlations, from now on we retain the first two addends in (2.6).
Let us discuss the general structure of the leading-order result for the correlation function. It is straightforward to realize how the Heaviside theta functions in (2.5) yield a representation in which the correlation function is expressed through cumulative distributions of the probability density P 4 . It turns out to be convenient to write the first two terms in the right hand side of (2.6) in the form with the jump of expectation values in pure phases given by ∆ σ i = σ i a − σ i b . The joint passage probability P 4 satisfies the following properties: P 3 (u 1 , y 1 ; u 2 , y 2 ; u 3 , y 3 ) =ˆR du 4 P 4 (u 1 , y 1 ; u 2 , y 2 ; u 3 , y 3 ; u 4 , y 4 ) (2.8) P 2 (u 1 , y 1 ; u 2 , y 2 ) =ˆR 2 du 3 du 4 P 4 (u 1 , y 1 ; u 2 , y 2 ; u 3 , y 3 ; u 4 , y 4 ) (2.9) P 1 (u 1 , y 1 ) =ˆR 3 du 2 du 3 du 4 P 4 (u 1 , y 1 ; u 2 , y 2 ; u 3 , y 3 ; u 4 , y 4 ) (2.10) 1 =ˆR 4 du 1 du 2 du 3 du 4 P 4 (u 1 , y 1 ; u 2 , y 2 ; u 3 , y 3 ; u 4 , y 4 ) ; (2.11) the first three equations give the marginal distributions and the last one is the normalization condition. Analogously to P 4 , the joint probabilities P n express the probability density for the crossing of n intervals. It has been shown in [43] that that P n can be expressed in terms of the n-variate normal distribution Π n [46,47] through with the following notations the parameter m, which is the kink mass in field theory, and the surface tension in the language of liquid state theories [37], is related to the subcritical bulk correlation length ξ via ξ = 1/(2m). Then, Π n (u 1 , . . . , u n |R 1...n ) is the n-variate normal distribution with (symmetric) correlation matrix R 1...n . The ij matrix element of the correlation matrix ρ ij = (R 1...n ) ij is given by 14) The elements in the lower triangular part follow by recalling that R 1...n is a symmetric matrix. Note that we are taking y i > y j for i < j, and so also τ i > τ j for i < j. Expectation values with measure Π n are defined by The probability density is standardized, meaning that E[u 2 i ] = 1 and The normalization condition then ensures that E[1] = 1. The last quantity we need to introduce is the cumulative function of the n-variate normal distribution The above allows us to cast integrals arising from the probabilistic interpretation in terms of cumulative functions. For instancê By proceeding along the above lines, the four-point correlation function (2.1) reads where all spin fields carry the same index i. The most general case in which different components of the order parameter occur in (2.19) can be treated along the same lines. It is worth pointing out how clustering properties of correlation functions can be derived from (2.19). For instance, by taking x 4 → ∓∞, we have the following limiting behavior lim (2.20) where the three-point correlation function agrees with field-theoretic calculations [43,48]. The remaining part of this section is devoted to a detailed case-by-case analysis of (2.19) with spin fields arranged as shown in Fig. 1. Such an analysis will allow us to find explicit expressions for (2.19) involving either single integrals or certain special functions related to the Gaussian distribution [49,50].

Parallel correlation function
For spin fields arranged in the configuration of Fig. 1a, we have x j = 0 with j = 1, . . . , 4 and −R/2 < y 4 < y 3 < y 2 < y 1 < R/2. Since χ j = 0 in (2.19), the calculation of the parallel correlation function involves orthant probabilities only. We introduce the following shorthand notation for orthant probabilities: Clearly, F 1 (0) = 1/2. We recall that R 1...n is a n × n symmetric matrix with 1 along the main diagonal, therefore it comprises n(n − 1)/2 nontrivial entries in the upper triangle. Thanks to the Markov property 3 , ρ ij ρ jk = ρ ik for i < j < k. It thus follows that the number of independent entries is actually lowered to n − 1. For this reason the orthant probability for the trivariate normal distribution depends on two independent variables and, analogously, the orthant probability for the quadrivariate normal distribution depends on three independent variables. Without loss of generality, we can take for the independent correlation coefficients those which appear in the superdiagonal entries, namely, ρ i,i+1 for i = 1, . . . , n − 1. Thanks to (2.19) and (2.22) the parallel correlation function reads The orthant probabilities Φ j for the bivariate (j = 2) and trivariate (j = 3) normal distributions are given by [46,47] where ρ 13 = ρ 12 ρ 23 because of the Markov property which characterizes the Brownian bridge. The expression of the orthant probability for the quadrivariate normal distribution is considerably more involved and reads [53] Φ 4 (ρ 12 , ρ 23 , ρ 34 ) = 1 16 (2.26) By plugging (2.24)-(2.26) into (2.23), and introducing 28) where C and G are the scaling functions and G (y 1 , y 2 , y 3 , y 4 ) = 4 π 2 sin −1 (ρ 12 ) sin −1 (ρ 34 ) (2.30) The above provides an exact result for the parallel correlation function when order parameter fields are placed along the interface, as shown in Fig. 1a.
In the following, we further specialize the parallel correlation function to spin fields symmetrically arranged with respect to the horizontal axis. We set y 1 = −y 4 , y 2 = −y 3 and observe that this choice implies ρ 12 = ρ 14 in (2.29) and (2.30). Thus, we introduce the above reads (2.32) with the scaling functions and (2.34) We defer the examination of some limiting cases in which the parallel correlation function reduces to the rhomboidal or rectangular one in Sec. 2.2 and Sec. 2.3, respectively.

Rhomboidal correlation function
The analysis for the correlation function with spin fields in rhomboidal arrangement does not proceed along the lines illustrated in the previous section. Firstly, because we don't provide a general expression for F 4 (X 1 , X 2 , X 3 , X 4 |R 1234 ) which we could use in the limit of interest, secondly because the correlation coefficient degenerates to unity, i.e., ρ 23 = 1; such a limit has to be handled carefully. The perfect correlation between the random variables u 2 and u 3 -as signaled by ρ 23 = 1 -is best treated by taking the limit ρ 23 → 1 in the passage probability density (2.12). The correlation function then follows by plugging the passage probability for the rhomboidal arrangement, hereinafter denoted P ( ) 4 , into (2.19). This is the calculation strategy which we are going to detail.
The passage probability for the rhomboidal arrangement is obtained from the following limiting procedure P ( ) As a result, the passage probability P ( ) 4 is characterized by one correlation coefficient, say α; note that ρ 14 = ρ 12 ρ 23 ρ 34 = α 2 by virtue of ρ 23 = 1. Recalling the relationship between P 4 and Π 4 given by (2.12), a simple calculation yields for the rescaled passage probability the following expression is correctly normalized. This fact can be confirmed by a direct calculation: as we anticipated and as consistency requires. The four-point correlation function is computed by plugging the passage probability (2.37) into the probabilistic representation (2.19). Leaving in Appendix A the lengthy mathematical manipulations, the result for the rhomboidal correlation function reads where C (η, τ ) and G (η, τ ) are the scaling functions where T is Owen's function (see (A.5) for its definition) and Y(η, r) is defined by where erf(z) = (2/ √ π)´z 0 dt e −t 2 is the error function [54], η = x/λ is the rescaled horizontal coordinate, and Let us discuss some general properties of the rhomboidal correlation function. By definition, G (x, y) is an even function under the exchange of x with −x. The parity G (x, y) = G (−x, y) is thus satisfied by both the scaling functions (2.43); recall that T ( √ 2η, r) = T (− √ 2η, r) (see (A.5)). The limit x → +∞ for the rhomboidal correlation function projects the order parameter field σ i (x, 0) deep into the b phase and, analogously, σ i (−x, 0) deep into the a phase. As a result, we obtain the following clustering relation The quantity enclosed in square brackets is the two-point correlation function with order parameter fields placed along the interface, i.e., σ i (0, y)σ i (0, −y) ab [34]. As a cautionary check, it is useful to compare the limit x → 0 of the rhomboidal correlation function with the limit y 2 → 0 of the parallel correlation function. Consistency requires that the above limits have to coincide, i.e. lim x→0 G (x, y) = lim In order to check that (2.47) is satisfied, we firstly examine the limit x → 0 of the scaling functions which appear in the correlation function G (x, y). For the scaling function C (η, τ ), we have Let us consider now the right hand side of (2.47). The limit y 2 → 0 of the scaling function C (sym) (τ, τ 2 ) gives (see (2.33)) lim which coincides (2.48), as it should. Finally, we consider the same limit discussed above but now for scaling function G (sym) (τ, τ 2 ). From (2.34), we have lim the last equality follows from the identity (A.14), whose proof is supplied in Appendix A. Summarizing, thanks to the results (2.48)-(2.51), we have proved the connection between the parallel and rhomboidal correlation functions given by (2.47). The scaling function G (η, τ ) is plotted in Fig. 2 as function of η and τ .

Rectangular correlation function
The correlation function with spin fields arranged as shown in Fig. 1c can be calculated by following the procedure outlined in Sec. 2.2. Leaving the mathematical details in Appendix A, the rectangular correlation function admits the following representation with the scaling functions and the rescaled coordinates η = x/λ, χ = η/κ, κ = √ 1 − τ 2 , τ = 2y/R. Analogously to the rhomboidal case, reflection symmetry is expected, namely: G (x, y) = G (−x, y). The limit x → +∞ in (2.52) entails a separation of the four spin fields into two infinitely separated clusters. By using the properties of Owen's T function listed below (A.5) and the limits which is the expected clustering behavior. Then, the limit x → 0 of the rectangular correlation function must reduce to the limit y 2 → y 1 ≡ y in the parallel correlation function, i.e., lim x→0 G (x, y) = lim which is indeed the case. The identity (2.57) can be easily verified by observing that the scaling functions which appear in the left hand side of (2.57) satisfy the limits while for the scaling functions in the right hand side, we have lim By virtue of the identity the right hand side of the first equations in (2.58) and (2.60) coincide, and the property (2.57) follows. The analytical examination of interfacial correlations can be carried out by taking the limit of small τ in (2.53). Among the two scaling functions in (2.53), we focus on G (η, τ ) because it is the one which characterizes the Ising model. The analysis we are going to provide for G (η, τ ) can be straightforwardly extended to the scaling function C (η, τ ). The analytical study is better achieved by writing the scaling function G (η, τ ) in the equivalent form which follows from the functional identity (A.8) satisfied by the T -function; in the above, erfc(z) = 1 − erf(z) is the complementary error function [54]. The scaling function G (η, τ ) is plotted in Fig. 3. As already anticipated by the study of clustering properties, G (η, τ ) exhibits a nontrivial dependence through the coordinates η and τ only within the interfacial region, i.e., when the rescaled horizontal coordinate η is not large and the vertical one does not vanish (τ > 0). For small vertical separations τ the scaling function (2.63) can be approximated with the expansion  Figure 3: The scaling function G (η, τ ). Dashed curves indicate contour lines corresponding to the values ranging from 0.9 to 0.1 with spacing 0.1 from outer to next to inner curve, and 0.01 for the inner curve.
up to higher-order corrections of order O(τ 3/2 ). The small-y expansion provided by (2.64) accurately encompasses any value of η both within and without the interfacial region provided ξ y R. The term √ τ ∝ y/R is the signature of long range correlations in the direction parallel to the interface. The term enclosed in square brackets provides a non-trivial spatial dependence which is specific for the rectangular arrangement of spin fields.
The T -functions in (2.63) may not allow for a direct visualization of the dependence through the horizontal coordinate η, the coordinate transverse to the interface. Alternatively, we can examine the slope in the horizontal direction, the latter is given by the first derivative of G (η, τ ) with respect to η The overall exponential envelope e −χ 2 suppresses interfacial correlations far away from the interfacial region (|χ| 1) and determines the confinement of the long-ranged character of density fluctuations in the region where |χ| 1, where it is more probable to find the interface. We conclude this section with a discussion on a general property satisfied by correlation functions. It is known from rigorous studies that fluctuations of the interface midpoint grow as √ R [55,56] and that interface fluctuations in Ising [57] and q-state Potts models [58] are characterized by the Brownian bridge property. This feature is also well visualized by the explicit form of the passage probability P 1 (x, 0) = π −1/2 λ −1 e −(x/λ) 2 , which gives for the mean square position of the interface midpoint In the limit R → ∞ with fixed x and y both the rescaled variables variables η = x/ √ Rξ and τ = 2y/R tend to zero. As a result, the scaling functions encountered in this paper satisfy the following relations which can be obtained by evaluating the scaling functions for η = τ = 0. The above limits imply (2.68) meaning that the limit R → ∞ of the four point correlation function yields an averaging over coexisting phases. The above limit is actually a particular case of a more general result which holds for any n and generic universality classes [43], and includes an analogous result for the Ising model [27] as a particular case.

Numerical results
The analytical predictions obtained in the previous sections are specialized to the Ising model and compared against Monte Carlo simulations. We set the notation by recalling the lattice Hamiltonian where s i = ±1 and the sum is restricted over nearest neighboring sites of a two-dimensional square lattice. The critical temperature is given by T c /J = 2/ log(1 + √ 2) = 2.269 185 . . . [59]. The correlation length 4 in the low-temperature phase reads with the dual coupling K defined by means of exp(−2K ) = tanh K with K = J/T [27]. Without loss of generality, we set the ferromagnetic coupling J = 1 in our simulations. We consider boundary conditions with a = −1, b = +1 and recall that Ising symmetry implies σ + = − σ − . The spontaneous magnetization M = σ + > 0 is given by [60,61] From the above boundary conditions it follows that σ = 0 and ∆ σ = −2M . Focusing on the leading-order result, the parallel, rhomboidal and rectangular correlation functions reduce to: 4) with G , G , and G the scaling functions (2.34), (2.43), and (2.53), respectively. We recall that ρ 12 and ρ 23 in (2.34) are given by (2.14) with τ 2 = −τ 3 . Then, η and τ in (3.4) are defined below (2.53). We emphasize that (3.4) are exact results at leading order in finite-size corrections. Interface structure corrections proportional to R −1/2 occur in general but actually vanish for the Ising model [43]. As a result, (3.4) are valid up to subleading corrections at order O(R −1 ).
The numerical simulations have been carried out on a rectangular lattice of horizontal length L and temperature T such that ξ R. The theory is defined for L → ∞ but in practice it is enough to take L/λ 7 as a criterion for the sizing of the simulation box. We summarize the details of the simulation scheme which we have already employed in a companion paper [48]. The simulations are performed by means of a hybrid Monte Carlo scheme (see, e.g. [62]) which combines the standard Metropolis algorithm and the Wolff cluster algorithm [63]. 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 [64], 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.
The first observable we consider is the parallel correlation function. The analytic result (2.34) and the numerical data are compared in Fig. 4 as function of τ 2 = 2y 2 /R for fixed τ 1 = 2y 1 /R. For small values of y 2 the correlation function approaches the asymptotic value given by lim The above result is actually confirmed by the numerical simulations. For each value of y 1 sampled in Fig. 4, τ 1 provides the maximum value allowed for τ 2 . This maximum value is reached when the order parameter field σ(0, y 2 ) approaches σ(0, y 1 ) from below; see Fig. 1a. In such a regime the correlation function approaches the limiting value lim a feature which is well visible in the plot of Fig. 4. It has to be emphasized that such a limiting result is obtained in field theory from the single-particle contribution which dominates the asymptotic of the large-separation behavior of correlation functions and that further corrections arising from interface structure effects are expected at order O(R −1 ) for the Ising model [43].
Therefore the limit y 2 → y 1 has to be understood in the sense that the separation y 1 − y 2 is small compared to R but it has to be large compared to the bulk correlation length, i.e., ξ y 1 − y 2 R. With this in mind, we actually observe an excellent agreement between theory and numerics. The next quantity we consider is the correlation function for the rhomboidal pattern. The comparison between theory and numerics is presented in Fig. 5. In Fig. 5a, we plot the correlation function G (η, τ ) as function of the rescaled coordinate η = x/λ for several values of τ = 2y/R; the vice versa is done in Fig. 5b. The analytic result for G (η, τ ) is provided by (2.43) with α = (1 − τ )/(1 + τ ). Focusing firstly on Fig. 5a, we observe that for small η the correlation function approaches the limiting value while for η → +∞ the asymptotic value is the opposite of the one attained for η → 0. The origin of this symmetry can be easily understood by noticing that in both cases the spin fields σ(−x, 0) and σ(x, 0) probe either the same phase, or oppositely magnetized phases. Another interesting feature is the rapidity upon which the above limiting values are attained. By taking a first derivative of G (η, τ ) with respect to η, we find meaning that the decay of correlations is exponentially fast along the x axis, in analogy with the rectangular correlation function; see (2.65). Such a behavior happens to be in sharp contrast with the algebraic form which characterizes the correlation along the direction parallel to the interface (y-direction); see the power law ∝ y/R in (2.64). The long-range form of interfacial correlations can be investigated analytically by performing an asymptotic analysis of the function Y(|η|, r) for small y/R. This task has been recently carried out in the examination of the threepoint correlation function σ(0, y)σ(x, 0)σ(0, −y) −+ [48]. In fact, it has been shown in [43] that the above mentioned three-point correlation function for the Ising model is proportional to the scaling function Y(|η|, r) which occurs in the rhomboidal correlation function. We thus refer the interested reader to Ref. [48] for a detailed mathematical account on this matter. The remarkable agreement between theory and numerics is confirmed by Fig. 5b where the correlation function is plotted as function of τ for fixed η. The deviations for small τ which occur in the plot of Fig. 5b are actually expected since that regime is at the boundary of the applicability domain of the analytic result at leading order in powers of (ξ/R) 1/2 . We refer to [34,43] for a detailed discussion about the theoretical treatment of interface structure corrections and to [45,48] for numerical simulations. Lastly, we discuss the comparison between theory and numerics for spin fields in the rectangular arrangement. The qualitative features discussed for the rhomboidal pattern occur also for the rectangular one. The analytic expression for G (η, τ ) given in (2.53) turns out to be in good agreement with the numerical results shown in Fig. 6. The small deviations which are visible in the plot of Fig. 6b have the same nature of those arising in Fig. 5b.
Although in the field-theoretical approach the bulk correlation length is supposed to be large compared to the lattice spacing a, and small compared to the system size, i.e. a ξ R, numerical data are found to be in quantitative agreement even when the constraint ξ a is not satisfied. For the temperature T = 2 considered in the simulations shown in this paper, the correlation length slightly exceeds two lattice spacings. An analogous agreement between theory and simulations for correlation functions within the interfacial region of an Ising interface has been observed for T = 2.15, corresponding to ξ ≈ 5 lattice spacings. Furthermore, the striking agreement between theory and numerics even when ξ does not exceeds few lattice spacings has been already observed in recent simulations for the q-state Potts model [45].

Conclusions
In this paper we have illustrated how to find exact closed-form expressions for certain four-point correlation functions for statistical systems exhibiting phase separation on a strip of finite width R ξ, with ξ the bulk correlation length. The analytic results derived in this paper are selfconsistent with the sole exception of equations (2.5), (2.6) and (2.12), corresponding to the exact probabilistic representation of n-point correlation functions which has been recently established from a field-theoretic calculation in Ref. [43]. We have examined correlation functions with order parameter fields arranged in three configurations: along the interface, in a rhomboidal pattern, and in a rectangular pattern; see Fig. 1. For each pattern the corresponding correlation function is characterized by two scaling functions whose analytic expression is provided in closed form. These scaling functions are universal in the sense that they are shared by models which exhibit phase separation through a single interface. The exact analytic expressions for the correlation functions provided in this paper allow for a direct examination of general properties of interfacial correlations. In particular, the long-range character of correlations in the direction parallel to the interface and their confinement within the interfacial region follow from a direct inspection of the analytic results we provided. In the second part of the paper we have specialized the analytic results to the Ising model. We have shown that Ising symmetry selects only one of the above mentioned scaling functions (see (3.4)) and that analytic results are confirmed by high-precision Monte Carlo simulations we performed.

A Mathematical details A.1 The correlation function G
The calculation of G involves certain integrals arising from the cumulative distribution for the passage probability Π . In this appendix, we adopt some tricks which enable us to shorten the calculation considerably. Firstly, we write the sharp magnetization profile as follows Then, we introduce the following shorthand notation for integrals with measure Π : where f is a function of u 1 , . . . , u 4 and f stands for its the expectation value. The result of (A.2) is a function of x and y but, as we shall see in a while, the dependence through x occurs by means of the rescaled coordinate η = x/λ while the dependence through y is encoded in α = (1 − τ )/(1 + τ ), with τ = 2y/R. Thanks to these notations, the correlation function we are interested in can be written as follows and the functional identity [49] T (h, a) + T (ah, 1/a) = 1 2 g(h) + 1 2 g(ah) − g(h)g(ah) , a 0 g(x) = 1 + erf(x/ √ 2) 2 .

A.2 The correlation function G
Let us consider now the rectangular correlation function. Since for this configuration y 1 = y 2 and y 3 = y 4 the corresponding correlation coefficients degenerate to unity, i.e., ρ 12 = 1 and ρ 34 = 1, respectively. The corresponding passage probability is found by taking the double limit ρ 12 → 1 and ρ 34 → 1, the latter yields where Π 2 (u 1 , u 4 |ρ 23 ) is a bivariate normal distribution with correlation coefficient ρ 23 . By inserting (A.9) into the probabilistic representation for the correlation function, we find has been adopted. By using the integral representation (A.7) of Owen's T function, a simple (but rather long) calculation entails (A.13) The results (2.52) and (2.53) given in the main body of the paper follow straightforwardly.

A.3 A useful identity
Here, we prove the following identity (A.14) with 0 x 1. By taking the first derivative with respect to x, we obtain the integral in the third term readŝ 16) thus, it follows that The identity (A.14) is established by integrating back with respect to x and fixing the integration constant by using the boundary condition, f (0) = 0, or f (1) = 1.