text:
# Singular CombinatoricsWork supported in part by the IST Programme of the EU under contract number IST-1999-14186 (ALCOM-FT).

## 1. Introduction

Large random combinatorial structures tend to exhibit great statistical regularity. For instance, an overwhelming proportion of the graphs of a given large size are connected, and a fixed pattern is almost surely contained in a long random string, with its number of occurrences satisfying central and local limit laws. The objects considered (typically, words, trees, graphs, or permutations) are given by construction rules of the kind classically studied by combinatorial analysts via *generating functions* (abbreviated as gfs). A fundamental problem is then to extract asymptotic information on coefficients of a gf either explicitly given by a formula or implicitly determined by a functional equation. In the univariate case, asymptotic counting estimates are derived; in the multivariate case, moments and limit probability laws of characteristic parameters will be obtained.

In what follows, given a combinatorial class $`𝒞`$, we let $`C_n`$ denote the number of objects in $`𝒞`$ of size $`n`$ and introduce the *ordinary* and *exponential* gf (ogf, egf),

$$\text{ogf}\text{:}C(z):=\underset{n0}{}C_nz^n,\text{egf}\text{:}\widehat{C}(z):=\underset{n0}{}C_n\frac{z^n}{n!}.$$

Generally, egfs and ogfs serve for the enumeration of labelled classes (atoms composing objects are distinguished by labels) and unlabelled classes, respectively. One writes $`C_n=[z^n]C(z)=n![z^n]\widehat{C}(z)`$, with $`[z^n]()`$ the coefficient extractor.

General rules for deriving gfs from combinatorial specifications have been widely developed by various schools starting from the 1970’s and these lie at the heart of contemporary combinatorial analysis. They are excellently surveyed in books of Foata & Schützenberger (1970), Comtet (1974), Goulden & Jackson (1983), Stanley (1986, 1998), Bergeron, Labelle & Leroux (1998). We shall retain here the following simplified scheme relating combinatorial constructions and operations over gfs:

Construction  Labelled case Unlabelled case    Disjoint union $`+𝒢`$ $`f\left(z\right)+g\left(z\right)`$ $`\widehat{f}\left(z\right)+\widehat{g}\left(z\right)`$  Product $`\times 𝒢`$, $`𝒢`$ $`f\left(z\right)g\left(z\right)`$ $`\widehat{f}\left(z\right)\widehat{g}\left(z\right)`$  Sequence $`𝔖\left\{\right\}`$ $`(1f\left(z\right)))^1`$ $`(1f\left(z\right)))^1`$  Set $`𝔓\left\{\right\}`$ $`\mathrm{exp}(f\left(z\right)`$) $`\mathrm{exp}\left(f\left(z\right)+\frac{1}{2}f\left(z^2\right)+\mathrm{}\right)`$  Cycle $`\left\{\right\}`$ $`\mathrm{log}(1f\left(z\right))^1`$ $`\mathrm{log}(1f\left(z\right))^1+\mathrm{}`$    (1.1)
Such operations on gfs yield a wide variety of analytic functions, either given explicitly or as solutions to functional equations in the case of recursively defined classes. It is precisely the goal of singularity analysis to provide means for extracting asymptotic informations. What we termed “singular combinatorics” aims at relating combinatorial form and asymptotic-probabilistic form by exploiting *complex-analytic* properties of generating functions. Classical approaches  are Tauberian theory and Darboux’s method, an offspring of elementary Fourier analysis largely developed by Pólya for his programme of combinatorial chemistry . The path followed here, called “singularity analysis” after , consists in developing a systematic correspondence between the local behaviour of a function near its singularities and the asymptotic form of its coefficients. (An excellent survey of central aspects of the theory is offered by Odlyzko in .)

## 2. Basic singularity analysis

Perhaps the simplest coefficient estimate is $`[z^n](1z)^\alpha n^{\alpha 1}/\mathrm{\Gamma }(\alpha )`$, a consequence of the binomial expansion and Stirling’s formula. For the basic scale

$$\sigma _{\alpha ,\beta }(z)=(1z)^\alpha \left(\frac{1}{z}\mathrm{log}(1z)^1\right)^\beta ,$$

much more is available and one has a fundamental translation mechanism :

###### Theorem 1 (Coefficients of the basic scale)

For $`\alpha _0`$, $`\beta `$, one has

$$[z^n]\sigma _{\alpha ,\beta }(z)\underset{n\mathrm{}}{}\frac{n^{\alpha 1}}{\mathrm{\Gamma }(\alpha )}\left(\mathrm{log}n\right)^\beta .$$
(2.1)
Proof. The estimate is derived starting from Cauchy’s coefficient formula,

$$[z^n]f(z)\frac{1}{2i\pi }_\gamma f(z)\frac{dz}{z^{n+1}},$$

instantiated with $`f=\sigma _{\alpha ,\beta }`$. The idea is then to select for $`\gamma `$ a contour $``$ that is of Hankel type and follows the half-line $`(1,+\mathrm{})`$ at distance exactly $`1/n`$ (Fig. 1). This superficially resembles a saddle-point contour, but with the integral normalizing to Hankel’s representation of the Gamma function, hence the factor $`\mathrm{\Gamma }(\alpha )^1`$. $`\mathrm{}`$

The method of proof is very flexible: it applies for instance to iterated logarithmic terms ($`\mathrm{log}\mathrm{log}`$) while providing full asymptotic expansions; see  for details.

A remarkable fact, illustrated by Theorem 1, is that larger functions near the singularity $`z=1`$ give rise to larger coefficients as $`n\mathrm{}`$. This is a general phenomenon under some suitable auxiliary conditions, expressed here in terms of *analytic continuation*: a $`\mathrm{\Delta }`$-domain is an indented disc defined by ($`r>1`$, $`\vartheta <\pi /2`$)

$$\mathrm{\Delta }(\vartheta ,r):=\left\{z\right||z|<r,\vartheta <\mathrm{Arg}(z1)<2\pi \vartheta ,z1\}.$$

###### Theorem 2 ($`O`$-transfer)

With $`f(z)`$ continuable to a $`\mathrm{\Delta }`$-domain and $`\alpha _0`$:

$$f(z)\underset{z1,z\mathrm{\Delta }}{=}O\left(\sigma _{\alpha ,\beta }(z)\right)[z^n]f(z)=O\left([z^n]\sigma _{\alpha ,\beta }(z)\right).$$

Proof. In Cauchy’s coefficient formula, adopt an integration contour $`(\vartheta )`$ passing at distance $`1/n`$ left of the singularity $`z=1`$, then escaping outside of the unit disc within $`\mathrm{\Delta }`$. Upon setting $`z=1+t/n`$, careful approximations yield the result . $`\mathrm{}`$

This theorem allows one to transfer error terms in the asymptotic expansion of a function at its singularity (here $`z=1`$) to asymptotics of the coefficients. The Hankel contour technique is quite versatile and a statement similar to Theorem 2 holds with $`o()`$-conditions. replacing $`O()`$-conditions. The case of $`\alpha `$ being a negative integer is covered by minor adjustments due to $`1/\mathrm{\Gamma }(\alpha )=0`$; see . In concrete terms: *Hankel contours combined with Cauchy coefficient integrals accurately “capture” the singular behaviour of a function.*

By Theorems 1 and 2, whenever a function admits an asymptotic expansion near $`z=1`$ in the basic scale, one has the implication, with $`\sigma \tau \mathrm{}\omega `$,

$$f(z)=\lambda \sigma (z)+\mu \tau (z)+\mathrm{}+O(\omega (z))f_n=\lambda \sigma _n+\mu \tau _n+\mathrm{}+O(\omega _n),$$

where $`f_n=[z^n]f(z)`$. In other words, a *dictionary* translates singular expansions of functions into the asymptotic forms of coefficients. Analytic continuation and validity of functions’ expansions outside of the unit circle is a logical necessity, but once granted, application of the method becomes quite mechanical.

For combinatorics, singularities need not be placed at $`z=1`$. But since $`[z^n]f(z)\rho ^n[z^n]f(z\rho )`$, the dictionary can be used for (dominant) singularities that lie anywhere in the complex plane. The case of finitely many dominant singularities can also be dealt with (via composite Hankel contours) to the effect that the translations of local singular expansions get composed additively. In summary one has from function to coefficients:

location+nature of singularity (fn.) $``$ exponential+polynomial asymptotics (coeff.)

Example 1. *2-Regular graphs* (Comtet, 1974). The class $`𝒢`$ of (labelled) 2-regular graphs can be specified as sets of unordered cycles each of length at least $`3`$. Symbolically:

$$𝒢𝔓\left\{\frac{1}{2}_3\left\{𝒵\right\}\right\}\text{so that}\widehat{G}\left(z\right)=\mathrm{exp}(\frac{1}{2}\mathrm{log}(1z)^1\frac{z}{2}\frac{z^2}{4})=\frac{e^{z/2z^2/4}}{\sqrt{1z}}.$$

$`𝒵`$ represents a single atomic node. The function $`\widehat{G}(z)`$ is singular at $`z=1`$, and

$$𝒢\left(z\right)\underset{z1}{}e^{3/4}\left(1z\right)^{1/2}\frac{G_n}{n!}\underset{n\mathrm{}}{}\frac{e^{3/4}}{\sqrt{\pi n}}.$$

This example can be alternatively treated by Darboux’s method.  $`\mathrm{}`$

Example 2. The *diversity index* of a tree (Flajolet, Sipala & Steyaert, 1990) is the number of non-isomorphic terminal subtrees, a quantity also equal to the size of maximally compact representation of the tree as a directed acyclic graph and related to common subexpression sharing in computer science applications. The mean index of a random binary tree of size $`2n+1`$ is asymptotic to $`Cn/\sqrt{\mathrm{log}n}`$, where $`C=\sqrt{8\mathrm{log}2/\pi }`$. This results from an exact gf obtained by inclusion-exclusion:

$$K\left(z\right)=\frac{1}{2z}\underset{k0}{}\frac{1}{k+1}\left(\genfrac{}{}{0pt}{}{2k}{k}\right)\left(\sqrt{14z+4z^{k+1}}\sqrt{14z}\right).$$

