Multipoint correlation functions at phase separation. Exact results from field theory

We consider near-critical two-dimensional statistical systems with boundary conditions inducing phase separation on the strip. By exploiting low-energy properties of two-dimensional field theories, we compute arbitrary $n$-point correlation of the order parameter field. Finite-size corrections and mixed correlations involving the stress tensor trace are also discussed. As an explicit illustration of the technique, we provide a closed-form expression for a three-point correlation function and illustrate the explicit form of the long-ranged interfacial fluctuations as well as their confinement within the interfacial region.


Introduction
The characterization of the fluctuating interface separating coexisting phases is a longstanding problem in classical statistical mechanics [1]. A conspicuous amount of investigations in the field of interfacial phenomena has been stimulated by the need of a theoretical understanding and also because of numerous technological applications triggered by capillary forces and wetting effects at the nanoscale [2]. The current understanding of interfacial behavior benefitted from several theoretical approaches based on microscopic descriptions formulated within lattice models, effective models, renormalization group, and numerical simulations; we refer to [2][3][4][5] for general reviews on the subject.
The two-dimensional case holds a central role because of the availability of non-perturbative techniques which lead to exact solvability of the strongly fluctuating regime. Among exactly solvable planar systems, the Ising model occupies a predominant position in the above framework due the existence of exact solutions based on the diagonalization of the transfer matrix for a wide class of boundary conditions, which, in turn, lead to the formation of interfaces in certain planar geometries [6][7][8][9]. Results based on the scaling limit of exact solutions for the planar Ising model [9] motivated the introduction of phenomenological approaches relying on the analogy between interfaces and random walks [10]. Such an analogy has been exploited in order to construct effective coarse-grained descriptions based on the so-called Solid-On-Solid (SOS) models [11,12]. In these models, spin interfaces are identified with fluctuating Onsager-Temperley strings [13]. The equilibrium statistical mechanical problem in two dimensions is thus equivalent to a quantum mechanical problem in one space dimension. Within such an analogy, the construction of the partition function by summing over interfacial configurations is thus mapped onto the evaluation of path summations in the quantum-mechanical picture [11,14,15].
In recent years, an additional analytical framework allowed for the full exploration of interfacial phenomena in near-critical systems belonging to a wider range of universality classes in two dimensions [16], ranging from the Ising model, the q-state Potts model, as well as the Ashkin-Teller model, and other models which exhibit interfacial wetting and phase separation through intermediate phases [17]. The versatility of the field-theoretical formalism then allowed also for the investigation of the interplay between interfacial fluctuations and entropic repulsion due to a flat wall [18], a defect line [19], a wedge-shaped boundary [20,21], and its corresponding wetting/filling transition. Some of the above mentioned exact findings have been successfully tested by means of high-precision Monte Carlo simulations [22].
In most of the cases, as the ones mentioned above, the knowledge of one-point correlation functions in certain geometries is informative enough and therefore it suffices for an adequate description. A more refined characterization of the interfacial behavior, however, demands the knowledge of higher-order correlation functions, with the pair correlation as the simplest nontrivial representative. The occurrence of long-range correlations within the interfacial region has been established within the context of theories of inhomogeneous fluids [23][24][25][26], effective models such as the so-called capillary wave theory [27], and full scale numerical simulations [28,29]. The fact that the above studies refer either to space dimension d 3, or to effective interface models arising from specific assumptions, prevents the straight application of them to the strongly fluctuating regime which characterizes the scenario in d = 2.
In order to avoid the introduction of ad-hoc assumptions -which intrinsically characterize any effective model -the exact investigation of the two-dimensional case has to be inevitably built on a formalism which is based on the truly fundamental degrees of freedom of the system. This first-principles-based viewpoint is at the basis of the results of [30]; there, it has been showed how field theory can be used in order to extract the exact form of interfacial correlations in real space. As a result, by examining the structure of the pair correlation function of the order parameter, it has been possible to exhibit the specific form of long-range correlations generated by phase separation.
In this paper, we show how the field-theoretic formalism [30] can be extended to the calculation of n-point correlation functions for arbitrary n. We will illustrate firstly the case of correlations of the order parameter field, for which we will single out those contributions which are originated by interfacial fluctuations from those which are genuinely due to bulk fluctuations. Then, we will identify the exact analytic form of order parameter correlations and interpret them according to a probabilistic picture. Closed form expressions for certain three-point correlation functions are written as an explicit illustration of the technique, which then is further specialized to the cases of the Ising and q-state Potts models. The theoretical analysis of correlation functions is then pushed in order to capture subleading finite-size corrections. It is then showed how the treatment of such corrections can be systematized in a power series of the small parameter (ξ b /R) 1/2 , with ξ b the bulk correlation length and R the separation between the interface endpoints. Specific results about corrections at order R −1/2 and R −1 are discussed in great detail and are related to the probabilistic picture, the latter amounts to interpret those terms as effects due to interface structure. The stress tensor trace Θ and its n-point correlations are also considered, as well as mixed correlators involving both Θ and spin fields.
This paper is organized as follows. In Sec. 2, we laid the basis for the calculation of npoint correlation functions of the order parameter profile. The calculation is thus broken into successive steps which are structured into subsections. The connected part of the correlation function, which is examined in Sec. 2.1, is expressed through n-body cluster functions, the latter are constructed explicitly in Sec. 2.2. The disconnected parts of the correlation function are computed in Sec. 2.3. The full result for the correlation function is thus supplied in Sec. 2.4 in general terms together with explicit applications to Ising and q-state Potts models. Secs. 3 and 4 deal with subleading corrections at orders R −1/2 and R −1 , respectively, and their emergence within the probabilistic picture. Conclusive remarks and a summary with a description of future perspectives is outlined in Sec. 5. Two appendices contain additional mathematical details related to the buildout of the material covered in Sec. 2.

Spin correlations: field theoretical derivation
We illustrate the calculation of the n-point correlation function of the spin field on the finite strip showed in Fig. 1. The quantities of interest are the correlation functions G n (x 1 , . . . , x n ) = σ 1 (x 1 ) · · · σ n (x n ) ab , (2.1) in which the notation · · · ab stands for the statistical average on the strip with ab boundary conditions (see Fig. 1) and σ j (x j , y j ) ≡ σ j is the spin field in the point x j = (x j , y j ) ∈ R × (−R/2, R/2) on the Euclidean plane. The following ordering is considered in (2.1): y 1 > y 2 > · · · > y n . The subscript j in σ j labels the j-th spin field; in the most generic case the spin field carries a color index such as in the field theory associated to the scaling q-state Potts model (see [31] and Sec. 3.4). The switching of boundary condition from a to b at x = 0 along the edges y = ±R/2 is . . . Figure 1: The strip geometry with ab boundary conditions. The spin fields which define the correlation function G n are illustrated with green circles.
implemented within the field-theoretical language through the boundary state formalism. The correlation function (2.1) is thus written as follows is the partition function and |B ab is the boundary state with the inhomogeneous boundary condition shown in Fig. 1. We refer to [32] for translationally invariant boundaries in the framework of massive integrable field theories. The boundary state is expanded in the basis of bulk excitations compatible with the topological charge of the boundary, namely where |K ab (θ) is the kink state corresponding to a topological particle with mass m which interpolates between vacua |Ω a and |Ω b , and θ is the rapidity variable. Within the dictionary of phase separation, the vacuum |Ω a is identified with the homogeneous system filled by phase a (pure phase). The second and subsequent terms on the right hand side of (2.4) stem from the propagation of multi-kink states. The simplest of them -the double-kink state |K ac K cbplays the dominant role only when the single-kink state |K ab is absent in the expansion of the boundary state. For those models and boundary conditions in which such an instance happens, the double-kink state is responsible for the formation of a double interface with an intermediate layer of phase c adsorbed on the ab interface [17]. We refer to [19] for a further characterization of models which exhibits such a phenomenology.
The amplitudes f ab (θ), f abc (θ, θ ′ ) depend on both the bulk and boundary universality classes and are known for certain integrable field theories [33][34][35]. For the purposes of this paper, however, it is sufficient to know only the infrared properties of the amplitude f ab (θ), which is responsible for the emission of a single-kink state from the boundary. To this end, we only need the Taylor expansion around θ = 0, which reads f ab (θ) = f ab (0) + O(θ 2 ) for small rapidities. The absence of the linear term in θ follows by reflection symmetry around the vertical axis, f ab (θ) = f ba (−θ), in conjunction with the symmetric role played by a and b, i.e., f ab (θ) = f ba (θ).
In the limit mR ≫ 1, which we consider from now on, the partition function (2.3) can be computed straightforwardly thanks to a saddle-point calculation. The result is up to higher order terms due to multi-kink states. From the boundary state formalism it is possible to identify the surface tension Σ ab associated to the creation of an interface separating the coexisting phases a and b. The surface tension is computed as the excess free energy 1 per unit length and is defined through the limit where Z a (R) stands for the partition function of the strip with uniform boundary condition a. The latter can be determined from the boundary state |B a for a uniform boundary. These states are known exactly for several models including, among the most relevant examples for the purposes of this paper, the Ising model with a surface field [32] and the q-state Potts model [36]. The state with lowest mass entering the expansion of the uniform boundary is actually the vacuum |Ω a ; hence, 2 Z a (R) = Ω a |Ω a ≃ 1 and the interfacial tension equals the kink mass, i.e., Σ ab = m . (2.7) In the symmetry broken phase the bulk correlation length ξ b is related to the kink mass m via ξ b = 1/(2m). Therefore, (2.7) implies Σ ab ξ b = 1/2. It is worth emphasizing that such a relationship is compatible with Widom scaling [37] in two dimensions. More interestingly, thanks to calculations based on the exact solution of the planar Ising model, the relation Σ ab ξ b = 1/2 has been proved to be valid for all subcritical temperatures. This result follows as an application of duality [9,38,39]; see also [40]. Let us move on the correlation function (2.2). The time ordering of spin fields we discussed above is actually implemented in a more strict sense since consecutive spin fields in (2.1) are separated by a distance much larger than the bulk correlation length ξ b and, in a similar fashion, spin fields are also taken to be far from the boundaries, hence, (x i − x j ) 2 + (y i − y j ) 2 ≫ 1/m, R/2 − y 1 ≫ 1/m and y n + R/2 ≫ 1/m. 1 In units of kBT . 2 The symbol ≃ stands for the omission of exponentially suppressed terms as in (2.5).
We apply a spectral decomposition which amounts to insert a resolution of the identity between each pair of spin operators. The identity operator is itself expanded in terms of kink states, something that we write as follows with kink states normalized according to By inserting the resolution of the identity between each spin field appearing in (2.2), and recalling translational invariance for bulk fields where H and P are the Hamiltonian and momentum operators in relativistic quantum field theory, i.e. H|K ab (θ) = m cosh θ|K ab (θ) , we can extract the space-time dependence of matrix elements and write (2.13) The dependence through the coordinates is encoded in the function (2.14) Matrix elements of the spin field M σ j ab are decomposed into a connected part and a disconnected one. The connected part is expressed in terms of the spin field two-particle form factor F σ j (2.15) By adopting a pictorial representation for matrix elements, (2.15) reads The disconnected part originates a Dirac delta corresponding to particle annihilation. The vacuum expectation value which multiplies the Dirac delta is σ j a if the two kinks are annihilated by passing right aside the spin field, as depicted in (2.16). Conversely, the overall vacuum expectation value is σ j b for the passage aside left. The right-left alternative is ultimately responsible for the presence of the annihilation pole 3 [31], whose behavior at small rapidity differences reads is the jump of vacuum expectation values across the ab interface. The coefficient c ab and the subsequent ones appearing in (2.17) are known for integrable field theories [31,42]. By plugging the expansion (2.16) into the product n j=1 M σ j ab (θ j |θ j+1 ), we obtain a decomposition of the npoint correlation function G n which comprises a connected part and a sequence of disconnected parts.
The rest of this section is organized as follows: in Sec. 2.1, we address the calculation of the connected part. Such a task will require the introduction of a certain class of special functions -which will be termed n-body cluster functions -whose definition and main properties will be provided in Sec. 2.2; there, we will also provide the result for the connected part of G n . The disconnected parts of the correlation function will be investigated in Sec. 2.3 and the full result for G n will be supplied Sec. 2.4.