Singularities accumulate geometrically to the right of $`1/4`$ while $`K(z)`$ is $`\mathrm{\Delta }`$–continuable. The unusual singularity type ($`1/\sqrt{X\mathrm{log}X}`$) precludes the use of Darboux’s method.  $`\mathrm{}`$

Rules like those of Table (1.1) preserve analyticity and analytic continuation. Accordingly, generating functions associated with combinatorial objects described by simple construction rules usually have gfs amenable to singularity analysis. The method is systematic enough, so that an implementation within computer algebra systems is even possible as was first demonstrated by Salvy .

## 3. Closure properties

In what follows, we say that a function is *amenable to singularity analysis*, or “*of S.A. type*” for short, if it is $`\mathrm{\Delta }`$-continuable and admits there a singular expansion in the scale $`𝒮=\{\sigma _{\alpha ,\beta }(z)\}`$. First, functions of S.A. type include polylogarithms:

###### Theorem 3

The generalized polylogarithms $`\mathrm{Li}_{\alpha ,k}`$ are of S.A. type, where $`\mathrm{Li}_{\alpha ,k}(z):=_{n1}(\mathrm{log}n)^kn^\alpha z^n`$, $`k_0`$.

The proof makes use of the Lindelöf representation

$$\underset{n1}{}\varphi (n)(z)^n=\frac{1}{2i\pi }_{1/2\mathrm{}}^{1/2+\mathrm{}}\varphi (s)z^s\frac{\pi }{\mathrm{sin}\pi s}𝑑s,$$

in a way already explored by Ford , with Mellin transform techniques providing validity of the singular expansion in a $`\mathrm{\Delta }`$-domain .

Example 3. Entropy computations. The gf of $`\{\mathrm{log}(k!)\}`$ is $`(1z)^1\mathrm{Li}_{0,1}(z)`$, which is of S.A. type. The entropy of the binomial distribution, $`\pi _{n,k}=\left(\genfrac{}{}{0pt}{}{n}{k}\right)p^k(1p)^{nk}`$, results:

$$H_n:=\pi _{n,k}\mathrm{log}\pi _{n,k}\underset{n\mathrm{}}{}\frac{1}{2}\mathrm{log}n+\frac{1}{2}+\mathrm{log}\sqrt{2\pi p\left(1p\right)}+\mathrm{}.$$

Such problems are of interest in information theory, where redundancy estimates precisely depend on higher order asymptotic properties (Jacquet-Szpankowski, 1998). Full expansions for functionals of the Bernoulli distribution are also obtained systematically .  $`\mathrm{}`$

As it is well-known, asymptotic expansions can be integrated while differentiation of asymptotic expansions is permissible for functions analytic in sectors:

###### Theorem 4

Functions of S.A. type are closed under differentiation and integration.

Finally, the Hadamard product of two series, $`fg`$ is defined as the termwise product: $`f(z)g(z)=_nf_ng_nz^n`$, if $`f(z)=_nf_nz^n`$, $`g(z)=_nf_nz^n`$. Hadamard (1898) proved that singularities get composed multiplicatively. Finer composition properties  result from an adaptation of Hankel contours to Hadamard’s formula

$$f(z)g(z)=\frac{1}{2i\pi }_\gamma f(t)g\left(\frac{w}{t}\right)\frac{dt}{t}.$$

###### Theorem 5

Functions of S.A. type are closed under Hadamard product.

Example 4. *Divide-and-conquer* algorithms solve recursively a problem of size $`n`$ by splitting it into two subproblems and recombining the partial solutions. Under the assumption of randomness preservation, the expected costs $`f_n`$ satisfy a “tree recurrence” of the form

$$f_n=t_n+\underset{k}{}\pi _{n,k}\left(f_k+f_{nka}\right),a\{0,1\},$$

where the “toll” sequence $`t_n`$ usually has a simple form (e.g., $`n^\beta ,\mathrm{log}n`$) and the $`\pi _{n,k}`$ characterize the stochastic splitting process. The corresponding gfs then satisfy an equation $`f(z)=t(z)+\left[f\right](z)`$, where the linear operator $``$ reflects the splitting probabilities. For instance, binary search trees and the Quicksort algorithm have $`[f(z)]=2_0^zf(x)𝑑x/(1x)`$. One then has in operator notation $`f(z)=(I)^1[t](z),`$ where the quasi-inverse acts as a “singularity transformer”. Closure theorems allow for an asymptotic classification of the cost functions induced by various tolls under various probabilistic models mirrored by the splitting probabilities .  $`\mathrm{}`$

## 4. Functional equations

*Algebraic functions* have expansions at singularities that are expressed by fractional power series (Newton-Puiseux). Consequently, they are of S.A. type with rational exponents; accordingly their coefficient expansions are linear combinations of algebraic elements of the form $`\omega ^nn^{r/s}`$, with $`\omega `$ and algebraic number and $`r/s`$. By Weierstrass preparation, such properties extend to many *implicit* gfs. For instance, gfs of combinatorial families of trees constrained to have degrees in a finite set have a branch point of type $`\sqrt{X}`$ at their dominant singularity, which in turn corresponds to the “universal” asymptotic form $`T_nc\omega ^nn^{3/2}`$ for coefficients (Pólya, 1937; Otter 1948; Meir-Moon, 1978).

Order constraints in labelled classes are known to correspond to integral operators and, in the recursive case, there result gfs determined by *ordinary differential equations*. Furthermore, moments (i.e., cumulative values) of additive functionals of combinatorial structures defined by recursion and order constraints have gfs that satisfy *linear differential equations*, for which there is a well-established classification theory going back to the nineteenth century. In particular, in the *Fuchsian case*, singularity analysis applies unconditionally, so that the resulting coefficient estimates are linear combinations of terms $`\omega ^nn^\alpha (\mathrm{log}n)^k`$, with $`\omega ,\alpha `$ algebraic and $`k_0`$. This covers a large subset of the class of “holonomic” functions, of which Zeilberger has extensively demonstrated the expressive power in combinatorial analysis .

Example 5. Quadtrees are a way to superimpose a hierarchical partitioning on sequences of points in $`d`$-dimensional space: the first point is taken as the root of the tree and it partitions the whole space in $`2^d`$ orthants in which successive points are placed and then made to refine the partition . The problem is expressed by a linear differential equation with coefficients in $`(z)`$. The average cost of finding a point knowing only $`k`$ of its coordinates is of the asymptotic form $`cn^\alpha `$ with $`\alpha =\alpha (k,d)`$ an algebraic number of degree $`d`$. For instance $`k=1`$ and $`d=2`$ yield a solution involving a $`{}_{2}{}^{}F_{1}^{}`$ hypergeometric function as well as $`\alpha =(\sqrt{17}3)/20.56155`$, in contrast to an exponent $`\frac{1}{2}`$ that would correspond to a perfect partitioning, i.e., a regular grid (Flajolet, Gonnet, Puech & Robson, 1993).  $`\mathrm{}`$

*Substitution equations* correspond to “balanced structures” of combinatorics. An important rôle in the development of the theory has been played by Odlyzko’s analysis  of 2-3 trees (such trees have internal nodes of degree 2 and 3 only and leaves are all at the same depth). The ogf satisfies the equation $`T(z)=z+T(\tau (z))`$, with $`\tau (z):=z^2+z^3,`$ and has a singularity at $`z=1/\varphi `$, a fixed point of $`\tau `$, ($`\varphi =(1+\sqrt{5})/2`$). The singular expansion involves periodic oscillations, corresponding to infinitely many singular exponents having a common real part. Singularity analysis extends to this case and the number of balanced 2-3 trees is found to be of the form $`T_n\frac{1}{n}\varphi ^n\mathrm{\Omega }(\mathrm{log}n),`$ for some nonconstant smooth periodic function $`\mathrm{\Omega }`$.

A similar problem of *singular iteration* arises in the analysis of the height of binary trees . The gf $`y_h(z)`$ of trees of height at most $`h`$ satisfies the Mandelbrot recurrence $`y_h=z+y_{h1}^2`$, with $`y_0=z`$. The fixed point is the gf of binary trees, that is, of Catalan numbers, $`y_{\mathrm{}}=(1\sqrt{14z})/2`$ which has its dominant singularity at $`1/4`$. The analysis of moments of the distribution of height turns out to be equivalent to developing *uniform* approximations to $`y_h(z)`$ as $`z1/4`$ and $`h\mathrm{}`$ simultaneously, this for $`z`$ in a $`\mathrm{\Delta }`$-domain. The end result, by singularity analysis and the moment method, is: the height of a random binary tree with $`n`$ external nodes when normalized by a factor of $`1/(2\sqrt{n})`$ converges in distribution to a theta law defined by the density $`4x_{k1}k^2(2k^2x23)e^{k^2x^2}.`$ The result extends to all simple families of trees in the sense of Meir and Moon and it provides pointwise estimates of the proportion of trees of given height, that is, a *local limit law*.

Generalized digital trees (Flajolet & Richmond, 1992) correspond to a *difference differential equation*, $`_z^k\phi (z)=t(z)+2e^{z/2}\phi (z/2)`$, whose solution involves basic hypergeometric functions. Catalan sums of the form $`_k\left(\genfrac{}{}{0pt}{}{2n}{nk}\right)\nu (k)`$, with $`\nu `$ an arithmetical function, arise in the statistics of “order” (also known as Horton-Strahler number) of trees (Flajolet & Prodinger, 1986). Both cases are first subjected to a Mellin transform analysis, which provides the relevant singular expansions. Periodic fluctuations similar to the case of balanced trees then result from singularity analysis.

A notable parallel to the paradigm of generating functions and singularity analysis has been developed by Vallée in a series of papers. In her framework, singularities of certain transfer operators (of Ruelle type) replace singularities of generating functions. See, *e.g.* for applications to Euclidean algorithms and statistics on sequences produced by a general model of dynamical sources.

## 5. Limit laws