Connected part and cluster functions
The connected part of the n-point correlation function G n (x 1 , . . . , x n ) is obtained from the product of the two-particle form factors F σ j aba (θ j − θ j+1 + iπ) in (2.12). Diagrammatically, this operation corresponds to stack the decomposition (2.16) for j = 1, . . . , n and retain the diagram in which all spin fields are connected, as illustrated in (2.19). The open necklace shown in the right hand side of (2.19) is the diagram which we will examine in this section.
According to (2.17), the leading low-energy behavior of the connected diagram on the right hand side of (2.19) is captured by the product of kinematical poles, therefore where O({θ} −n+1 ) stands for terms which are homogeneous functions of order −n + 1 in the rapidity variables. Subsequent terms in the above expansion lead to subleading finite-size corrections of the correlation function whose systematic analysis will be carried out in Sec. 3. The computation of the integrals in (2.12) proceeds by expanding the function U n at small rapidities. To this end it is convenient to rescale rapidities through the change of variables θ j → 2/(mR)θ j . The function U n becomes U n (θ 1 , . . . , θ n+1 ) → e −mR Y n (θ 1 , . . . , θ n+1 ), where with τ 0 ≡ 1, τ n+1 ≡ −1 and In order to ease the calculations, we introduce a compact notation for the evaluation of (n + 1)fold integrals with respect rapidities; we define where Ψ(θ 1 , . . . , θ n+1 ) is a function of the rapidities. The symbol − stands for the principal value of the integral. The need for the principal value follows from the fact that spin field matrix elements exhibit a kinematical pole [16]. The result of the integrations in (2.23) is, in general, a function of the rescaled coordinates {η j , τ j } which are indicated as a subscript. For the sake of simplicity, we will omit the subscripts when there is no ambiguity and we shall write Ψ in place of (2.23). The connected part of the correlation function is computed as follows As it has been illustrated in [16] and [30], a simple route for the calculation of spin field matrix elements is to take a first derivative with respect to the horizontal coordinate. Thanks to this procedure it is possible to get rid of the kinematical pole 1/(θ j − θ j+1 ) by taking the first derivative with respect to η j . Hence, by applying the differential operator ∂ x 1 · · · ∂ xn , we get rid of kinematical poles and we are left with The quantity indicated with 1 amounts to compute a (n + 1)-fold gaussian integral. Anticipating some results, it is convenient to write the outcome of the integration in the following way 1 η 1 ,τ 1 ;...;ηn,τn = λ n P n (x 1 , y 1 ; . . . ; x n , y n ) . (2.26) We will show in Sec. 4 that P n is the joint probability density of a Brownian bridge. This means that P n dx 1 · · · dx n is the probability for the interface -regarded as the trajectory of a Brownian particle -to pass through all intervals (x j , x j + dx j ) at y = y j , for j = 1, . . . , n. Leaving the details in Appendix A, the passage probability reads P n (x 1 , y 1 ; . . . ; x n , y n ) = 2 n/2 λ n n j=1 κ j Π n ( √ 2χ 1 , . . . , √ 2χ n |R 1...n ) .
(2.27) Some comments are in order. The dependence through the coordinates is encoded in the rescaled coordinates χ j and τ j , which are defined by for j = 1, . . . , n. Then, Π n indicates the multivariate normal distribution with correlation matrix R 1...n . Further details on the multivariate normal distribution are collected in Appendix A. The correlation matrix is a n × n symmetric matrix with 1 along the main diagonal and with entries in the upper triangle. The joint passage probability is normalized such that R n n j=1 dx j P n (x 1 , y 1 ; . . . ; x n , y n ) = 1 . (2.30) Marginal passage probabilities are obtained upon integration with respect to a subset of the n coordinates. For instance, (2.31) Note that P k (x 1 , y 1 ; . . . ; x k , y k ) is characterized by a k × k correlation matrix obtained by removing the last n − k rows and columns of R 1...n . Analogously, by integrating P n with respect to x k , we obtain a joint passage probability in the variables with labels 1, . . . , k − 1, k + 1, . . . , n and whose correlation matrix is obtained by removing the k th row and the k th column of R 1...n . Coming back to the correlation function, (2.25) reads (−∆ σ j ) P n (x 1 , y 1 ; . . . ; x n , y n ) . (2.32) The above equation is actually satisfied also by the "full" correlation function G n and not necessarily by its connected part. The reason is that, as we are going to show, the action of ∂ x 1 · · · ∂ xn on disconnected terms gives zero. Such a property follows from the fact that disconnected terms depend only on a subset of coordinates x 1 , . . . , x n . In order to make further progresses it is convenient to adopt a systematic notation for connected correlation functions. We introduce cluster functions of order n by means of the following integral representation −i θ j − θ j+1 η 1 ,τ 1 ;...;ηn,τn , (2.33) thanks to which it is possible to write the connected correlation function as follows with the half jump of vacuum expectation values across the ab interface. The product enclosed by square brackets in (2.34) contains the jumps of vacuum expectation values, which are intrinsically model-dependent quantities. Conversely, the cluster function G n is universal in the sense that it is shared by all models in which a single interface separates coexisting phases a and b. We further anticipate that, thanks to the normalization in (2.33), each cluster function tends to +1 when all the arguments are sent to +∞. By combining the above equations, (2.25) becomes In order to find the cluster function G n , we integrate back with respect to x 1 , . . . , x n . This procedure however must be followed carefully because (2.36) defines cluster functions up to arbitrary functions of a subset of coordinates. For instance, the cluster function G 2 would be determined up to functions of x 1 and x 2 . In order to fix the cluster function in a unique fashion, we impose a set of constraints which ensure the clustering property of correlation functions. The above discussion can be rephrased under a slightly different angle by using the identity which allows us to replace each simple pole in (2.33) with an integration with respect to an auxiliary variable conjugated to a rapidity difference. Such a variable can be identified by noting that x j , or its rescaled counterpart, η j , is coupled to the rapidity difference θ j − θ j+1 in the function Y n . By inserting (2.37) into (2.33) and carrying out the integrations with respect to the rapidities, we find (2.38) The lower integration extremes in the auxiliary variables η j are conventionally set to −∞. The residual term R n (x 1 , . . . , x n ), which satisfies ∂ x 1 · · · ∂ xn R n (x 1 , . . . , x n ) = 0 because of the property (2.36), will be identified in the following section.

Construction of n-body cluster functions
By expressing the passage probability P n in terms of the multivariate normal distribution Π n , (2.38) becomes We introduce the cumulative distribution function (CDF) of the multivariate normal distribution and denote it as follows (2.40) Thanks to (2.40), (2.39) becomes The functions R n can be fixed by requiring that G n satisfies the clustering property of correlation functions when at least one of its arguments is sent to infinity. The correct clustering is achieved provided that G n (x 1 , . . . , It is important to stress that such an information emerges from the analysis of the full correlation function G n , which includes both the connected part and the disconnected ones. Thus, in order to facilitate the exposition, we shall use the above input in order to construct the n-body cluster functions and we will check a posteriori that such a prescription is indeed the correct one. This consistency check is actually the content of Theorem 2.5, which will be enunciated at the end of Sec. 2.4. In order to illustrate the approach in a constructive fashion, we consider the simplest case, n = 1, which will guide our further considerations towards the case of arbitrary n. The function R 1 is inevitably a constant since it has to satisfy ∂ x 1 R 1 = 0. The value of such constant is determined by imposing G 1 (x 1 ) → ±1 for x → ±∞. Thus, the one-body cluster function is is the error function [43], and therefore G 1 (x 1 ) = erf(χ 1 ). The analysis of the case n = 2 reveals that Analogously to the case n = 1, the 2-body cluster function G 2 can be expressed in closed form by introducing a suitable set of special functions -Owen's T function [44,45]; we refer to [30] for a detailed account on this aspect. The results (2.42) and (2.43) already suggest what the mathematical structure of the n-body cluster function for arbitrary n should be. In order to construct a formal expression for G n , we need some preparatory definitions; we begin with the following one Definition 2.1 (block functions). Let p be an integer such that 0 p n. For p 1 the (p, n) block function B p,n is defined by the sum in the above runs over the set of ordered p-tuples with respect to the natural ordering (<) of integers. For p = 0, B 0,n = 1, then, B n,n = Φ n .
It is useful to write some explicit examples. For n = 2: . Note that the dependence on the correlation matrix, and so on the correlation coefficients, occurs for p 2.
The structure of the n-body cluster functions is formalized by the following theorem Theorem 2.2 (n-body cluster functions). The n-body cluster function is expressed in terms of block functions by means of The functions R n are identified by means of the next corollary Corollary 2.2.1. Since the term with p = 0 in (2.46) gives the first term in the right hand side of (2.41), it follows that R n is identified as the sum of the terms with p = 1, . . . , n in (2.46).
We stress that G n depends on both the horizontal and vertical coordinates, x j and y j . Such a dependence is codified by the rescaled coordinates χ j and by the rescaled vertical coordinates τ j , both introduced in (2.22). In particular, the dependence on τ j occurs also through the correlation coefficients ρ ij = (R 1...n ) ij ; see (2.29).
It is useful to introduce a graphical notation. We represent the CDF of the n-variate normal distribution by means of the following block diagram Thanks to the diagrammatic representation provided by (2.47) and (2.48), cluster functions admit the following graphical rewriting For mathematical convenience, we can take G 0 = G 0 ≡ 1 as the seed for the recursive hierarchy of cluster functions. The explicit form of G 3 reads (2.50) The correlation matrix R 123 is characterized by three independent correlation coefficients: ρ 12 , ρ 13 and ρ 23 ; the remaining one is obtained by virtue of the Markov property (see e.g. [46,47]): ρ 13 = ρ 12 ρ 23 . The function Φ 3 can be expressed in closed-form in terms of Steck's S and Owen's T functions 4 . Consequently, even for n = 3 it is possible to write the cluster function in an analytic form which involves single integrals instead of three-fold ones [48]. Carrying on the above procedure, we can write the cluster function for the four-point correlation function, namely (2.51) or equivalently, within the pictorial form, we have (2.52) We observe that for n = 4 the independent correlation coefficients are ρ 12 , ρ 23 and ρ 34 . The other correlation coefficients follow from the Markov property: ρ 13 = ρ 12 ρ 23 , ρ 14 = ρ 12 ρ 23 ρ 34 , and ρ 24 = ρ 23 ρ 34 .
The next task we need to carry out is to prove Theorem 2.2. To this end, we need to recall the asymptotic properties satisfied by block functions. Lemma 2.3. The block function B p,n vanishes when at least one of its argument is sent to −∞, e.g., lim For x n → +∞, the block function B 1,n satisfies the following property , The limits in which x j → ±∞ with j = n are treated along the same lines by replacing R 1...n−1 with the correlation matrix R 1...ĵ...n , where R 1...ĵ...n is obtained by removing the j-th row and the j-th column of R 1...n , andĵ stands for the removed label.
Proof. The derivation of the asymptotic relations listed in Lemma 2.3 follows by using elementary properties of cumulative distribution functions. Let us consider the first of the properties listed in Lemma 2.3. The limit in which x n → −∞ is established thanks to because the lower integration extrema in the CDFs are −∞. Let us consider the case p = n.
The limit x n → +∞ can be analyzed by using the identity which is a natural consequence of the marginalization property of the probability distribution P n ; see (2.31). The properties with 1 p n − 1 follow straightforwardly.
Theorem 2.4 (clustering). The n-body cluster function satisfies the following clustering properties lim The proof of Theorem 2.4 follows as a straightforward application of Lemma 2.3.
As an explicit illustration of the asymptotic properties, we consider the asymptotic properties of the block functions B p,3 , which constitute the building blocks of the three-point correlation function. For p = 3, we have lim which is actually the asymptotic property of the cumulative distribution function Φ 3 . Then, for p = 2 one finds lim (2.58) We conclude this section by adding some considerations about the construction of block functions. According to Definition 2.1 the (p, n) block function is constructed by summing cumulative functions Φ p in which the p arguments are ordered p-tuples drawn from of the set A n = {1, . . . , n}. This observation allows us to rationalize the construction of block functions by putting them in touch with Hasse diagrams [49,50] in discrete mathematics.
In order to proceed along this direction, we recall the definition of power set. The power set of A n , denoted P(A n ), is the set which contains all subsets of Subsets of A n are naturally ordered by set inclusion (⊆). The partially ordered set (poset) (P(A n ), ⊆) can be visualized by means of a graph in which the largest element is placed at the top, the smallest at the bottom, and other elements are allocated in between according; the notion of large/small has to be interpreted in terms of the cardinality. Two vertices are connected by an edge if the elements are ordered by set inclusion (⊆) [49]. Coming back to the example quoted one moment ago, the Hasse diagram for the power set of A 2 and A 3 are shown in Fig. 2a and Fig. 2b, respectively.
Hasse diagrams H n comprise n+1 levels in which elements share the same cardinality. Levels with cardinality p are indicated with a dash-dotted notation in Fig. 2. It then follows a level with cardinality p contains n p elements. Then, the level p contains the elements which define the block function B p,n . The asymptotic clustering properties listed in Lemma 2.3 can be viewed in terms of Hasse diagrams. For instance, the limit x 3 → ∞ amounts to remove the node with label 3 and those bonds attached to it in the diagram H 3 . Hasse diagrams can be used also to classify the disconnected diagrams arising from disconnected parts of matrix elements, as will be clear in the next section.

Disconnected parts
We can now illustrate how to compute the contribution of disconnected matrix elements to the n-point correlation function of the spin field. The disconnected parts of matrix elements which appear in the right hand side of (2.19) are constructed by contracting legs according to the procedure outlined for n = 1 in (2.16). The decomposition of matrix elements given in (2.19) is then written as follows where the sum runs over the number m of disconnected spin fields. Thus, D ab stands for the fully connected matrix element corresponding to the diagram in the right hand side of (2.19) and whose contribution yields the connected part of the n-point correlation function. Conversely, ab is the fully disconnected matrix element. The decomposition (2.59) of matrix elements induces an analogous expansion of the n-point correlation function, which reads  (2.19) are replaced by disconnected lines as those shown in the second term on the right hand side of (2.16). Correspondingly, matrix elements which contribute to D (m) ab contain the product of m Dirac deltas of type δ(θ j − θ j+1 ) and the product of n − m two-particle form factors F σ j aba (θ j − θ j+1 + iπ). Let us consider some illustrative examples. The diagrams which appear in the disconnected matrix elements D (1) ab admit the graphical depiction shown in (2.61). The term "perm" in (2.61) indicates the sum over permutations of diagrams formed by detaching -one at the time -the spin field with label j = 1, . . . , n − 1.
The graphical construction which gives the disconnected matrix elements formed by detaching two spin fields proceeds in an analogous manner. Thus, D (2) ab admits the graphical decomposition shown in (2.62). +perm.
(2.62) By applying the above rules it is straightforward to construct diagrams corresponding to disconnected matrix elements with an arbitrary number of disconnected spin fields.
We can now address the calculation of (G n (x 1 , . . . , x n )) [D (m) ab ] . In order to simplify the exposition, we show the calculation for a particular type of disconnected diagrams, denoted D (m) ab , which are obtained by disconnecting the last m spin fields with labels j = n−m+1, . . . , n. The diagrams which contribute toD (1) ab are precisely those depicted in the right hand side of (2.61). Analogously, the diagrams which contribute toD (2) ab are those depicted in the right hand side of (2.62). There is actually no loss of generality in this choice since an arbitrary disconnected diagram can be obtained by a permutation of the spin and coordinate labels. Therefore, we shall focus on the following matrix element The anatomy of (2.63) follows by noting that the first product is originated by those spin fields which form the connected part of the diagram, while the second product stems by tying together those legs which are detached from the disconnected spin fields. The occurrence of σ j follows since the arithmetic average between diagrams obtained within the left and right annihilations has to be performed [30].
By keeping the leading-order term in the small rapidity expansion, we can write the matrix element in the factorized formD  It is interesting to observe how the pictorial representation of matrix elements provides insights on the structure of G n . A matrix element represented by a diagram in which the spin fields σ 1 . . . , σ n−m are connected yields a cluster function G n−m (x 1 , . . . , x n−m ) which depends on the spatial coordinates carried by those spin fields which constitute the connected part of the diagram. The spatial coordinates relative to disconnected spin fields do not report in the resulting cluster function.
It is now evident how to construct all the disconnected diagrams belonging to the family ab . Firstly, we observe that the number of such diagram is #D (m) ab = 2 m n m because each spin field can be disconnected either passing left or right aside it, hence the factor 2 m follows. Thus, up to left/right combinatorics, there is one ( n 0 ) fully connected diagram, there are n 1 = n diagrams with one disconnected spin field, and so on. Note that, the total number of diagrams is n m=0 n m = 2 n . It is thus clear how Hasse diagrams can be used in order to classify the disconnected diagrams too.

Full result and specific cases
We are now in the position to construct the full correlation function G n . An arbitrary disconnected diagram can be identified by specifying the labels of those spin fields which are disconnected. Let V m be the ordered set of vertices which composes the connected part of the diagram with m disconnected spin fields and let us denote its vertices with the labels j 1 , . . . , j n−m . Analogously, let V m be the set of vertices which compose the disconnected part. Clearly, for any m, . . , n}. The contribution stemming from the diagrams with m disconnected spin fields reads we recall that σ j is given by (2.35). Thanks to (2.60), the n-point spin correlator G n is given by It is instructive to consider some examples. The case n = 1 gives the magnetization profile (2.72) The above agrees with the result of [16]. In the last line follows by using the expression of the one-body cluster function G 1 given below (2.42). We also observe that (2.72) retrieves the known magnetization profile for the Ising model ; the correction at order R −1/2 vanishes (see (4.9)). Let us consider the case n = 2 corresponding to the pair correlation function of the order parameter. The connected part yields σ 1 σ 2 G 2 (x 1 , x 2 ). The disconnected parts with m = 1 give − σ 1 σ 2 G 1 (x 2 )− σ 2 σ 1 G 1 (x 1 ) and the fully disconnected part contributes with σ 1 σ 2 . Collecting the various pieces, we find 73) which perfectly matches with the findings of [30]. As a further example, the three-point correlation function is given by