One of the important features of singularity analysis, in contrast with Darboux’s method or (real) Tauberian theory, is to allow for *uniform* estimates. This makes it possible to analyse asymptotically coefficients of multivariate generating functions, $`f(z,u)`$, where the auxiliary variable $`u`$ marks some combinatorial parameter $`\chi `$. One *first* proceeds to extract $`f_n(u):=[z^n]f(z,u)`$ by considering $`f(z,u)`$ as a parameterized family of univariate gfs to which singularity analysis is applied. (The coefficients $`f_n(u)`$ are, up to normalization, probability generating functions of $`\chi `$.) A *second* level of inversion is then achieved by the standard continuity theorems for probability characteristic functions (equivalently Fourier transforms). Technically, consideration of a (small) neighbourhood of $`u=1`$ is normally sufficient for extracting central limit laws.

Two important cases are those of a smoothly varying singularity and of a smoothly varying exponent. In the first case, $`f(z,u)`$ has a constant singular exponent $`\alpha _0`$ and one has $`f(z,u)c(u)(1z/\rho (u))^{\alpha _0}`$. Then, uniformity of singularity analysis implies the estimate $`f_n(u)/f_n(1)(\rho (1)/\rho (u))^n`$. In other words, the probability generating function of $`\chi `$ over objects of size $`n`$ is analytically similar to the gf of a sum of independent random variables—this situation is described as a *“quasi-powers” approximation*. A Gaussian limit law for $`\chi `$ results from the continuity theorem, with mean and variance that grow in proportion to $`n`$. The other case of a smoothly varying exponent is dealt with similarly: one has $`f(z,u)c(u)(1z/\rho )^{\alpha (u)}`$ implying $`f_n(u)/f_n(1)n^{\alpha (u)\alpha (1)}`$; this is once more a quasi-power approximation, but with the parameter now in the scale of $`\mathrm{log}n`$. (See Gao & Richmond, 1992, for hybrid cases.)

The technology above builds on early works of Bender , continued by Flajolet & Soria , and H. K. Hwang . In particular, under general conditions, the following hold: a local limit law expresses convergence to the Gaussian density; speed of convergence estimates result from the Berry-Esseen inequalities; large deviation estimates derive from singularity analysis applied at fixed real values $`u1`$.

Example 6. *Polynomials over finite fields*. Consider the family $`𝒫`$ of all polynomials with coefficients in the Galois field $`𝔽_q`$. A polynomial being determined by its sequence of coefficients, the gf $`P(z)`$ of all polynomials has a polar singularity. Furthermore, the unique factorization property implies that $`𝒫`$ is isomorphic to the class of all multisets ($`𝔐`$) of the irreducible polynomials $``$: $`𝒫𝔐\{\}`$. Since taking multisets corresponds to exponentiating singularities of gfs, the singularity of the gf $`I(z)`$ is logarithmic. By singularity analysis, the number of irreducible polynomials is asymptotic to $`q^n/n`$—this is an analogue of the prime number theorem, which was already known to Gauß. The bivariate gf of the number of irreducible factors in polynomials turns out to be of the singular type $`(1qz)^u`$, with a smooth variable exponent, so that: *the number of irreducible factors of a random polynomial over $`𝔽_q`$ is asymptotically Gaussian*. This parallels the Erdős-Kac theorem for integers. Similar developments lead to a complete analysis of a major polynomial factorization algorithm (Flajolet, Gourdon & Panario, 2001).  $`\mathrm{}`$

Movable singularities and exponents occur frequently in the analysis of parameters defined by recursion, leading to algebraic or differential equations, which “normally” admit a smooth perturbative analysis.

Example 7. *Patterns in random strings.* Let $`\mathrm{\Omega }`$ be the total number of occurrences of a fixed pattern (as a contiguous block) in a random string over a finite alphabet. For either the Bernoulli model, where letters are independently identically distributed, or the Markov model, the bivariate gf, with $`z`$ marking the length of the random string and $`u`$ the number $`\mathrm{\Omega }`$ of occurrences, is a *rational function*, as it corresponds to a finite-state device. Perron-Frobenius properties apply for positive $`u`$. Therefore the bivariate gf viewed as a function of $`z`$ has a simple dominant pole at some $`\rho (u)`$ that is an algebraic (and holomorphic) function of $`u`$, for $`u>0`$. Quasi-powers approximations therefore hold and the limit law of $`\mathrm{\Omega }`$ in random strings of length $`n`$ is Gaussian. Such facts holds for very general notions of patterns and are developed systematically in Szpankowski’s book .  $`\mathrm{}`$

Example 8. *Non-crossing graphs.* Consider graphs with vertex set the $`n`$th roots of unity, constrained to have only non-crossing edges; let the parameter $`\chi `$ be the number of connected components. The bivariate gf $`G(z,u)`$ is an *algebraic function* satisfying

$$G^3+\left(2w^3z^23w^2z+w3\right)G^2+\left(3w^2z2w+3\right)G+w1=0.$$

$`G(z,1)`$ has a dominant singularity at $`\rho (1)=3/2\sqrt{2}`$ which gets smoothly perturbed to $`\rho (u)`$ for $`u`$ near 1. The singularity type is consistently of the form $`(1z/\rho (u))^{1/2}`$. A central limit law results for the number of components in such graphs (Flajolet–Noy, 1999). Drmota has given general conditions ensuring Gaussian laws for problems similarly modelled by multivariate algebraic functions .  $`\mathrm{}`$

Example 9. *Profile of quadtrees*. Refer to Example 4.. The bivariate gf $`f(z,u)`$ of node levels in quadtrees satisfies an equation, which, for dimension $`d=3`$ reads

$$f(z,u)=1+2^3u_0^z\frac{dx_1}{x_1\left(1x_1\right)}_0^{x_1}\frac{dx_2}{x_2\left(1x_2\right)}_0^{x_2}f(x_3,u)\frac{dx_3}{1x_3}.$$

This corresponds to a *linear differential equation* with coefficients in $`(z,u)`$ and a fixed singularity at $`z=1`$. The indicial equation is an algebraic one parameterized by $`u`$ and, when $`u1`$, there is a unique largest branch $`\alpha (u)`$ that determines the dominant regime of the form $`(1z)^{\alpha (u)}`$. This is the case of a movable exponent inducing a central limit law: *The level profile of a $`d`$-dimensional quadtree is asymptotically Gaussian.* Such properties are expected in general for models that are perturbations of linear differential equations with a fixed Fuchsian singularity (Flajolet & Lafforgue, 1994).  $`\mathrm{}`$

Finally, singularity analysis also intervenes by making it possible to “pump” moments of combinatorial distributions. Examples include the height of trees discussed earlier, as well as tree path length (Louchard 1983, Takács 1991) and the construction cost of hashing tables (Flajolet, Poblete & Viola, 1998). The latter problems were first shown in this way to converge to Brownian Excursion Area.

## 6. Conclusions

Elementary combinatorial structures are enumerated by generating functions that satisfy a rich variety of functional relations. However, the singular types that are observed are usually somewhat restricted, and *driven by combinatorics*. In simpler cases, the generating functions are explicit combinations of a standard set of special functions. Next, implicitly defined functions (associated with recursion) have singularities that arise from failures of the implicit function theorem and are consequently of the algebraic type, often with exponent $`\frac{1}{2}`$. Linear differential equations have a well-established classification theory that, in the Fuchsian case, leads to algebraic-logarithmic singularities. In all such cases, the singular expansion is known to be valid outside of the original disc of convergence of the generating function. This means that singularity analysis is *automatically* applicable, and precise asymptotic expansions of coefficients result.

Parameters of combinatorial structures, provided they remain “simple” enough, lead to local deformations (via an auxiliary variable $`u`$ considered near 1) of the functional relations defining univariate counting generating functions. Under fairly general conditions, such deformations are amenable to perturbation theory and admit of uniform expansions near singularities. Then, since the singularity analysis process preserves uniformity, limit laws result via the continuity theorem for characteristic functions. In this way, the behaviour of a large number of parameters of elementary combinatorial structures becomes predictable. (The theory of functions of several complex variables is thus bypassed. See Pemantle’s recent work  based on this theory for a global characterization of all the asymptotic regimes involved.)

The generality of the singular approach makes it even possible to discuss *combinatorial schemas* at a fair level of generality , the case of polynomial factorization (Ex. 5.) being typical. Roughly, combinatorial constructions viewed as “singularity transformers” dictate asymptotic regimes and probabilistic laws. Analytic combinatorics then represents an attractive alternative to probabilistic methods, whenever a strong analytic structure is present—this is the case for most combinatorial problems that are “decomposable” and amenable to the generating function methodology. Very precise asymptotic information on the randomness properties of large random combinatorial objects results from there. This in turn has useful implications in the analysis of many fundamental algorithms and data structures of computer science, following the steps of Knuth’s pioneering works .


format: text

text:
# Wetting under non-equilibrium conditions

## I Introduction

Wetting phenomena occur in a large variety of experiments, where a planar substrate is exposed to a gas phase under thermal equilibrium conditions. Generally, ‘wetting’ refers to a situation where a bulk phase $`(\alpha )`$ in contact with a substrate coexists with a layer of a different phase ($`\beta `$) which is preferentially attracted to the surface of the substrate. By changing physical parameters such as temperature and chemical potential, the system may undergo a wetting transition from a non-wet phase, where the thickness of the layer stays finite, to a wet phase, where the layer becomes macroscopic.

The phase diagram associated with the surface layer could be rather complex exhibiting a variety of surface phase transitions, prewetting phenomena and multicritical behavior WettingReview ; NakanishiFisher . For example, by increasing the temperature $`T`$ while moving along the $`\alpha `$-$`\beta `$ coexistence curve, a wetting transition may take place at a temperature $`T_W`$, beyond which the thickness of the layer becomes infinite. Usually this transition is first order, although in certain models the transition is continuous, and is then referred to as continuous wetting. On the other hand, when the chemical potential difference between the two phases is varied, moving towards the coexistence curve at $`T>T_W`$, a different type of transition takes place in which the thickness of the layer diverges. This phenomenon is known as complete wetting.