76)
where σ j (x j ) stands for the omission of σ j (x j ) in the correlation function.
The proof follows by using the results of the (clustering) Theorem 2.4.

The limit R → ∞
In this case, τ j = 2y y /R → 0, meaning that all correlation coefficients tend to unity, i.e., Correspondingly, the correlation matrix reduces to a n × n matrix whose entries consists of all 1s; we denote such a matrix with J n . Analogously, the limit R → ∞ projects the variables χ j , which encode the dependence through the coordinates x j and y j , to the origin, i.e, lim R→∞ χ j = 0 . (2.78) In both the limits (2.77) and (2.78) the coordinates x j , y j with j = 1, . . . , n are fixed. The limit (2.78) implies that the cumulative distribution functions which appear in the block functions are evaluated at the origin. Said differently, the cumulative distribution functions become the so-called orthant probabilities [51,52]. The calculation of orthant probabilities is a notoriously difficult problem. The first few orthant probabilities are:ˆ0 (2.79) The above result may indicate a general pattern for the orthant of the n-variate normal distribution. This is actually not the case, as it is revealed by the orthant of the quadrivariate normal distribution [53].
The limit (2.77) comes in our rescue. Although the orthant of the n-variate normal distribution is a complicated function of the correlation coefficients ρ ij , the orthant drastically simplifies when ρ ij = 1, and the corresponding result is simply 1 /2. With this in mind, the cumulative function given in (2.47) reduces to 1 /2 and the block function (2.48) reduces to 1 /2 times the number of elements in the level p 1 of the Hasse diagram H n , therefore  (2.83) The interface separating phases a and b is characterized by midpoint fluctuations of order R 1/2 along the x-axis. For R → ∞ the unbounded interfacial fluctuations yield the averaging over the phases a and b given by (2.83). For the Ising model, the averaging property (2.83) is known from rigorous result [9]. The above derivation shows that the averaging (2.83) is actually a more general feature.

Large-R expansion
In this section, we examine the correlation function G n including the leading-order corrections in finite size. By extending the approach of [16,30] the correction at order R −1/2 is interpreted in terms of a probabilistic picture.

Correction at order R −1/2 : connected part
To be definite, we begin by examining the connected part of G n . The treatment of finite-size corrections for disconnected parts will be facilitated by the treatment of the connected part. The large-R expansion of G n can be studied in a systematic fashion by expanding the numerator of (2.2), as well as the partition function Z ab (R), in powers of the small parameter R −1/2 . The large-R expansion of the numerator in (2.2) proceeds by Taylor expanding the integrand in (2.12) at small rapidities. Retaining the connected part, the integrand in (2.12) reads The function C(θ 1 , . . . , θ n+1 ) is expanded at small rapidities and the corresponding result is organized as follows where C ∆ is a homogeneous function of order ∆, i.e., C ∆ (αθ 1 , . . . , αθ n+1 ) = α ∆ C ∆ (θ 1 , . . . , θ n+1 ), α > 0. It is simple to check that terms in the aforementioned series with homogeneity exponent ∆ contributes to the correlation function G n (x 1 , . . . , x n ) at order R −(∆+n)/2 . The leading order term in the large-R expansion is thus generated by the function C −n and its corresponding expression is provided in (2.20). The function C ∆ with ∆ = −n + 1 gives the first subleading correction which occurs at order R −1/2 . Regarding the large-R expansion of the partition function, we denote the leading order expression (2.5) with Z (0) ab (R). Therefore The correction at order O(R −1/2 ) is absent since the low-energy expansion of f ab (θ) does not exhibit odd powers of θ. As a result, the first subleading correction for the n-point correlation function is, in general, proportional to O(R −1/2 ) and is entirely originated by the matrix element C −n+1 . The latter is obtained by multiplying n−1 kinematical poles with one of the c (σ k ) ab factors appearing in the low-energy expansion (2.17) and then summing over permutations of labels. We have with I k,n an overall factor which depends on the vacuum expectation values, and D k,n the following function of the rapidities Note how D k,n can be obtained from the connected matrix element simply by removing one annihilation pole. The removal of the pole 1/(θ j − θ j+1 ) in the matrix element M n can be achieved thanks to a differentiation with respect to the horizontal coordinate x j conjugated to the rapidity difference θ j − θ j+1 ; hence, The last equality brings in touch the first subleading correction to connected matrix elements with the cluster functions introduced in Sec. 2.
We are now in the position to write the correction at order R −1/2 . The large-R expansion of the n-point connected correlation function can be written as follows with G CP n ℓ = O(R −ℓ/2 ). The term with ℓ = 0 is the one computed in Sec. 2. According to the above discussion, the term with ℓ = 1 reads

Correction at order R −1/2 : disconnected parts and full result
Once we have established how the large-R expansion is implemented for the connected part, the analysis of the disconnected parts follows from the diagrammatic construction of matrix elements. Let us consider n = 1 as the first example. In this case the disconnected term coincides with the fully disconnected one, the latter simplifies with the partition function up to the factor σ 1 . The correction is thus entirely due to the connected part. As a result, the expansion of the magnetization profile reads with the leading-order (ℓ = 0) connected part given by According to (3.10), the first subleading correction is given by The above perfectly matches the results established in [16,17]. Let us consider the case n = 2. The diagrams with the field σ 1 disconnected contribute to G 2 with σ 1 G CP 1 (x 2 ) 1 and analogously when the labels 1 and 2 are interchanged. The large-R expansion of the pair correlation function reads where G 2 (x 1 , x 2 ) 0 is the expression given in (2.73). The terms due to disconnected matrix elements are those multiplied by factors σ j . Subleading corrections for n = 2 are computed from (3.10) and expressed as follows This result can be written in a more explicit way by carrying first derivatives with respect to x 1 and x 2 of the cluster function G 2 . It is simple to show that (3. 16) Grouping together corrections at order R −1/2 stemming from both connected and disconnected parts, we find (3.17) The above expressions coincide with the results given in [30]. We further stress how the clustering of the two-point correlation function is satisfied at order R −1/2 . This can be easily inspected by considering the following limits lim 18) which are in agreement with (3.13). By carrying on the above procedure, the three-point correlation function expands as follows (3.19) with the leading-order term G 3 (x 1 , x 2 , x 3 ) 0 given by (2.74) and the terms at order R −1/2 given by (3.10). The result for arbitrary n follows along the same lines.

Probabilistic interpretation
It is possible to reconstruct the n-point correlation function within a probabilistic interpretation in which the interface is regarded as a fluctuating line with fixed extremities. Let P n (x 1 , y 1 ; . . . ; x n , y n )dx 1 · · · dx n (3.20) be the probability that the interface crosses the intervals (x j , x j + dx j ) at time y j . Then, let be the magnetization profile at in the point x j and u j the abscissa in which the interface crosses the horizontal line y = y j . The first two terms in the right hand side of (3.21) account for coexisting phases sharply separated by a structureless interface. Endowing the interface with interface structure amounts to the subsequent terms beyond the sharp picture, as indicated in (3.21). The sum over interfacial configurations which define the n-point correlation function is formulated as follows G n (x 1 , . . . , x n ) =ˆR n du 1 . . . du n P n (u 1 , y 1 ; . . . ; u n , y n ) n j=1 σ j (x j |u j ) . (3.22) The fact that P n occurring in (3.22) is the expression found in field theory can be established by matching (3.22) with the field-theoretical calculation for arbitrary n. This is what we will do in the following. We begin by focusing on the leading order in the large-R expansion which is captured by the first two terms in (3.21). The development of the product appearing in (3.22) yields 2 n terms whose integral with respect to u 1 , . . . , u n reproduces the cluster functions introduced in Sec. 2. Proving the above statement is a simple matter. We denote the n-fold integral with respect to horizontal coordinates {u j } j=1,...,n with the compact notation f (u 1 , . . . , u n ) ≡ˆR n du 1 . . . du n f (u 1 , . . . , u n )P n (u 1 , y 1 ; . . . ; u n , y n ) . (3.23) Then, we employ the following abbreviation for the sign function s j ≡ sign(x j − u j ) = 2ϑ(x j − u j ) − 1 and ϑ j ≡ ϑ(x j − u j ), with ϑ(· · · ) Heaviside theta function. Block functions are expressed as follows

(3.24)
Consequently, cluster functions admit the following representation For m = 0 the normalization condition gives G 0 = 1 = 1, as we also stipulated Sec. 2. The matching between the field theoretic calculation and the probabilistic interpretation is thus completely characterized at the leading order in the large-R expansion.
We can now establish the matching at order R −1/2 . In order to do this, we focus on the connected part of the correlation function, [G CP n (x 1 , . . . , x n )] 1 . The latter can be extracted from the probabilistic picture (3.22) by removing the offset values σ j in (3.21), namely (3.27) By applying the multiple derivative ∂ x 1 · · · ∂ xn , one finds y 1 ; . . . ; x n , y n ) .

(3.29)
By matching (3.28) and (3.29), we readily extract the structure amplitudes As consistency requires, (3.30) agrees with calculations based on n = 1 and n = 2, respectively in [16] and [30]. It is actually possible to cast the above results within the diagrammatic framework which we employed in the previous section. We note the following property which expresses the cluster function in terms of the matrix element generated by the product of n kinematical poles. Then, the cluster function G n admits the following diagrammatic representation in terms of block diagrams In turn, the relationships (3.31) and (3.32) allow us to connect the calculation of the (n + 1)-fold integral with respect to rapidities M n to a diagrammatic expansion. Such a relationship turns out to be extendable to the instance in which one of the simple poles of M n is replaced by 1, which is precisely the construction which leads to the matrix element D k,n . Note that each circle appearing in the block diagrams is in one-to-one correspondence with a Heaviside theta function, as we have shown.
The diagrammatic representation of block function is extended by introducing the modified diagrams in which one of the Heaviside theta is replaced by a Dirac delta and the latter is represented with an orange diamond; thus where δ j ≡ δ(x j − u j ). The second equalities in (3.33) follow by virtue of ∂ x j ϑ j = δ j . As an example, the connected pair correlation function at order R −1/2 reads (3.34) The diagrams appearing in the above can be easily computed. Focusing on the first two diagrams, the results are: and analogous results are obtained for the other diagrams. It is thus evident that (3.34) can be written in the following explicit form this result perfectly matches with the connected part of the expression (3.17) obtained from the field theoretical calculation. Such a connected part can be selected simply by removing the terms proportional to the offsets σ 1 and σ 2 .