In many experimental situations it is reasonable to assume that a wetting stationary layer is in thermal equilibrium. In fact, methods of equilibrium statistical mechanics turned out to be very successful in a large variety of theoretical and experimental studies (for a review, see Ref. WettingReview ). Within this approach, a wetting transition is usually considered as the unbinding of an interface from a wall. The interface configuration is described by a function $`h(\text{x})`$ which gives the height of the interface at point x on the substrate. One then introduces an effective Hamiltonian of the form EquilibriumFieldTheory

$$=d^{d1}x\left[\frac{\sigma }{2}(h)^2+V\left(h(\text{x})\right)\right],$$
(1)
where $`\sigma `$ is the effective surface tension of the $`\alpha `$-$`\beta `$ interface, $`V(h)`$ is a potential accounting for the interaction between the wall and the interface, and $`d1`$ is the interface dimension. In the non-wet phase the potential $`V`$ contains an attractive component which binds the interface to the wall. Assuming thermal equilibrium, the probability of finding the interface in a certain configuration is then given by the canonical distribution

$$P[h]\mathrm{exp}\left(\beta [h]\right).$$
(2)
As the parameters describing the system are varied, the attractive component of the potential may become weaker so that it is no longer able to bind the interface, leading to a wetting transition.

In order to study wetting phenomena under thermal equilibrium conditions one usually introduces a stochastic Langevin equation corresponding to the effective Hamiltonian (1). This Langevin dynamics should reproduce the equilibrium distribution (2) in the limit $`t\mathrm{}`$. Since many different dynamical rules may approach the same stationary state, this condition does not fully determine the form of the Langevin equation. However, assuming short-range interactions and keeping only the most relevant terms in the renormalization group sense, one is led to the Edwards-Wilkinson equation with a potential BarabasiStanley

$$\frac{h(\text{x},t)}{t}=\sigma ^2h(\text{x},t)\frac{V\left(h(\text{x},t)\right)}{h(\text{x},t)}+\zeta (\text{x},t),$$
(3)
where $`\zeta (\text{x},t)`$ is a zero-average Gaussian noise field with a variance

$$\zeta (\text{x},t)\zeta (\text{x}^{},t^{})=2\mathrm{\Gamma }\delta ^{d1}(\text{x}\text{x}^{})\delta (tt^{}),$$
(4)
and a noise amplitude $`\mathrm{\Gamma }=k_BT`$. This Langevin equation has the same symmetry properties as the Hamiltonian (1), namely, it is invariant under translations, rotations, and reflections in space. Apart from the potential term, the equation is also invariant under shifts $`hh+a`$ and reflections $`hh`$. Moreover, it can be shown that this type of Langevin dynamics obeys detailed balance and relaxes towards the equilibrium distribution (2) in the bound phase.

Wetting phenomena may also take place in many systems under non-equilibrium conditions. For example, in growth processes such as molecular beam epitaxy or others, a layer is grown on a substrate, whose properties depend on the growth conditions. By varying these conditions one expects wetting phenomena to take place. Here, unlike the equilibrium case, the dynamics does not obey detailed balance, leading to a rather different class of wetting phenomena.

The simplest way to study nonequilibrium wetting on the level of the Langevin equation is to introduce a nonlinear term in Eq. (3), leading to a Kardar-Parisi-Zhang (KPZ) equation with a potential KPZWall

$$\frac{h(\text{x},t)}{t}=\sigma ^2h(\text{x},t)\frac{V\left(h(\text{x},t)\right)}{h(\text{x},t)}+\lambda \left(h(\text{x},t)\right)^2+\zeta (\text{x},t).$$
(5)
It is important to note that this nonlinear term is a relevant perturbation of the underlying field theory, i.e., even if $`\lambda `$ is very small, it will be amplified under renormalization group transformations, driving the system away from thermal equilibrium.

Recently, a simple solid-on-solid (SOS) model for non-equilibrium wetting in $`1+1`$ dimensions was introduced Alon ; Wetting1 . The model is controlled by an adsorption rate $`q`$ and a desorption rate $`p`$ and exhibits a continuous wetting transition at a critical growth rate $`q_c(p)`$. The wetting transition is related to the unpinning process of an interface from a substrate and may be described by the KPZ equation (5). The model has then been generalized to include a short-range interaction between the interface and the substrate Wetting2 . This was done by introducing a modified growth rate $`q_0`$ at the substrate level. This results in a contact interaction between the interface and the substrate, which is attractive for $`q_0<q`$ and repulsive for $`q_0>q`$. It was found that sufficiently strong attractive interaction modifies the nature of the wetting transition, making it first order. In addition, it has been demonstrated that there exists an extended region in the phase diagram, where the pinned and the moving phases coexist in the sense that the transition time from the pinned to the moving phase grows exponentially with the system size so that the two phases become stable in the thermodynamic limit. It should be emphasized that this kind of phase coexistence, which has also been observed in the past in other models Toom80 ; BennettGrinstein85 ; Bridge95 ; Mueller97 , can only occur in nonequilibrium systems.

Some of these results have been confirmed by numerical and mean field studies of KPZ type models GiadaMarsili ; Santos . In particular, by adding a repulsive interaction between the interface and the substrate, a crossover from a continuous to a first order wetting transition was found. Non-equilibrium wetting phenomena have also been studied recently in models of growing magnetic domains CandiaAlbano .

In this paper we present a detailed account of the phase diagram of the SOS model in 1+1 dimensions. In order to get some indication on the behavior of the model in higher dimensions we introduce a mean field approximation for the model, and study the resulting phase diagram and the dynamical behavior of the interface. It is found that the mean field approximation reproduces the main qualitative features of the phase diagram obtained in 1+1 dimensions. In addition, a new feature is found in the pinned phase. By studying the height distribution of the pinned interface, two distinct types of behavior are discovered upon varying the dynamical parameters of the model: in one regime the distribution decays superexponentially with the height, while in the other the decay is exponential. The two regimes are separated by a prewetting transition line. Relevance of this phase diagram to wetting phenomena in dimensions higher than $`1+1`$ is considered.

The paper is organized as follows. In Sec. II we recall the definition of the SOS model and summarize its properties. The special case $`p=1`$, where the model is exactly soluble, is discussed in detail in Sec. III. The mean field approximation, which is the main focus of the present work, is presented and analyzed in Sec. IV. The paper ends with concluding remarks in Sec. V.

## II Definition and properties of the model

The SOS model defined in Wetting1  is probably the simplest model which follows the spirit of Eq. (5). It is defined on a one-dimensional lattice with $`N`$ sites and periodic boundary conditions, where each site $`i`$ is associated with an integer variable $`h_i`$ describing the local height of the interface. The repulsive part of the potential $`V(h)`$ is implemented as a hard-core wall at zero height, restricting the heights $`h_i`$ to be non-negative. Thus $`h_i=0,1,2,\mathrm{},`$ . Moreover, an effective surface tension is introduced by imposing the restricted solid-on-solid (RSOS) condition

$$|h_ih_{i\pm 1}|1.$$
(6)
The model evolves random-sequentially by choosing a random site and carrying out one of the following processes (see Fig. 1):

* Deposition of an atom on the substrate with rate $`q_0`$:

$$h_i=0h_i=1$$
(7)
* Deposition of an atom on top of already deposited islands with rate $`q`$:

$$h_ih_i+1\text{if }h_i1$$
(8)
* Evaporation of an atom at the edge of a terrace with rate $`r`$:

$$h_i\mathrm{min}(h_{i1},h_i,h_{i+1})$$
(9)
* Evaporation of an atom in the middle of a plateau with rate $`p`$:

$$h_ih_i1\text{if}h_{i1}=h_i=h_{i+1}>0$$
(10)

If the selected process would violate the restrictions $`h_i0`$ or $`|h_ih_{i\pm 1}|1`$, the attempted move is abandoned and a new site $`i`$ is selected. Each attempted update corresponds to a time increment $`\mathrm{\Delta }t=1/N`$. Since one of the four rates can be chosen freely by rescaling time, we set $`r=1`$. Thus the model is controlled by three parameters, namely, a growth rate $`q`$, a desorption rate $`p`$, and a special growth rate $`q_0`$ at the substrate which accounts for an additional short-range interaction between the substrate and the wetting layer.

The model can be easily generalized to higher dimensions. Note that in this case it is possible to introduce different evaporation rates for various types of edge sites, e.g. linear edges and corner sites. For simplicity we will assume that all rates for evaporation at edges are equal to $`1`$. In this case the model can simply be generalized to higher dimensions by including all nearest neighbors in Eqs. (9) and (10).

### II.1 Properties for $`q_0=q`$

Let us first consider the case without interactions between substrate and wetting layer, i.e., $`q_0=q`$. In this case the presence of a hard-core wall at zero height leads to a continuous phase transition between a bound and a moving phase (see Fig. 2). The phase transition line $`q_c(p)`$ is determined by a vanishing propagation velocity of a freely evolving interface. For $`q<q_c`$ the interface moves downward until it fluctuates close to the wall while for $`q>q_c`$ the propagation velocity is positive and the interface detaches from the wall. Obviously this transition takes place even in finite systems with a critical threshold depending on $`N`$. We note that for $`p=0`$, where evaporation from the middle of plateaus is forbidden, the interface velocity in the bulk cannot be negative so that in this case the transition relies on a different mechanism Alon .

One can easily verify that the dynamical rules (7)-(10) do not generally satisfy detailed balance. To this end we consider the closed cycle of deposition and evaporation processes, as shown in Fig. 3. Let the statistical weights of these configurations in the steady state be $`P_A,P_B,P_C`$, and $`P_D`$. Obviously detailed balance implies

$$q^2P_A=qP_B=P_C=\frac{q}{p}P_D=\frac{q^2}{p}P_A.$$
(11)
These equations can only be satisfied if $`p=1`$ (as already anticipated, $`r`$ is assumed to be equal to 1 without any loss of generality). This special case will be discussed in more detail in Sec. III.

The two parameters $`p`$ and $`q`$ can be used to control $`V/h`$ and $`\lambda `$ in the KPZ equation (5). In the case of detailed balance, where a bound interface is thermally equilibrated, the coefficient $`\lambda `$ is expected to vanish. This can be verified by comparing the propagation velocities of a horizontal and an artificially tilted interface far away from the wall. The line, where both velocities coincide, is shown as a dashed line in Fig. 2. As expected, it crosses the phase transition line exactly at the point $`p=q=1`$, where detailed balance is satisfied.