Triplet correlations
As an application of the formal results derived in the previous sections, here we consider the explicit form of the three-point correlation function The integrals which define the corresponding cluster functions can be evaluated in closed form.
We leave the technical calculations in Appendix B and here we recall the main notations. Thus, we introduce the rescaled coordinates η = x/λ, τ = 2y/R, the correlation coefficient 38) and the parameter r(τ ) = ρ 12 (τ )/ 1 − ρ 2 12 (τ ). The correlation function (3.37) reads where Y(η, τ ) and K(η, τ ) are the functions given by and K(η, τ ) = 2 π sin −1 (ρ 2 12 ) + 8T ( √ 2η, r(τ )) , (3.41) where T is Owen's function [45] T ( √ 2η, r) = 1 2πˆr 0 du e −(1+u 2 )η 2 1 + u 2 . (3.42) Let us comment on some general properties. In the limit x → ±∞ the triplet correlation function reduces to the pair correlation function with spin fields placed along the interface and symmetrically displaced with respect to the horizontal axis. In that limit, one finds lim x→∓∞ σ(0, y)σ(x, 0)σ(0, −y) ab = σ a[b] σ(0, y)σ(0, −y) ab ; (3.43) the quantity in the right hand side is the two-point correlation function with spin fields along the interface Interestingly enough, for τ = 1 /3 -corresponding to y = R/6, and r = 1 -the special functions in (3.39) reduce to powers of the error function, in particular Lastly, the correlation function with the three spins placed along the line which joins the pinning points reads where [g 3 (0, y)] 1 = O(R −1/2 ) is the subdominant correction due to interface structure effects. The result for [g 3 (0, y)] 1 can be obtained by taking y 1 = y, y 2 = 0 and y 3 = −y in thee expression for the correction at order R −1/2 of the correlation function σ c (0, y 1 )σ c (0, y 2 )σ c (0, y 3 ) ab , which is calculated in Appendix B.2. By taking the above limiting case in (B.10), we obtain the subdominant correction The occurrence of long-range interfacial correlations can be verified by expanding the correlation coefficient ρ 12 and κ for small τ . Focusing on the leading-order tern, we find the asymptotic behavior with ξ b ≪ y ≪ R. An analogous expansion can be performed for the interface structure correction given by (3.47).
The term proportional to √ y in (3.48) is the signature of long range interfacial correlations.
This power-law behavior in the direction parallel to the interface has to be compared with the exponential decay of correlations which characterizes the transverse direction. This feature can be neatly appreciated simply by evaluating the derivative of g 3 (x, y) with respect to x; a simple calculation gives The exponential factor e −η 2 produces the confinement within the interfacial region of the longrange fluctuations of the order parameter. Let us us consider now the Ising model. Denoting the spontaneous magnetization with M = σ + , the leading-order form of the triplet correlation becomes (3.50) The above vanishes for x = 0. However, going away from x = 0 the decay of correlations along the direction parallel to the interface exhibits a long-range character analogous to (3.48) as well as an additional dependence on x with the anisotropic features (i.e., dependence on x and y) discussed above. A detailed asymptotic analysis of (3.50) and the comparison with results obtained with Monte Carlo simulations is carried out in a forthcoming publication [48]. Finally, we discuss the case of the q-state Potts model [54]. For ferromagnetic interactions and with q 4 the model exhibits a continuous phase transition [55]. In the low-temperature phase there are q degenerate ground states and, in the scaling limit, phase separation between them is described by field theory [16]. Thanks to permutational symmetry, the vacuum expectation values of the order parameter field satisfy (3.52) We observe that when c equals one of the two boundary colors, e.g. c = b (with c = a), the correlation function (3.52) reduces to a particularly simple expression It has to be observed how the result corresponding to the Ising model given in (3.50) is retrieved in the limit q → 2.
For the q-state Potts model with q = 3 and q = 4 it is possible to consider the correlations between non-boundary colors. By taking c = a, b the first three terms in the right hand side of (3.52) vanish and one finds a term proportional to 1/(q − 1) 3 up to corrections proportional to R −1/2 . This feature is actually expected because the non-boundary color contributes in a nontrivial way to the magnetization profile at order R −1/2 , the same happens for correlation functions. As we are going to show in an explicit fashion, the term proportional to R −1/2 depends on y. We can compute the correction at order R −1/2 in (3.52) by adopting the probabilistic interpretation illustrated in Sec. 3.3. For the sake of simplicity we show the specific form of these corrections for the special case x = 0. The subleading correction is given by (3.47) with the structure amplitude for the q-state Potts field theory given by where B 3 = 1/(2 √ 3) and B 4 = 2/(3 √ 3) [16]. By inserting (3.54) into (3.47), we find the subleading correction (3.55) By performing a small-y expansion it is possible to show that (3.55) exhibits power-law correlations which are analogous to those obtained at the leading order in (3.48). For q = 3 and q = 4 the correlations of the non-boundary color are characterized by a non-vanishing amplitude. For q = 2 the color c must coincide either with a or b and the amplitude vanishes. In this case, we expect the first correction to occur at order R −1 . All these features are actually shared by the interface structure correction of the magnetization profile and can be interpreted as the formation of isolated droplets of phase c adsorbed along the ab interface [16].

Corrections at order R −1
We have seen that finite-size corrections proportional to R −1/2 computed within field theory match with a calculation based on the probabilistic interpretation. We show in this section that corrections at order R −1 for the magnetization profile can be interpreted within the probabilistic picture by allowing certain structure amplitudes to be y-dependent. However, by using such an information gained for n = 1, the case n = 2 does not necessarily lead to a matching between the two formulations. We will cover these aspects by focusing on the explicit example of the Ising model.