Moving away from this point, we can therefore study the crossover from equilibrium to non-equilibrium wetting. Numerical simulations suggested that the transitions for $`p<1`$ and $`p>1`$ are associated with different sets of critical exponents Wetting1 . These findings are in accordance with results obtained by Munõz and Hwa KPZWall , who showed that the scaling properties of a KPZ interface interacting with a wall depend on the sign of $`\lambda `$.

### II.2 Properties for $`q_0<q`$

In most experimental applications the wetting layer interacts with the substrate, giving rise to an additional short-range force at the bottom layer which may be either attractive or repulsive. In the present model such a force can be taken into account by introducing a different growth rate $`q_0`$ for deposition at zero height. The influence of this parameter was studied in detail in Ref. Wetting2 . Since $`q_0`$ does not influence the propagation velocity of the interface far away from the wall, the location of the transition line, along which the moving phase is stable, remains unchanged. However, if $`q_0`$ is smaller than a certain threshold $`q_0^{}`$, the attractive interaction is so strong as to stabilize the bound phase even when the free interface would grow. In other words, for very low values of $`q_0`$, there exists an extended region in the phase diagram where the bound and the moving phase coexist in the sense that the transition time from the bound to the moving phase grows exponentially with the systems size. The upper boundary of the coexistence region depends on $`q_0`$, as shown in Fig. 4.

A thermodynamically stable coexistence of the bound and the moving phase requires a robust mechanism which eliminates large protruding islands in the bound phase. In the present model this mechanism works as follows. Once an island has been formed by fluctuations, the detached part of the interface quickly grows since $`q>q_c`$. In the phase coexistence region, where the coefficient $`\lambda `$ in the KPZ equation is negative, the island will grow until the slope at the edges exceeds a critical value, where the growth is compensated by the nonlinear term. Afterwards, the pyramidial island shrinks linearly with time until it is eliminated. Therefore, phase coexistence can only occur under non-equilibrium conditions in those regions of the phase diagram where $`\lambda `$ is negative. Very recently, the tricritical point and the critical behavior at the upper boundary of the phase coexistence region has been investigated in a discretized KPZ equation with a potential Ginelli .

## III Exactly soluble case: $`p=1`$

### III.1 Detailed balance and transfer matrix approach

We first consider the special case $`p=1`$, where detailed balance is satisfied. In this case the stationary probability distribution of the bound interface configurations $`\{h_1,\mathrm{},h_N\}`$ is given by the canonical ensemble expression corresponding to an energy functional $``$

$$P(h_1,\mathrm{},h_N)=\frac{1}{Z_N}\mathrm{exp}\left[(h_1,\mathrm{},h_N)\right].$$
(12)
The partition sum

$$Z_N=\underset{h_1,\mathrm{},h_N}{}\mathrm{exp}\left[(h_1,\mathrm{},h_N)\right]$$
(13)
runs over all interface configurations obeying the RSOS constraint (6). The energy $``$ is given by

$$(h_1,\mathrm{},h_N)=\underset{i=1}{\overset{N}{}}V(h_i),$$
(14)
where $`V(h)`$ is a potential of the form

$$V(h)=\{\begin{array}{cc}\mathrm{}\hfill & \text{if }h<0\hfill \\ \mathrm{ln}(q/q_0)\hfill & \text{if }h=0\hfill \\ h\mathrm{ln}(q)\hfill & \text{if }h>0\hfill \end{array}$$
(15)
as sketched in Fig. 5. Therefore, Eq. (12) may be rewritten as

$$P(h_1,\mathrm{},h_N)=Z_N^1q^{(_{i=1}^Nh_i)}(q/q_0)^{(_{i=1}^N\delta _{h_i,0})}.$$
(16)
Note that in this expression $`\mathrm{ln}(q)`$ plays the role of the chemical potential difference $`\mathrm{\Delta }\mu /k_BT`$ between the wetting layer and the gas phase. Obviously the transition takes place at $`q_c=1`$.

To verify the validity of this probability distribution, it is sufficient to demonstrate that the dynamical rules obey detailed balance with respect to it. In fact, the deposition process (8) reduces the probability $`P`$ by a factor of $`q`$ while the evaporation processes (9) and (10), which both take place with the same rate, increase $`P`$ by a factor of $`1/q`$. Consequently, the probability currents between pairs of configurations compensate each other so that detailed balance is satisfied, proving the validity of the equilibrium ensemble (16). Similarly one can show that detailed balance is also satisfied at the bottom layer. We note that these considerations are only valid in the bound phase $`q<1`$, where the probability distribution is stationary. Once the system enters the moving phase, the process is out of equilibrium.

The canonical ensemble can be used to compute the density profile of a bound interface in the case of detailed balance. To this end we use a transfer matrix formalism introduced in Hilhorst ; Burkhardt , writing the Boltzmann factor $`\mathrm{exp}()`$ in Eqs. (12) and (16) as a product

$$P(h_1,\mathrm{},h_N)=\frac{1}{Z_N}\underset{i=1}{\overset{N}{}}T_{h_i,h_{i+1}},$$
(17)
where $`Z_N=\text{Tr}(T^N)`$ and

$$T_{h,h^{}}=\{\begin{array}{cc}q^{(h+h^{})/2}(q/q_0)^{(\delta _{h,0}+\delta _{h^{},0})/2}\hfill & \text{ if }|hh^{}|1\hfill \\ 0\hfill & \text{ otherwise. }\hfill \end{array}$$
(18)
The transfer matrix $`T`$ is infinite-dimensional, acts in spatial direction, and yields the contribution to the Boltzmann factor between adjacent sites with the heights $`h`$ and $`h^{}`$. Because of the RSOS condition (6), it has a tridiagonal structure and reads

$$T=\left(\begin{array}{ccccccc}q/q_0& q/q_0^{1/2}& & & & & \\ q/q_0^{1/2}& q& q^{3/2}& & & & \\ & q^{3/2}& q^2& q^{5/2}& & & \\ & & q^{5/2}& q^3& q^{7/2}& & \\ & & & \mathrm{}& \mathrm{}& \mathrm{}& \\ & & & & \mathrm{}& \mathrm{}& \mathrm{}\end{array}\right).$$
(19)
Using the transfer matrix formalism, the stationary density $`\rho (h)`$ of sites at height $`h`$ can be expressed as

$$\rho (h)=Z_N^1h|T^N|h,$$
(20)
where $`\{|h\}`$, $`\{h|\}`$ denote canonical basis vectors in height space. For $`N\mathrm{}`$ this expression is governed by the largest eigenvalue $`\mathrm{\Lambda }`$ of the transfer matrix

$$T|\varphi =\mathrm{\Lambda }|\varphi ,$$
(21)
where $`|\varphi `$ denotes the corresponding eigenvector. Thus, in an infinite system, the stationary densities may also be written as

$$\rho (h)=\frac{|h|\varphi |^2}{\varphi |\varphi }.$$
(22)
Note that the numerator in this expression is quadratic in $`|\varphi `$, just as in a quantum-mechanical problem.

### III.2 The case $`q_0=q`$

In order to understand how the interface detaches from the wall, it is useful to study the scaling behavior of the density of sites at the bottom layer close to the transition point. Let us first consider the case $`q_0=q`$, where interactions between substrate and bottom layer are absent. For $`h>0`$ the eigenvalue problem reads

$$q^{h1/2}\varphi _{h1}+q^h\varphi _h+q^{h+1/2}\varphi _{h+1}=\mathrm{\Lambda }\varphi _h,$$
(23)
where $`\mathrm{\Lambda }`$ is the largest eigenvalue of $`T`$ and $`\varphi _h`$ are the components of the corresponding eigenvector $`|\varphi `$ representing the stationary state. Close to criticality we can carry out the continuum limit $`\varphi _h\varphi (\stackrel{~}{h})`$, replacing discrete heights $`h`$ by real-valued heights $`\stackrel{~}{h}`$. In this limit the above eigenvalue problem turns into a differential equation which, to leading order in $`ϵ=1q>0`$, is given by

$$\left(\frac{^2}{\stackrel{~}{h}^2}+(3\mathrm{\Lambda })3ϵ\stackrel{~}{h}\right)\varphi (\stackrel{~}{h})=0.$$
(24)
This differential equation is solved by an Airy function

$$\varphi (\stackrel{~}{h})=\text{Ai}\left(\frac{3ϵ\stackrel{~}{h}+\mathrm{\Lambda }3}{(3ϵ)^{2/3}}\right).$$
(25)
Since $`\varphi (\stackrel{~}{h})`$ has to vanish for $`\stackrel{~}{h}0`$, we obtain the eigenvalue $`\mathrm{\Lambda }=3`$ so that $`\varphi (\stackrel{~}{h})=\text{Ai}([3ϵ]^{1/3}\stackrel{~}{h})`$. Therefore, the heights, in particular the average height and the interface width, scale as

$$\stackrel{~}{h}wϵ^{1/3}.$$
(26)
Next, we determine the density of exposed sites at the substrate $`\rho (0)`$. In the continuum limit of Eq. (22) $`0|\varphi `$ is proportional to $`\varphi ^{}(0)`$, hence

$$\rho (0)=\frac{1}{𝒩}[\varphi ^{}(0)]^2,$$
(27)
where $`𝒩=_0^{\mathrm{}}𝑑\stackrel{~}{h}\varphi ^2(\stackrel{~}{h})`$ is a normalization factor. Since $`\varphi ^{}(0)ϵ^{1/3}`$ and $`𝒩ϵ^{1/3}`$ one obtains a linear scaling law

$$\rho (0)ϵ.$$
(28)
Thus the density of exposed sites at the bottom layer scales linearly with the distance from criticality, proving that the wetting transition at $`p=q_c=1`$ is continuous.

Imposing fixed boundary conditions $`h_1=h_N=0`$, it is also possible to study finite-size scaling at the critical point $`q_0=q=1`$. At criticality the transfer matrix has a simple structure and can be thought of as generating a simple random walk near a wall so that the mean height and the bottom layer density scale as

$$\stackrel{~}{h}N^{1/2},\rho (0)N^{3/2}.$$
(29)
The transfer matrix does not provide any information regarding dynamical properties. Numerical simulations (whose details are not shown here) suggest that an initially flat interface at zero height roughens with time (for $`tN^2`$) as

$$\stackrel{~}{h}t^{1/4},\rho (0)t^{3/4}.$$
(30)
Assuming standard power law scaling, we can combine Eqs. (26), (28), (29), and (30) in the scaling forms

$$\begin{array}{cc}\hfill \stackrel{~}{h}(ϵ,N,t)& =N^{1/2}f(t/N^2,ϵN^{3/2})\hfill \\ \hfill \rho _0(ϵ,N,t)& =N^{3/2}g(t/N^2,ϵN^{3/2})\hfill \end{array}$$
(31)
where $`f`$ and $`g`$ are scaling functions with an appropriate asymptotic behavior. As expected these scaling forms are consistent with the critical exponents $`z=2,\alpha =1/2`$ of the Edwards-Wilkinson universality class BarabasiStanley .

Another interesting aspect is the stationary distribution $`P(\mathrm{})`$ of island sizes $`\mathrm{}`$ in the bound phase (see Fig. 6). In the case of detailed balance this distribution can be computed exactly. To this end we introduce the projection operator $`P=1|00|`$ which projects onto states with nonzero height. Moreover, let $`Q=PTP`$ be a transfer matrix describing an interface that does not touch the bottom layer. Obviously the distribution $`P(\mathrm{})`$, which may be interpreted as a first-return probability of the interface to the bottom layer, is given by

$$P(\mathrm{})=\frac{0|TQ^\mathrm{}2T|0}{0|T^{\mathrm{}}|0}\frac{0|TQ^\mathrm{}2T|0}{\mathrm{\Lambda }^{\mathrm{}}},$$
(32)
Note that the leftmost column and the topmost row of the restricted transfer matrix $`Q`$ are zero, while all other matrix elements are the same as in Eq. (19). Because of this simple structure, the spectrum of $`Q`$ is just the spectrum of $`T`$ multiplied by $`q`$ combined with a zero mode so that the largest eigenvalue of $`Q`$, which dominates the matrix product in Eq. (32), is $`q\mathrm{\Lambda }`$. In the limit $`\mathrm{}\mathrm{}`$ we therefore obtain an exponential distribution of the form

$$P(\mathrm{})q^{\mathrm{}}.$$
(33)
Therefore, the average island size scales as

$$\overline{\mathrm{}}\frac{1}{\mathrm{ln}q}\frac{1}{ϵ}$$
(34)
and diverges at the transition. Since $`\rho (0)=1/\overline{\mathrm{}}`$, this result is in agreement with Eq. (28).

In order to determine the stationary interface profile $`\rho (h)`$ in the limit $`h\mathrm{}`$, let us go back to Eq. (23). By assuming that the second and the third term on the left hand side can be neglected, one finds that

$$\varphi _h\frac{1}{\mathrm{\Lambda }}q^{h1/2}\varphi _{h1},$$
(35)
with the corresponding asymptotic solution

$$\varphi _h\mathrm{\Lambda }^hq^{h^2/2}.$$
(36)
Therefore, in the case of detailed balance, the profile of a bound interface decays as a Gaussian for large values of $`h`$. Even for $`p1`$, where detailed balance is violated, numerical simulations (see Fig. 7) suggest that in 1+1 dimensions the entire bound phase is characterized by Gaussian density profiles.

### III.3 First-order phase transition for small values of $`q_0`$

Let us now consider the influence of an attractive force between substrate and wetting layer by taking $`q_0<q`$. Obviously, the interface is bound to the wall for $`q<1`$, while for $`q>1`$ it will tunnel through the potential barrier. This means that the critical point $`q_c=1`$ remains unchanged. However, if $`q_0`$ is decreased below a certain threshold $`q_0^{}`$, the attraction is sufficiently strong such that the transition becomes first order.

To demonstrate the crossover to a discontinuous transition in the case of detailed balance, we look for a localized pinned interface solution at the transition point $`q=1`$. Assuming an exponential interface profile

$$\varphi _h=z^h\text{ for }h1$$
(37)
with some $`z<1`$, the eigenvalue problem then reduces to three independent equations

$$\begin{array}{cc}\hfill q_0^1\varphi _0+q_0^{1/2}z& =\mathrm{\Lambda }\varphi _0\hfill \\ \hfill q_0^{1/2}\varphi _0+z+z^2& =\mathrm{\Lambda }z\hfill \\ \hfill z^1+1+z& =\mathrm{\Lambda }\hfill \end{array}$$
(38)
which have the (unnormalized non-negative) solution

$$\varphi _0=q_0^{1/2},z=\frac{\sqrt{1+2q_03q_0^2}}{2(1q_0)}\frac{1}{2},\mathrm{\Lambda }=\frac{z+1}{q_0}.$$
(39)
Consequently, the stationary density of exposed sites at the bottom layer is

$$\rho (0)=\frac{\varphi _0^2}{_{h=0}^{\mathrm{}}\varphi _h^2}=\frac{q_0}{q_0+\frac{z^2}{1z^2}}=\frac{1+q_06q_0^2+\sqrt{1+2q_03q_0^2}}{2+4q_06q_0^2},$$
(40)
while the densities at higher levels are given by

$$\rho (h)=\frac{\varphi _h^2}{_{h=0}^{\mathrm{}}\varphi _h^2}=\frac{z^{2h}}{q_0+\frac{z^2}{1z^2}}.$$
(41)
It turns out that the bottom layer density $`\rho (0)`$ is positive for $`q_0<2/3`$ and vanishes at $`q_0=2/3`$. For $`q_0>2/3`$, however, one has $`z>1`$, so that the exponential ansatz $`\varphi _h=z^h`$ is no longer physically meaningful. Hence, in the case of detailed balance, the transition becomes first order at the tricritical point

$$p=q_c=1,q_0^{}=\frac{2}{3}.$$
(42)
For $`q_0<\frac{2}{3}`$, the potential well is deep enough to bind the critical interface to the wall, leading to an exponentially decaying interface profile. For $`q_0>\frac{2}{3}`$, such a localized solution does not exist and the transition becomes continuous.

### III.4 Scaling properties near the tricritical point

The phase diagram for $`p=1`$ is shown in Fig. 8. Using the transfer matrix approach, we can show that the critical behavior along the second order transition line is always the same as the one discussed in Sec. III.2. However, in the vicinity of the tricritical point the scaling properties are different. Moreover, they depend on the direction from which the tricritical point is approached.

Let us first consider the case $`q=1`$, approaching the tricritical point horizontally from the left along the first-order phase transition line. Using the expressions (40) and (41) one can compute the interface height

$$\overline{h}=\underset{h=0}{\overset{\mathrm{}}{}}h\rho (h)=\frac{2q_0}{1+q_06q_0^2+\sqrt{1+2q_03q_0^2}}$$
(43)
and the squared interface width

$$w^2=\underset{h=0}{\overset{\mathrm{}}{}}(h\overline{h})^2\rho (h)=\frac{2\left(1+4q_03q_0^3+(q_0^23q_01)\sqrt{1+2q_03q_0^2}\right)}{q_0^2\left(1+3q_03\sqrt{1+2q_03q_0^2}\right)^2}.$$
(44)
Approaching the tricritical point from the left by increasing $`q_0`$, these quantities scale to lowest order in $`\delta =q_0^{}q_0`$ as

$$\begin{array}{cc}& \overline{h}=w=\frac{1}{6\delta }\hfill \\ & \rho (0)=4\delta \hfill \end{array}$$
(45)
On the other hand, if the tricritical point is approached vertically keeping $`q_0=q_0^{}=2/3`$ fixed, a numerical diagonalization of the transfer matrix suggests that the asymptotic interface profile crosses over from an exponential to a Gaussian decay. In this case the height, the width, and the bottom layer density are found to scale as

$$\begin{array}{cc}& \overline{h}wϵ^{1/3}\hfill \\ & \rho (0)ϵ^{1/3},\hfill \end{array}$$
(46)
where $`ϵ=1q`$. Moreover, by keeping the boundary sites fixed at zero height, we can study finite-size scaling at the tricritical point. Evaluating products of the transfer matrix numerically we find that

$$\begin{array}{cc}& \overline{h}N^{1/2}\hfill \\ & \rho (0)N^{1/2}.\hfill \end{array}$$
(47)
Finally, numerical Monte Carlo simulations at the tricritical point (see Fig. 9) suggest that an initially flat interface roughens in such a way that

$$\begin{array}{cc}& \overline{h}t^{1/4}\hfill \\ & \rho (0)t^{1/4},\hfill \end{array}$$
(48)
Assuming standard power law scaling all these results can be combined in the scaling forms

$$\begin{array}{cc}\hfill \stackrel{~}{h}(ϵ,\delta ,N,t)& =N^{1/2}F(t/N^2,ϵN^{3/2},\delta N^{1/2})\hfill \\ \hfill \rho _0(ϵ,\delta ,N,t)& =N^{1/2}G(t/N^2,ϵN^{3/2},\delta N^{1/2})\hfill \end{array}$$
(49)
where $`F`$ and $`G`$ are scaling functions with an appropriate asymptotic behavior. Again these scaling forms are consistent with the critical exponents $`z=2,\alpha =1/2`$ of the Edwards-Wilkinson universality class.

## IV Mean field theory for non-equilibrium wetting

Mean field theories describe a system in an approximate way by ignoring spatial correlations. For a given model there are many possible types of mean field theories, depending on the microscopic level one is trying to describe. Previous mean field approaches to non-equilibrium wetting considered the average height as the dynamical variable and studied the mean field approximation of its dynamics. For example, in a study by Giada and Marsili GiadaMarsili  the KPZ equation (5) was mapped by a Hopf-Cole transformation to a Langevin equation with multiplicative noise, discretizing space and replacing nearest-neighbor interactions by global couplings. Using a Morse potential with a potential well at zero height they were able to reproduce second- and first order transitions as well as phase coexistence. More recently, Santos et al. Santos  extended these studies, suggesting that the mean field phase diagram does not change if fluctuations are taken into account. Moreover, they identified a narrow domain close to the borderline between phase coexistence region and wet phase, where the system exhibits spatio-temporal intermittency.

In this Section we construct mean field equations describing the temporal evolution of the height distribution $`\psi _n`$ for height $`h=n`$. This approach is expected to yield more detailed information on the structure of the interface than previous mean field theories.

### IV.1 Mean field equations

In order to construct the mean field equations for the height distribution, let us first consider the case of a freely evolving interface without a wall. The rate equations describing the temporal change of $`\psi _n`$ consist of several terms corresponding to different processes. Let us, for example, consider the probability that an attempted update leads to the desorption of an atom from the interior of a plateau at level $`n`$. This probability is proportional to $`\psi _n`$, which is the probability to find a randomly selected site at height $`n`$. Moreover, it depends on the heights of the nearest neighbors, which are restricted to take the values $`\{n1,n,n+1\}`$. For simplicity we assume that each site has two nearest neighbor sites. Clearly one can generalize it to the case where each site has $`2(d1)`$ nearest neighbors. In this case one may have several types of edge of corner sites, requiring a larger number of growth rate parameters. To avoid this complication, we restrict ourselves to the case of two nearest neighbors. Thus, neglecting spatial correlations, the probability to find the two nearest neighbor sites at the heights $`l`$,$`m\{n1,n,n+1\}`$ is assumed to be proportional to $`\psi _l\psi _m`$ divided by $`(\psi _{n1}+\psi _n+\psi _{n+1})^2`$. For example, for evaporation from a plateau we have $`n=l=m`$ and thus the contribution of the process (10) to the dynamical equation of $`\psi _n`$ is

$$p\frac{\psi _n^3}{(\psi _{n1}+\psi _n+\psi _{n+1})^2}.$$
(50)
Similarly the loss of probability due to evaporation at the edges (9) is given by

$$\frac{2\psi _n^2\psi _{n1}+\psi _n\psi _{n1}^2}{(\psi _{n1}+\psi _n+\psi _{n+1})^2}$$
(51)
while the depositon process (8) leads to a loss term of the form

$$q\frac{\psi _n^3+2\psi _n^2\psi _{n+1}+\psi _n\psi _{n+1}^2}{(\psi _{n1}+\psi _n+\psi _{n+1})^2}.$$
(52)
In addition, there are corresponding gain contributions at the neighboring levels such that the total probability is conserved. Collecting all terms, the mean field equations can be written as

$$\frac{d\psi _n}{dt}=A_nA_{n1},$$
(53)
where

$$A_n=\frac{p\psi _{n+1}^3+2\psi _n\psi _{n+1}^2+\psi _n^2\psi _{n+1}}{(\psi _n+\psi _{n+1}+\psi _{n+2})^2}q\frac{\psi _n^3+2\psi _n^2\psi _{n+1}+\psi _n\psi _{n+1}^2}{(\psi _{n1}+\psi _n+\psi _{n+1})^2}.$$
(54)
The hard-core wall at the bottom layer can be taken into account by formally setting $`\psi _1=A_1=0`$ and replacing $`q`$ with $`q_0`$ in the expression for $`A_0`$, which then can be written as

$$A_0=\frac{p\psi _1^3+2\psi _0\psi _1^2+\psi _0^2\psi _1}{(\psi _0+\psi _1+\psi _2)^2}q_0\psi _0.$$
(55)
### IV.2 Mean field phase diagram

The main result of the calculations, which will be presented in detail below, is the mean-field phase diagram shown in Fig. 10 for the case $`q_0=q`$. As in the $`1+1`$ dimensional model, there is a continuous phase transition from a bound to a moving phase. In calculating the mean field height distribution in the moving phase it is found that the interface is localized around its average height at any given time, representing a smooth growing interface. This is expected in high dimensions, for which mean field represents a reasonable approximation. Clearly, in $`d=1+1`$ dimensions this is not the case as discussed before.

In the bound region, two types of phases were found. In the larger part of the $`p,q`$-plane (denoted as SE in Fig. 10) the height profile decays superexponentially at large height. In particular it is found that for large $`n`$ the profile takes the form

$$\psi _nq^n\mathrm{exp}(2^{n\alpha }),$$
(56)
where $`\alpha `$ is a constant. In addition, another region is found (denoted as E in Fig. 10), where the height distribution decays exponentially

$$\psi _ne^{an},$$
(57)
where $`a>0`$ is a constant.

The superexponential behavior may be understood by considering the equilibrium case (i.e. $`p=1`$) in which a $`d1`$-dimensional manifold is attracted by a gravitational force to a hard wall. The energy of a fluctuation reaching a height $`n`$ scales as $`n^d`$, leading to the height distribution

$$\psi _ne^{\beta n^d}$$
(58)
where $`\beta >0`$ is a constant. This behavior is seems to be valid in an extended non-equilibrium region of the phase diagram, as shown in Fig. 7, where numerical simulations in $`d=2+1`$ dimensions are presented.

Numerical attempts to find signatures of a possible exponential phase in $`2`$ and $`3`$ dimensions failed since it is very difficult to obtain a reliable statistics in the tail of the height distribution, especially in higher dimensions, where finite size effects become increasingly relevant. The conjecture, that the exponential phase might be related to the roughening transition of KPZ interfaces in $`d>2`$ could not be substantiated by numerical simulations.

### IV.3 Mean field equations

In the stationary state one has $`A_0=0`$ and $`A_nA_{n1}=0`$ so that $`A_n`$ vanishes for all $`n0`$. Therefore, the stationary mean field equations read

$`0`$ $`=`$ $`p\psi _1^3+2\psi _0\psi _1^2+\psi _0^2\psi _1q_0\psi _0(\psi _0+\psi _1+\psi _2)^2,`$ (59)
$`0`$ $`=`$ $`p\psi _{n+1}^3+2\psi _n\psi _{n+1}^2+\psi _n^2\psi _{n+1}q\psi _n(\psi _n+\psi _{n+1})^2\left({\displaystyle \frac{\psi _n+\psi _{n+1}+\psi _{n+2}}{\psi _{n1}+\psi _n+\psi _{n+1}}}\right)^2,`$ (60)
where $`n=1,2,\mathrm{},\mathrm{}`$. In order to solve this equation by iteration, it is convenient to consider quotients of successive densities

$$x_n=\frac{\psi _n}{\psi _{n1}}.$$
(61)
In terms of these variables, the stationary mean field equations take the form

$`0`$ $`=`$ $`px_1^3+2x_1^2+x_1q_0(1+x_1+x_1x_2)^2,`$ (62)
$`0`$ $`=`$ $`px_{n+1}^3+2x_{n+1}^2+x_{n+1}q(1+x_{n+1})^2x_n^2\left({\displaystyle \frac{1+x_{n+1}+x_{n+1}x_{n+2}}{1+x_n+x_nx_{n+1}}}\right)^2.`$ (63)
The bulk equation (63) can be interpreted as an iterative map $`(x_n,x_{n+1})(x_{n+1},x_{n+2})`$, where

$$x_{n+2}f(x_n,x_{n+1})=1\frac{1}{x_{n+1}}+\frac{1+x_n+x_nx_{n+1}}{x_n(1+x_{n+1})}\sqrt{\frac{1+2x_{n+1}+px_{n+1}^2}{qx_{n+1}}}$$
(64)
for $`n=1,2,\mathrm{}\mathrm{}`$. In addition the initial condition for $`x_1,x_2`$ is given by the bottom layer equation

$$x_2=1\frac{1}{x_1}+\sqrt{\frac{1+x_1(2+px_1)}{q_0x_1}}.$$
(65)
### IV.4 Stationary solutions

In order to evaluate the stationary height distribution for given $`p`$ and $`q`$, one has to look for a point on the line of possible initial conditions (65) and iterate the map (64) such that it reaches a real fixed point with $`x^{}<1`$. This trajectory then corresponds to a physical height distribution.

To proceed, we first analyze the fixed points of the map. The fixed point equation

$$x\left(qx^3+(2qp)x^2+(q2)x1\right)=0$$
(66)
has four solutions. Two of the solutions, $`x_0^{}=0`$ and $`x_1^{}`$ are real in the entire $`p,q`$-plane, while other two solutions $`x_2^{},x_3^{}`$ may either be both real or complex conjugate to each other. We denote the region in the $`p,q`$-plane, where $`x_2^{},x_3^{}`$ are real, by region III (see Fig. 11). We further divide the complementary region to III into two regions, I where $`x_1^{}<1`$, and II, where $`x_1^{}>1`$. Since the bracket in Eq. (66) is equal to $`1`$ at $`x=0`$ and positive for $`q>0`$ in the limit of large $`x`$ the fixed point $`x_1^{}`$ has to be positive.

Analyzing the map (64) with the initial condition (65) we find that the only fixed points, which correspond to physical height distributions, are either $`x_0^{}`$ or $`x_2^{}`$ when it is real. Here $`x_2^{}`$ is the solution of Eq. (66) which satisfies $`x_2^{}<1<x_3^{}`$ in the region III. In particular we find that below the transition line, the height distribution is controlled by the fixed point $`x_0^{}=0`$ in region I and by the fixed point $`x_2^{}`$ in region III. Details of the analysis, which led to this result, are given in the Appendix.

We now study the stationary height profile in the two regions. In region I the map flows to the hyperbolic fixed point $`x_0^{}=0`$ along its stable trajectory. Expanding the map for small values of $`x`$ to lowest order, this stable manifold is described by the nonlinear relation

$$x_{n+1}qx_n^2,$$
(67)
in the limit $`x0`$. Along this manifold, the map approaches the fixed point $`x_0^{}=0`$ super-exponentially as

$$x_n\frac{1}{q}\mathrm{exp}(2^{n\alpha })$$
(68)
yielding the height profile in Eq. (56).

In region III stationary solutions are controlled by the fixed point $`x_2^{}>0`$ (see Appendix). Linearizing the map around this fixed point one obtains an exponential behavior

$$x_ne^{an},$$
(69)
where $`a>0`$ is a constant, leading to the exponential height profile (57).

To complete the analysis of the phase diagram, one has to locate the wetting transition line. This is done by simulating the dynamical mean field equations (53) and determining the point from where on the interface detaches. This analysis leads to the transition line shown in Fig. 10.

So far all mean field results were obtained for $`q_0=q`$, i.e., without an attractive force between the substrate and the wetting layer. Lowering $`q_0`$ changes the bottom layer equation (65) and thereby the possible starting points of the iteration. As shown in the Appendix, the mean field approximation reproduces the phenomenon of phase coexistence. However, unlike in the model in $`1+1`$ dimensions, it emerges everywhere along the phase transition line as soon as $`q_0<q`$; hence, there is no tricritical point. Moreover, in the limit $`q_00`$ the threshold for the growth rate $`q`$, where the interface detaches, tends to infinity. Therefore, the region of phase coexistence is not bounded from above as in the $`1+1`$ dimensional model.

## V Conclusions

In this paper we have given a detailed account of a recently introduced solid-on-solid model for non-equilbrium wetting, which is defined in the spirit of a KPZ equation in a potential. Introducing a hard-core wall at zero height, the model exhibits a contininuous wetting transition from a bound to a moving phase. The model is controlled by two paramters $`p`$ and $`q`$, which effectively determine the asympotic slope of the potential and the coefficient of the nonlinear term in the KPZ equation.

For $`p=1`$ the dynamical rules of the model obey detailed balance. In this case the stationary distribution of a bound interface is given by a Boltzmann ensemble, which allows one to derive various quantities exactly. Moving away from this line, the model crosses over to a non-equilibrium behavior, that is characterized by different critical properties. The model can be generalized further by including an attractive interaction between the substrate and the wetting layer. If this force is strong enough it may turn the continuous transition into a discontinuous one. Moreover, the bound and the moving phase may coexist in regions where the coeffient of the nonlinear term in the KPZ equation is negative.

In order to assess the behavior of the model in higher dimensions, we have proposed a mean field approximation, which is based on rate equations for the densities at different heights. The phase diagram turns out to be surprisingly rich. It turns out that the mean field approximation reproduces the properties of the original model. However, in the moving phase the interface remains smooth as it is expected for KPZ-type growth in higher dimensions.

As a new feature, the mean field approximation predicts the existence of two different regions in the bound phase. In one of these regions the interace profile decays superexponentially with increasing height while in the other region an exponential decay is observed. This can be explained by classifying the fixed points of an iterative map for quotients of the densities. The question to what extent this crossover from superexponential to exponential profiles can be observed in the full model is still open.

Acknowledgments
The support of the Israel Science Foundation (ISF) and the Einstein center of the Minerva foundation is gratefully acknowledged. The simulations were partly performed on the ALiCE parallel computer at the IAI in Wuppertal. H.H would like to thank the INFN and the Weizmann Institute for hospitality, where parts of this work have been done.

## Appendix A Stationary solutions of the mean-field equations

### A.1 Fixed points

Besides $`x_0^{}=0`$, Eq. (66) has three other fixed points which are the roots of the polynomial equation

$$qx^3+(2qp)x^2+(q2)x1=0.$$
(70)
As shown in Fig. 11, the $`p,q`$-plane can be divided into three different regions. In regions I and II one fixed point is real and two of them are complex conjugate, while in region III all fixed points are real. The boundary of region III is characterized by the existence of a two-fold degenerate fixed point. Comparing Eq. (70) with a polynomial of the form $`(xa)^2(xb)=0`$, the boundary of region III can be given in a parameter representation by

$$q=\frac{2}{aa^2},p=\frac{1+3a+4a^2}{a^2a^3},b=\frac{1a}{2a},$$
(71)
where $`0<a<1`$. The boundary has the form of a wedge, shown as a dashed line in Fig. 11. The lower branch for $`0<a<1/2`$ and the upper branch for $`1/2<a<1`$ terminate in the triple point $`(p,q)=(8,28)`$, where Eq. (70) has a three-fold degenerate fixed point $`x^{}=1/2`$.

Physically meaningful stationary solutions of the iterative map must start from a point on the line given by Eq. (62) and have to flow towards a real fixed point with $`0x^{}<1`$. For this reason we divided the complement of region III into two parts, namely, region I with $`x_1^{}>1`$ and region II with $`x_1^{}<1`$. Both regions are separated by a straight line $`q=(p+3)/4`$, where $`x_1^{}=1`$.

### A.2 Stationary solutions in regions I and II

We first show that in regions I and II a physically meaningful stationary solution always flows to the fixed point $`x_0^{}=0`$. To this end we show that the other real fixed point $`x_1^{}`$ is either larger than $`1`$ or unstable.

As shown in the upper inset of Fig. 11, the phase transition line starts at $`p=0`$ in region I, crosses into region II at the point $`q=p=1`$ and then crosses back into region I at the point $`(p,q)(5.380,2.095)`$. Obviously the fixed point $`x_1^{}`$ can only be physically meaningful between these two crossing points, where $`x_1^{}<1`$. However, in this interval $`x_1^{}`$ turns out to be unstable. To demonstrate this point we consider the eigenvalues of the Jacobian of the map $`x_{n+2}=f(x_n,x_{n+1})`$ in Eq. (64)

$$\lambda _{1,2}=\frac{1}{2}\left(\frac{f}{x_{n+1}}\pm \sqrt{\left(\frac{f}{x_{n+1}}\right)^2+4\frac{f}{x_n}}\right)|_{x_n=x_{n+1}=x^{}}.$$
(72)
Using Eq. (70) the partial derivatives can be expressed as

$$\begin{array}{cc}\hfill \frac{f}{x_n}|_{x_n=x_{n+1}=x^{}}& =\frac{1}{x_{}^{}{}_{}{}^{2}},\hfill \\ \hfill \frac{f}{x_{n+1}}|_{x_n=x_{n+1}=x^{}}& =\frac{qx_{}^{}{}_{}{}^{4}+2qx_{}^{}{}_{}{}^{3}+2(q1)x_{}^{}{}_{}{}^{2}+(3q2)x^{}2}{2qx_{}^{}{}_{}{}^{3}(x^{}+1)}.\hfill \end{array}$$
(73)
Along the dotted line in Fig. 11, where $`x^{}=1`$, the two eigenvalues

$$\lambda _{1,2}|_{x^{}=1}=\mathrm{\hspace{0.33em}1}\frac{1}{4q}(3\pm \sqrt{924q})$$
(74)
are complex conjugate in the interval between the two crossing points $`1<q<5.380`$. As can be verified numerically, they are also complex conjugate in a neighborhood of this line, which includes the phase transition line. Since the determinant of the Jacobian

$$\lambda _1\lambda _2=\frac{1}{x_{}^{}{}_{}{}^{2}}$$
(75)
is larger than $`1`$ for $`x^{}<1`$, the fixed point $`x_1^{}`$ is found to be unstable. Thus we can conclude that physically meaningful stationary solutions in regions I and II are always controlled by the fixed point $`x_0^{}=0`$.

Since for $`x^{}0`$ the two eigenvalues tend to $`\lambda _1=0`$ and $`\lambda _2=\mathrm{}`$, the fixed point $`x_0^{}=0`$ is nonlinear and hyperbolic. It has superstable manifold which in the limit of $`x0`$ is given by $`x_{n+1}=qx_n^2`$. This manifold is shown in Fig. 12 for different values of $`q`$ below, at, and above the critical point. As it can be seen, a stationary solution exists if the manifold intersects the line of possible initial conditions (65).

### A.3 Stationary solutions in region III

In the wedge-shaped region III, there are three real fixed points $`x_1^{}<x_2^{}<x_3^{}`$. The first one is smaller than $`1`$ and unstable, while the second one is smaller than $`1`$ and hyperbolic. The properties of $`x_3^{}`$ depend on $`p`$ and $`q`$. Below the dotted line in Fig. 11 it is larger than $`1`$ and stable, while it is smaller than $`1`$ and unstable above.

Fig. 13 illustrates typical situations at $`p=60`$ for various values of $`q`$, crossing the region III from bottom to top. Below the wedge in region I (panel A), the stationary solution can be constructed in the same way as before. Entering region III from below (panel B), two new real fixed points $`x_1^{}<x_2^{}<1`$ emerge, which are fully unstable and hyperbolic, respectively. The superstable manifold of the fixed point $`x_0^{}=0`$ now originates in $`x_1^{}`$, while it is the linearly stable manifold of the hyperbolic fixed point $`x_2^{}`$ which intersects the dashed line of possible initial conditions (65). As before, the intersection point moves continuously to infinity as $`q`$ approaches the critical point (panel C). Thus the unbinding transition manifests itself in the same way as in the regions I and II, the only difference being that the physically relevant stable manifold is now controlled by the hyperbolic fixed point $`x_2>0`$ instead of $`x_0^{}=0`$. The linear stability of $`x_2^{}`$ is responsible for the purely exponential profile observed in region $`E`$.

### A.4 Phase coexistence in regions I and II

Lowering $`q_0`$ changes the line of possible initial conditions (65), while the stable manifold of the fixed point $`x_0^{}=0`$ remains the same. The typical situation for $`p=0.8`$ is shown in Fig. 14. The left panel shows the stable manifold (solid line) and the bottom layer equation (dashed line) without attractive force at the critical point $`q_0=q=q_c0.943`$. Both curves approach each other smoothly and intersect at infinity so that the transition is continuous. The right panel shows the same situation in the presence of an attractive force for $`q_0=0.5`$. Accordingly, the effect of lowering $`q_0`$ is to move the dashed line of possible intial conditions upward such that it intersects the critical stable manifold at a certain finite point. This means that the bottom layer density $`\psi _0`$ is finite at the critical point, making the transition first order. Moreover, a stationary solution still exists even if $`q`$ is increased, proving the possibility of phase coexistence in the mean field equations. Increasing $`q`$ beyond a certain threshold, a stationary solution no longer exists. This defines the upper boundary of the phase coexistence region in the mean field phase diagram.