Two-point correlation function
The calculation proceeds as follows with I 2n−2 a homogeneous function of order 2n − 2. The expansion yields . (4.13) The spin-spin correlation function expands as follows , (4.14) thus In order to simplify the analysis, we restrict ourselves to the parallel correlation function, G 2 | ≡ G 2 (x, y; x, −y). As a further simplification for our considerations, we take the double derivative ∂ x 1 ∂ x 2 and evaluate it for spin field in the parallel arrangement defined above. Therefore, we examine the subscript means that x 1 = x 2 ≡ x, y 1 = −y 2 ≡ y are set afterwards the application of ∂ η 1 ∂ η 2 . Now, we look more closely to ∂ 2 η 1 η 2 I n , for n = −2 and n = 0. We have with 1 = λ 2 P 2 (x 1 , y 1 ; x 2 , y 2 ) and The second term can be written as follows The above admits a simple expression for spin fields in parallel arrangement, i.e., with x 1 = x 2 ≡ x, y 1 = −y 2 ≡ y, in particular with The other term is arranged in a similar way. By writing Summing up all the pieces, we obtain (4.24) The factor λ −2 in the O-symbol is due to the differentiation with respect to x 1 = λη 1 and x 2 = λη 2 . The calculation within the probabilistic interpretation follows straightforwardly and reads 25) the superscript "prob." stresses that such an expression has been derived within the probabilistic interpretation. The second term can be written in a form which is similar to (4.24) (4.26) We see that the field-theoretic calculation (4.24) and the one carried out within the probabilistic interpretation, (4.26), agree at the leading order but disagree at the order R −1 . The disagreement is actually caused by the term proportional to f 2 and is ultimately originated by the features of the boundary field theory. On the other hand, the term proportional to the factor 2 /3, which emerges from the low-energy properties of the bulk form factor, does not originate the disagreement.
The above analysis suffices in order to provide an example which exhibits an explicit breakdown of the matching between probabilistic interpretation and field theory. As we have already proved, the probabilistic approach has to be limited to corrections at order R −1/2 in which the low-energy behavior of the boundary features does not report.

Stress tensor trace and mixed correlation functions
The calculation method illustrated in the previous sections can be applied to mixed correlation functions involving both the trace of the stress tensor and spin fields. In this conclusive section, we give an account on this aspect. Denoting the stress tensor trace with Θ(x), the counterpart of the matrix element decomposition (2.16) reads 27) or equivalently, Contrary to the spin field, for the stress tensor Θ a = Θ b ≡ Θ ; thus, its two-particle form factor can be expanded as follows The normalization of Θ actually implies F Θ aba (iπ) = 2πm 2 [56]. Let us consider the n-point correlation function of the stress tensor trace, Θ(x 1 ) · · · Θ(x n ) ab . The connected part follows straightforwardly Θ(x 1 ) · · · Θ(x n ) CP ab = F Θ aba (iπ) n m n P n (x 1 , y 1 ; . . . x n , y n ) , (4.30) up to subdominant large-R corrections. Therefore, the joint n-intervals passage probability is proportional to the connected part of the n-point correlation function of the stress tensor. It has to be observed that (4.30) scales as R −n/2 . As a first consequence of (4.29), the leading term in the large-R expansion is the one which counts the maximum number of disconnected pieces. The fully disconnected term cancels exactly the partition function at the denominator and yields a spatially-independent offset given by Θ n . The next-to-leading term comes at order R −1/2 and it is due to the contraction of n − 1 disconnected pieces with one connected matrix element; such terms are captured by the matrix element The corresponding result for the n-point correlations of Θ reads We can now consider a mixed correlation function which involves n − m spin fields. The leading-order term follows by contracting the product of (n − m) kinematical poles with m Dirac deltas stemming from the stress tensor matrix element (4.28). The corresponding result reads (4.33) As a consistency check we consider two limiting cases. For m = 0 the above reduces to the connected n-point spin correlation function and the corresponding result (2.34) is found as a limiting case. For m = n we retrieve the connected n-point correlator of Θ given by the first term in the right hand side of (4.32).

Conclusions
In this paper, we considered the scaling limit of a generic two-dimensional ferromagnetic system at phase coexistence near a second order phase transition point. We showed how field theory provides exact results for n-point spin and stress tensor trace correlation functions in presence of a fluctuating interface. More specifically, the system we considered is defined on an infinite strip of width R much larger than the bulk correlation length. Boundary conditions are used in order to enforce phase separation through an interface which spans between the two edges and whose endpoints are pinned.
By extending the field-theoretical technique developed for one-and two-point correlation functions, respectively in [16] and [30], we have been able to find the exact analytic form of order parameter and stress tensor correlation functions. Analogously to the case of the twopoint correlation function, we have showed that, as long as R is finite, the n-point correlation function is characterized by long-range correlations in the direction parallel to the interface. The spatial extent of the interface midpoint fluctuations grows as R 1/2 and, for R = ∞, these unbounded fluctuations lead to an exponential decay of bulk correlations averaged over the two coexisting phases separated by the interface.
More technically, these results follow by exploiting general low-energy properties of twodimensional field theory whose excitations -in two dimensions -are topological (kink) particles. We found that the leading asymptotic form of correlation functions involving spin fields is completely codified by the kinematical pole singularity exhibited by matrix element of the order parameter field.
Among our findings, the dominant asymptotic form of n-point correlation functions is expressed in terms of n-body cluster functions which are constructed out of cumulative distribution functions of the n-variate gaussian distribution. The first subleading finite-size correction, which is proportional to R −1/2 and arises from effects due to interface structure, depends on the bulk universality class only. Specificities related to the boundaries, which are incorporated in the low-energy behavior of matrix elements of boundary changing operators, do not report at order R −1/2 , but appear at order R −1 . Both the leading term and the first subleading corrections can be interpreted within a probabilistic picture in which the interface is regarded as the worldline of a particle which propagates randomly between the pinning points by undergoing a Brownian bridge. By using the Ising model as a specific example, we also show that the subleading correction at order R −1 does not necessarily emerge from the probabilistic description. We identify the origin of the mismatch as a specificity arising from matrix elements of the boundary condition changing operators.
Throughout this manuscript we have introduced a diagrammatic notation (block diagrams) for n-body cluster functions which facilitates the handling of expressions at both the leading and first subleading orders. Such a notation proved to be useful also in establishing a graphical connection between disconnected matrix elements and their contribution to the correlation function.
We conclude by discussing some interesting perspectives. The reconstruction of n-point correlation function through the probabilistic interpretation, which we have shown to be correct at both the leading order and including corrections at order R −1/2 , can be used in order to find exact results in closed form once the passage probability is known. This is indeed the case for n = 3 [48] and n = 4 [57] in which numerical simulations confirm the analytic results. The extension of the techniques developed in this paper has been merged with the techniques of [18] in a companion paper for the study of correlations in the half-plane. There, explicit results for the spin-spin correlation function on the half-plane with boundary conditions enforcing a droplet have been found and successfully tested by means of high-precision Monte Carlo simulations [58].
where D is the diffusion coefficient, (x 0 , t 0 ) defines the initial state and (x 1 , t 1 ) the final one. Since (A.1) is a probability density,´R dx 1 W (x 1 , t 1 |x 0 , t 0 ) = 1. Let I j = (x j , x j + dx j ) be a space interval at time t j as shown in Fig. 3.