Critical and Tricritical Wetting in the Two-Dimensional Blume–Capel model

The two-dimensional Blume–Capel model with free surfaces where a surface field H1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$H_1$$\end{document} acts and the “crystal field” (controlling the density of the vacancies) takes a value Ds\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$D _s$$\end{document} different from the value D\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$D$$\end{document} in the bulk, is studied by Monte Carlo methods. Using a recently developed finite size scaling method that studies thin films in a L×M\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$L \times M$$\end{document} geometry with antisymmetric surface fields (HL=-H1)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(H_L=-H_1)$$\end{document} and keeps a generalized aspect ratio c=L2/M\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$c = L^2/M$$\end{document} constant, surface phase diagrams are computed for several typical choices of the parameters. It is shown that both second order and first order wetting transitions occur, separated by tricritical wetting behavior. The special role of vacancies near the surface is investigated in detail.


Introduction
In recent years the development of nanoscience has led to new techniques to produce special patterns of nanostructured surface layers and devices in the nanoscopic size range [1,2], where surface and interfacial effects, such as wetting phenomena [3][4][5][6][7][8][9][10][11] may play an important role, and hence have become a subject of intense research of statistical physics.While for the twodimensional Ising model with nearest neighbor exchange J and a free surface at which a boundary field H 1 acts an exact solution [4,[12][13][14] is available and has provided very valuable insight, no exact solutions are available for extensions such as wetting transitions in the Blume-Capel model [15,16] and then the use of Monte Carlo simulation methods [17,18] is appropriate [19].This extension is of interest, since depending on the choice of parameters (exchange constant J and "crystal field" D that controls the density of vacancies in this model) its order-disorder transition may be of first or second order, and a bulk tricritical point [20] occurs.Applying an extension [19] of finite size scaling methods [21] that is appropriate for the anisotropic critical phenomena associated with continuous wetting transitions, the wetting phase diagram near the tricritical point in the bulk has been studied [19].
In the present work, we extend this study by allowing for the possibility that the "crystal field" in the surface layer differs from the bulk, D s = D.Such a choice is very natural physically, since surface sites of the square lattice no longer experience the full lattice symmetry.We show that then the possibility of first order wetting transitions and wetting tricritical points [5,6,8] arises.Note that mean field theory of wetting with short range surface forces [5][6][7][8]22,23] shows that both first order and second order wetting transitions may result.Of course, in d = 2 dimensions Ginzburg-Landau type [22,23] theories are very unreliable [24], due to their neglect of statistical fluctuations; but the tricritical wetting behavior in d = 2 beyond mean field was only occasionally considered [25].
The outline of the present paper is as follows: in Sect. 2 we introduce the model, briefly summarize the theoretical background and the simulation methodology.Section 3 gives our numerical results, while Sect. 4 summarizes our conclusions.

The Blume-Capel Model in Thin Film Geometry with Free Surfaces and Competing Surface Fields
As is well-known, the Blume-Capel model [15,16] is a lattice model where spins S i reside on lattice sites i and may take three states, S i = 0 and S i = ±1.In the bulk, the Hamiltonian is where we assume a ferromagnetic exchange J > 0, the summation i, j runs over all pairs of nearest neighbors on the square lattice once, and the "crystal field" D controls the average density of the vacancies (i.e., sites with S i = 0).For D → −∞ vacancies are suppressed, and the model reduces to the exactly solved Ising model [26].In the absence of magnetic fields and with a constant crystal field D, which leads to a homogeneous distribution of nonmagnetic impurities, the two-dimensional Blume Capel model exhibits a nontrivial bulk phase diagram with a tricritical point [27].For the square lattice, this tricritical point occurs at D t /J = 1.965 and k B T t /J = 0.609 [27,28].For D > D t the transition at T cb (D) is of first order, while for D < D t the transition becomes of second order.For studying phase transitions in the bulk, one needs to consider the limit where the system is infinite in both lattice directions.However, when one wishes to study wetting transitions, the proper geometry is a "semi-infinite" system (Fig. 1), i.e. we have a free surface at the first row (extending in x−direction) that is located at y = 1 and is labeled by the index Fig. 1 Sketch of the geometry of the model and its state at a temperature T slightly below the temperature T w (H 1 ) where in the thermodynamic limit critical wetting transitions occur.The irregular line indicates a configuration of the fluctuating interface, separating a domain of the spin pointing up phase (above the interface, indicated by the double arrow pointing upwards) from a domain of spin down phase (double arrow pointing downwards).For T ≤ T w (H 1 ) the interface still is (weakly) bound to one of the walls (which are shown shaded); without loss of generality we have assumed a state where the interface is bound to the lower wall, so in this state the magnetization m of the L × M strip is positive.Along the x−direction, periodic boundary conditions are used, while in y−direction the walls are represented by missing spins in the rows n = 0, n = L + 1, while boundary fields −|H 1 | and +|H 1 | act at the rows n = 1 and n = L, respectively.In addition, at both boundary rows the surface crystal field (D s ) differs from its bulk value (D).The interfacial fluctuations are characterized by anisotropic fluctuations, the correlation length ξ || along the interface is much larger than the correlation length ξ ⊥ perpendicular to it.Local fluctuations of the magnetization inside the domain are not shown n = 1.However, systems accessible to computer simulation must be finite in all directions [17,18].So the linear dimension in the perpendicular x−direction is taken to be M lattice spacings, (the lattice spacing is put to unity henceforth), and periodic boundary conditions in x−direction are used.At this first (n = 1) boundary row a surface magnetic field H 1 acts, and the crystal field parameter D s is assumed to differ from its bulk value D (unlike our previous work [19]).Now the linear dimension L in y−direction also must be finite, and the crucial question is what boundary condition is used at the last row (n = L).In early simulations of wetting transitions in the d = 2 and d = 3 Ising model both the lower and the upper boundary were treated fully symmetrically [29][30][31].However, this choice has the disadvantage that for any finite value of L phase coexistence between the "spin up" and "spin down" phases no longer occurs for zero bulk field, but at a shifted value (this effect is nothing but the well-known capillary condensation [8,32]).This shift of "bulk" coexistence is avoided when one chooses the boundary condition of strictly antisymmetric surface fields [33].Thus, the Hamiltonian becomes (2) (Fig. 1).Note that it is natural to choose the "surface crystal field" symmetric, D 1 = D L ≡ D s .

Brief Details of the Monte Carlo Simulation Method
Monte Carlo simulations are performed by using the standard Metropolis algorithm (see e.g.[17] for a review).Typical runs are performed over a length of 10 7 Monte Carlo steps per lattice site (MCS), disregarding the first 2×10 6 MCS to allow the system to reach equilibrium.Note that for systems far below bulk criticality exposed to boundary fields cluster algorithms do not present any advantage [34].The square lattice with a rectangular geometry of size L × M is used in the simulations.According to the description of Fig. 1 we assumed free boundary conditions along the M−direction (horizontal) where the surface fields are applied, while periodic boundary conditions are taken along the perpendicular (vertical) direction.Along the simulations we recorded the magnetization m per lattice site (N = L M), given by and we compute its thermal expectation value m T .Furthermore, relevant moments of the nagnetization, i.e. m 2 T , m 4 T , and the cumulant With the statisitcal effort quoted above, the relative error of our estimates typically is of the order of 2 % (or less ).Therefore error bars in the figures that will follow typically are smaller than the size of the symbols, and therefore not shown.More details of the simulation method used in this work can be found in reference [19].

Finite Size Scaling Methodology
Now in order to locate any phase transition of a model system within the framework of Monte Carlo simulations, it is necessary to carry out an extrapolation to the thermodynamic limit (L → ∞, M → ∞) via a finite size scaling analysis [17][18][19]21,35,36].For this L × M, geometry, normally there are three ways considered to take the thermodynamic limit: (i) Keeping the aspect ratio M/L of the system constant [37].This is the appropriate choice if one wishes to study the bulk phase transition of the model (which belongs to the universality class of the Ising model, of course, where correlations are isotropic at criticality [38]).(ii) Keeping L finite and taking the limit M → ∞.However, this system is quasi-onedimensional and can exhibit a phase transition at T = 0 only [39][40][41].(iii) Keeping the generalized aspect ratio c = L 2 /M constant [19].As discussed in detail in our previous paper [19], this is the appropriate choice to study critical wetting in d = 2, since the correlation lengths ξ || , ξ ⊥ of interfacial fluctuations (Fig. 1) are known [4][5][6]25] to diverge with different critical exponents ν || , ν ⊥ , namely where t is the reduced distance in temperature (t = T /T w (H and U = f 4 (Lt, c), (7) where f 1 , f 2 and f 4 are scaling functions, respectively.From Eqs. 5-7 it follows that when one considers any of these quantities as function of t, one obtains for t = 0 just constants (which depend on c only, but no longer on L).As a consequence, by plotting m , m 2 and U vs. t and several choices of L (or M, which is equivalent since using L = √ Mc one can always substitute M for L), these curves must intersect at t = 0; yielding an "intersection method" to locate where the critical wetting transition occurs.It is worth to recall that the method works for physical observables that are increasing (decreasing) functions of the sample size for t < 0 (t > 0), or viceversa, in the neighborhood of the critical point, as e.g. the quantities given by Eqs. ( 5)- (7).This method of finite size scaling analysis has also been used successfully for a slightly different model, where instead of boundary fields a fixed spin boundary condition was used [43], as well as for the case of inhomogenoeus short-range surface fields [44].The intersection method gives rather accurate critical points for modest lattice sizes, as it follows e.g. by comparing exact results derived by Abraham [12] for the critical wetting of the Ising model (Eq.( 2) with D = −∞) with simulation results [19,42].Based on this observation and the fact that the aim of the manuscript is to provide a general overview of the critical and tricritical wetting of the Blume-Capel model in a four dimensional space of parameters, rather than to determine accurate locations of critical points, we have made no attempts in order to either introduce scaling corrections to Eqs. ( 5)- (7) or to perform extrapolations to the limit of large lattice sizes.

Ground State Surface Phase Diagram
At zero temperature, a semi-infinite (L → ∞) Blume-Capel model where at the free surface (the row n = 1 in Eq. 2) the surface field H 1 and the surface crystal field D s acts, may be in three different surface phases.Without loss of generality, but for D/J < 2, we assume that the bulk of the system is in the state with spins S i = −1.So the surface phases are: (i) The surface is nonwet (and ordered).Having assumed that in the bulk the magnetization is negative the layer magnetization (m n = −1 for large layer index n), in the ground state then actually m n = −1 for all n, including the layer n = 1.(ii) The surface is wet (and ordered).Then m 1 = +1 although m n = −1 in the bulk.Since in the ground state there is a single straight interface in the system, separating rows of up spins oriented parallel to the free surface from the rows of down spins that extend into the bulk, there is a degeneracy with respect to : for the model of Eq. 2, all choices = 1, 2, 3, . . .are energetically degenerate.Of course, at nonzero temperature entropy requires → ∞ in the wet state, in order to avoid any constraint on the interfacial fluctuations.However, at T = 0 the surface excess energy of the wet state of the surface can be computed from the simple case = 1, a single row of up spins is separated by an interface (that involves an energy cost of 2J per spin) from the domain of the down spins.(iii) The surface is nonwet and disordered, all spins in the layer n = 1 take the value S i = 0, so H 1 has no effect.In this case there is no degeneracy, the interface must be between layer 1 (where m 1 = S i i∈n=1 = 0 and layer 2 (where m 2 = −1, like in all layers n with n ≥ 2).This "nonwet disordered state" of the surface has previously been identified as a "prewetting" ground state by Fytas and Selke [43] in studies of the wetting behav-ior of the Blume-Capel model where, instead of considering antisymmetric boundary fields, spins pointing up and down are fixed at opposite walls where a reduced coupling constant (α) is assumed.Under these assumptions and for D/J ≥ 2 − α [43] there is also a single line of spins S i = 0 at the surface layer.On the other hand, the prewetting of a single layer of zero spins included between a microscopic interface of +/− spins has also been discussed by Hryniv et al. [45] in a study of the surface tension of the Blume-Capel model at low temperatures (see also reference [46].) Note that in the groundstate of the model there cannot be any vacancies (for D < 2J since at D = 2J the transition in the bulk occurs) in any interior row.While the surface field H 1 can create a domain of up spins that is rows thick (with arbitrarily large) separated from the bulk phase where in the groundstate all spins S i = −1, the surface crystal field D s can only stabilize a single row (n = 1) that has S i = 0, it can not induce a "wetting-like" layer of multiple rows with S i = 0 at T = 0. Thus, there is a qualitative difference between the effects that the surface field H 1 and surface crystal field D s can cause.
We now establish the ground state phase diagram to clarify where the three surface states described above are stable.Noting from Eq. 1 that in the bulk the energy per spin is U b = −2J + D, we find by a "bond-counting" argument easily the surface excess energies of the three surface states nonwet and ordered (i), which we will denote by the index "nwo", wet (ii), and nonwet and disordered (iii), denoted as "nwd": and The total surface excess internal energies (due to both surfaces) are defined simply from the expectation value of the total Hamiltonian per spin, subtracting the bulk energy per spin, and multiplying by the number L of layers, for proper normalization.Eqs.( 8)-( 10) refer to the surface excess energy at T = 0 due to a single surface only.Surface excess free energies are defined analogously.The transitions between the above discussed states then occur at fields H 1c that are simply found from equating these energies, and wet → nwd : These results are summarized in the groundstate phase diagram (Fig. 2).Note that we predict that at nonzero temperature the transition nwo-nwd is always rounded.The transitions wet-nwo and wet-nwd then join into a single transition near the special point H 1c /J = 1, D s /J = 1, where all transition lines at T = 0 meet.Of course, at finite temperature the two phases, nwo and nwd, which the wet state at T = 0 transforms to, are no longer distinct, but part of the same phase, whose characteristics gradually change when the line nwo-nwd (where only at T = 0 a transition occurs) is crossed.Therefore we predict on the basis of these qualitative considerations, that at constant but low temperature the surface phase diagram

Monte Carlo Results and Their Analysis
In our previous work [19] we have studied the present model for the special case D s = D, and have found that for the regime D/J < D tri /J , where the phase transition in the bulk is of second order (while at the tricritical point D/J = D tri /J ≈ 1.965(5) the order of the transition changes from second order to first order [20]), the wetting transition is always continuous.For instance, for the special case D/J = 0 it was found that for the choice of H 1 /J = 0.55 this critical wetting transition occurred for k B T w (H 1 )/J = 1.393 ± 0.004.
We now ask what is the effect of choosing now D s = D in such a situation, leaving all other parameters unchanged.As an example, Fig. 3 presents an extreme case, D s /J = −19.65.This is close to the Ising limit, for D s → −∞ vacancies at the walls would be completely suppressed.Therefore the boundary excess free energies of both walls get somewhat larger in magnitude (the chosen field H 1 /J = 0.55 leads to a larger wall excess free energy, since it has an effect at practically all lattice sites adjacent to the wall, while for D s /J = 0 some of these sites will be taken by vacancies, and then on these sites the surface field has no effect).Since the location of the critical wetting transition depends on the wall excess free energy difference f between semi-infinite systems with positive (+) and negative (−) spontaneous magnetization (both exposed to a positive boundary field (H 1 /J = 0.55 in the present example)), we expect now the transition to occur at a somewhat lower temperature: wetting occurs, according to Young's criterion [3][4][5][6][7][8][9][10][11], when f is the interfacial tension between bulk coexisting phases.When we increase f s , as we do by choosing D s strongly negative, wetting can only occur at a lower temperature, since f int (T ) (that monotonically grows from f int (T c ) = 0 to larger values as T is lowered) must become correspondingly larger, which happens only at a lower temperature.Indeed, the data of Fig. 3 indicate that k B T w (H 1 = 0.55J )/J = 1.191(4), consistent with the predicted trend.As it should be, the three sets of curves exhibit common intersection points at the same An unexpected feature, however, is the observation that the critical values |m| crit , m 2 crit , and U crit where these intersections occur for D s < D are distinctly larger than their counterparts for D s = D (see [19]).This finding illustrates once more the important fact that in finite size scaling there occurs much less universality than one might have expected: the finite size scaling "invariants" f 1 (o, c), f 2 (0, c) , and f 4 (0, c) (Eqs.5 -7) at the transition not only depend on the generalized aspect ratio c, as the notation of Eqs.5-7 suggests, but they do also depend strongly on the details of the boundary conditions used.Different choices of D s /J then must be considered as different boundary conditions at criticality.
We have studied the wetting transition for this choice of D s /J also for many other choices of D/J .While for a range of values of D/J roughly when D/J < 1 (e.g. for H 1 /J = 0.55) the behavior is qualitatively similar as in the case shown in Fig. 3, for larger values of D/J an unexpected and qualitatively different behavior occurs, namely discontinuous wetting transitions.An example for a wetting transition that is strongly of first order is shown in Fig. 4a.Obviously the data are not very sensitive to the chosen initial condition, so there is little hysteresis, and one can locate the transition rather precisely (T w (H 1 )/J ≈ 0.585 in Fig. 4).Furthermore, the absence of strong finite size effects is also observed.On the other hand, first order wetting transitions can also be obtained by scanning different set of parameters, as e.g. in plots of |m| T versus the surface field (cf.Fig. 4b), see also below the phase diagram shown in Fig. 10b).Note that the rather large value of |m| T in the wet phase (characterized by a fluctuating interface in the center of the strip, see Fig. 5) imply that in the wet phase huge capillary-wave type fluctuations of the delocalized interface in the center of the strip occur.Due to our choice of elongated systems with constant aspect ratio (L 2 /M = 9/8) this situation of complete wetting is a kind of critical state, and thus we do not find that |m| T tends to zero with increasing L in this region.
The snapshot configurations of Fig. 5 show the abrupt change in the composition of the system due to a tiny variation of the temperature, as expected for first order transitions (compare Fig. 5a, b).Also, the enrichment of vacancies at the interface already reported for critical wetting [19,43] is observed here for complete wetting.In a recent paper Fytas and Selke [43] have provided evidence on the critical behavior of the interfacial adsorption of vacancies, which follows from a careful study of the temperature derivative of the excess density of vacancies caused by the presence of the wetting interface.Of course, it would be worthwhile to analyze the possible onset of criticality upon interfacial adsorption in the case of complete wetting reported in this paper, however this task requires a huge computational effort [43], so that it is beyond the aim of the present paper that was aimed to explore the occurrence of tricritical wetting in a four dimensional space of parameters of the Blume Capel model.Spins pointing up are shown in blue (dark-grey when the figure is displayed in black/white tones), spins pointing down in white, and spins S i = 0 in red color (gray when the figure is displayed in black/white tones), respectively.Three temperatures are shown, k B T /J = 0.579 (left), 0.591 (middle), and 0.603 (right).Note that vacancies occur in both "spin up" and "spin down" phases in the bulk, and are strongly enriched at the fluctuating, meandering interface, while they are completely expelled from the surface (Color figure online) Figure 6 shows then the location of these wetting transitions in relation to the bulk phase diagram, and to the line of critical wetting transitions that occur for D s = D (taken from [19]).Interestingly, the present model where vacancies are expelled from the boundaries (Fig. 5) exhibits a case of two-dimensional tricritical wetting.To our knowledge, this is the first time that a specific model with this behavior is discussed; of course, from the general renormalization group theory of wetting in d = 2 dimensions the possibility of tricritical wetting has already been pointed out [25], and the tricritical exponents were proposed.While for D s = D the wetting transition stays second order until the line T w (H 1 )/J hits the bulk first order phase boundary T cb (D)/J where it then ends, for at least sufficiently negative values of D s as shown in Fig. 6 a different behavior occurs.However, a systematic variation of D s and D to clarify when tricritical wetting transitions first appear is beyond the scope of the present paper.
The fact that reducing the vacancy concentration near the boundaries is enough to turn at low temperatures the wetting transition from second order to first order certainly is unexpected: a naive expectation would be the argument that removing vacancies from the boundaries just "renormalizes" the boundary field H 1 to a larger value (since at sites taken by vacancies, occurring for D s = D also at the boundaries, the boundary field has no effect).Thus it is of interest to study the effect of the lack of vacancies in the boundary rows in more detail when one changes D/J at fixed H 1 /J and D s /J , comparing profiles of the magnetization and density of vacancies across the film, when one moves (in the wet phase where the interface is centered in the middle of the strip) from the regime where second order transitions occur (Fig. 7a, b).One can see that in the second order case the magnetization profiles show an almost linear variation across the strip, while in the first order case (Fig. 7c, d) a pronounced S-shape occurs; while in the second order-case the vacancy distribution (distinctly nonzero vacancy densities are only found in the rows n = 2 to n = L − 1, while in the rows 1, L where the surface crystal field acts the vacancy density is essentially zero) shows only a small dependence on D. In the range where first order transition occurs vacancies are much more localized in the middle of the strip, reflecting the interfacial enrichment of vacancies also seen in the snapshot (Fig. 5).
For completeness, we briefly discuss the behavior when D s /J is positive and sufficiently large.Already from the analysis of the groundstate of the model, Fig. 2, one has to expect that no wetting transition can occur, and this expectation is confirmed by the simulations.Figure 8a and b give typical examples of the behavior of the profiles of magnetization and vacancy density when the temperature is scanned through the transition point of the bulk, respectively.Irrespective of temperature, the rows n = 1 and n = L for large D/J are exclusively taken by vacancies, and hence have zero magnetization.Effectively, we have a strip of linear dimension L − 2 across the system, neither the value of H 1 /J nor the precise value of D s /J matter.The rows n = 2, n = L − 1 behave similarly as in a system with free surfaces and no surface field and no surface crystal field.Hence the magnetization in these rows is reduced in comparison with the bulk, and the density of vacancies is enhanced.Since the data shown in Fig. 8a are taken close to the bulk critical point at all temperatures the correlation length is rather large in all cases, and hence a completely flat behavior of the profiles is only reached deeply in the interior of the strip.Note, however, that the data of Fig. 8a refer to a rather elongated system, L × M = 36 × 576.For such elongated geometries it is well known that the system at temperatures slightly below criticality typically does not develop a monodomain configuration, but rather several domains with opposite magnetization separated by boundaries that are oriented perpendicular to the external wall exist.These domains are readily recognized when one examines configuration snapnshots (not shown here since they are very similar to the corresponding ones in the Ising model [37]).These domain states are responsible for significant finite size effects (Fig. 8c) at temperatures somewhat lower than the critical temperature.Nevertheless a finite size scaling analysis of the phase transition from the ordered to the disordered phase is fairly straightforward: critical correlations grow isotropically as the transition is approached, irrespective of system geometry one has for the correlation length ξ b ∝ |1 − T /T cb | −ν with ν = 1 the standard Ising model critical exponent [26].In fact, the universality principle [38] puts the Blume-Capel model in the same universality class as the Ising model, and also the exponent of the spontaneous magnetization, β = 1/8 [26] (recall that m ∝ (1 − T /T cb ) β in the thermodynamic limit).Since the correlations grow isotropically, and Eq. 4 does not apply here, the standard aspect ratio L/M (and not the generalized aspect ratio L 2 /M appropriate for a study of wetting transitions in d = 2 dimensions) needs to be kept constant in the finite size scaling "data collapsing" analysis (Fig. 8).Indeed, one finds that the data points fall on a master curve reasonably well.Only for |m| L β/ν > 1, however, is this master curve the expected straight line on the log-log plot (necessary to recover the exponent β, as quoted above).The rapid fall-off of the scaling function for |m| L β/ν < 1 (which then stabilizes for T = T cb at a rather small value of about 0.15) reflects the breakup of the uniform magnetization into a state with at least 2 domains (for number of domains for L/M = 1/16 actually is 4); due to the periodic boundary condition in the x−direction along the free boundaries, only domain states with an even number of domains can occur.The situation depicted in Fig. 8 is qualitatively similar to early studies of the Ising model with periodic boundary conditions [37] (although there the magnetization profile across the strip is strictly constant) and to studies of capillary condensation [47,48] (where one uses symmetric rather than antisymmetric boundary fields).We emphasize that in the present work we use antisymmetric boundary fields also for the profiles shown in Fig. 8, but a symmetric magnetization profile nevertheless results there because the large vacancy density in the rows n = 1 and n = L (that is practically indistinguishable from its saturation value, unity) means that the effects of the surface fields are effectively removed completely.Of course, Fig. 8 shows an extreme case, where by choosing a large positive value of D s /J we have canceled the effect of the surface almost completely, and the magnetization profile becomes practically strictly symmetric, m(n) = m(L + 1 − n), n = 1, 2, . . ., but when we choose |D s |/J not so large, the vacancy density at the walls does not take its extreme values zero (when D s is negative but |D s |/J 1) or unity (when D s /J 1), but rather takes on intermediate values, and then an interesting interplay with the surface field is possible (Fig. 9).Consider for instance a state with negative magnetization in the bulk, for conditions where the system is in the non-wet state.For the same choice of H 1 /J < 0, the magnetization in row n = 1 can be more negative (for D s /J = −3.93 and −1.965, in the example shown in Fig. 9a) or less negative (for D s /J = −0.9825and −0.4913, cf Fig. 9b) than in the bulk, while at the other boundary of the strip (n = 24 in Fig. 9a, b), where the sign of the boundary field is positive, the magnetization always is distinctly less negative, for the cases shown.In the examples shown in Fig. 9a and b, the symmetry between m(n) and m(L + 1 − n) that was present in Fig. 8a clearly is lost, as expected.Also the density profile of the vacancies no longer exhibits an analogous symmetry strictly (Fig. 9b, c).It is clearly seen that the density of vacancies in the boundary layers is very small for strongly negative values of D s , and then it is also smaller than in the bulk (for the example shown here; for D s /J = −1.4738and D/J = 0 the vacancy density shows some enhancement already in the rows n = 2 and n = L − 1 adjacent to the boundary rows, however).But then occurs a relatively small range of values for D s (still for D s < 0) where the vacancy density increases.For the sake of completeness we note that for a large positive value of D s (D S = 1.965 in Fig. 9d) the density of vacancies increases rather strongly at the surfaces but it remains constant in the bulk.
Of course, such rather rapid variation of the vacancy density at the boundary can be considered as a rounded remnant of the sharp transition that occurs at T = 0 when the nwo-nwd boundary is crossed (Fig. 2). Figure 10a shows the location of critical wetting transitions in the H 1 − D s plane for the same, rather high, temperature k B T /J = 1.393 as used in Fig. 9. Also, Fig. 10b shows a similar phase diagram as in (a) but obtained at k B T /J = 0.60.However, note that at this lower temperature one has a line of first order wetting transitions (see also Fig. 4b) that meets a second order one in a tricritical point close to H tri 1 /J 0.945, D tri s /J −0.2.While for sufficiently negative values of D s the transition value of |H 1c /J | is indeed almost independent of D s (similar to the groundstate phase diagram, cf.Fig. 2, see also Fig. 4b), this value at a high temperature is much lower (e.g.|H 1c |/J ≈ 0.37 for k B T /J = 1.393, rather than |H 1c |/J = 1.0 for T = 0), but becomes rather close to the groundstate value at lower temperatures (|H 1c |/J 0.945 for k B T /J = 0.60).Also, while the increase of |H 1c |/J at T = 0 starts discontinuously at D s /J = 1, at nonzero temperature the increase starts earlier, and gradually.
From the variation of the surface vacancy density with D s we can still roughly estimate some gradual effective "transition" into the nonwet surface disordered regime already identified as a true phase transition in the groundstate phase diagram of Fig. 2. Of course, while at T = 0 the density of vacancies at the surface layer is unity for this regime, also identified as a prewetting state in an early publication [43], this density decreases at non zero temperatures.We locate these "transition" points just by determining the maximum of the fluctuations of the surface vacancy density (σ SV , cf.Fig. 11).In fact, here one observes rather broad peaks at finite temperatures that leads to a δ−function just at T = 0. Furthermore, the location of the peaks close to D s /J 1.77(2) is only slightly shifted from the exact value at T = 0, namely D s /J = 1.70 for H 1 /J = 0.3 (see Eq. 12).By using the above described procedure, we found that in Fig. 10a the "transition" line runs almost vertically along the strip with 1.4 < D S /J < 2.2.Only on the left hand side of this "transition" region is the magnetization in the boundary layers significantly different from zero; while at the right hand side one essentially has a bulk ordered state with a surface disordered state (B-O & S-D).Summing up, the transition found at T = 0 becomes rounded and shifted at higher temperatures since the boundary row is a one-dimensional object that can no longer exhibit a sharp transition for T > 0. After inspection of Fig. 10a and b it follows that for D s < 0 there should be a line of tricritical surface fields almost independent of D s that could be located just by scanning a suitable temperature interval.As an example to illustrate that scenery, Fig. 12 shows that a plot of the dependence of H 1c on temperature, for D/J = 1.5, exhibits a line of first order wetting transitions that meets a second order line at the tricritical point H tri 1 /J 0.225, k B T tri /J 0.925.
It also is of interest to study the location of the wetting transition k B T w (H 1 )/J at a fixed value of |H 1 |/J in the plane of variables k B T /J and D s /J (Fig. 13).We see that for D/J = 1.5 (Fig. 13a) and negative values of D s the wetting transition is of first order, and T w (H 1 ) is almost independent of D s .However, near D s /J = 1 the wetting transition becomes of second order and the wetting temperature starts to rise.So, a tricritical point is Finally, it is of interest to consider the approach of the wetting transition field H 1c (T ) as the temperature approaches the bulk transition T cb (Fig. 14).In our previous work [19] this problem was already considered for the case D s = D, and it was shown that H 1c (T ) ∝ (T cb (D) − T ) 1 , where the exponent 1 = 1/2.When we vary D s at constant D, T cb (D) does not change, and hence we expect a family of curves which all merge on the abscissa in the same point, and this is in fact what we see (Fig. 14).From the universality principle we furthermore expect that also the exponent 1 is independent of D s , so it is only the amplitude prefactor in the quoted power law that can depend on D s .The inset of Fig. 14 shows that our data are indeed compatible with this expectation, although the data are not taken close enough to T cb (D) to make this point really convincing.

Conclusions
In this paper, we have studied the surface effects (in particular, wetting phenomena) in a generic model exhibiting both second order and first order transitions between ordered and disordered phases in the bulk, namely the Blume-Capel model on the square lattice.Due to the reduced symmetry of lattice sites in the boundary rows, we have assumed that the Hamiltonian contains two types of surface perturbations, namely both a surface magnetic field H 1 , and a surface crystal field D s different from the bulk.The surface field H 1 (at zero bulk magnetic field, for the phase with magnetization direction opposite to H 1 ) drives the wetting transition.It is found that for D s different from the bulk crystal field D both first and second order wetting transitions occur, and hence tricritical wetting transitions can be located, while for D s = D only critical wetting occurs (as long as D is chosen such that one has a second order transition in the bulk).While at zero temperature a further surface transition in the nonwet state occurs when D s is varied, namely the boundary row undergoes a transition from ferromagnetic order (the "nwo"-state in Fig. 2) to a disordered state ("nwd"), at finite temperatures the only remainder of this transition is a large fluctuation in the vacancy density (Fig. 11).This "rounding" of the transition is expected, of course, since the boundary row of a two-dimensional system is a onedimensional system, which cannot exhibit a sharp phase transition at nonzero temperature.Of course, it will be interesting to consider the generalization of the present model for d = 3 dimensions, where finite temperature transitions in the boundary phase no longer can be generally excluded a priori; however, such a study is left for future work.
A plausible expectation on the action of changing D s on the wetting transition is that it "renormalizes" the strength of the surface magnetic field (that has no effect on sites taken by vacancies of course) and hence it could only cause a shift of the wetting transition but not change its order.However, our results clearly show that such a naive idea about "coarsegraining" our model is not true, and a physical interpretation why the order of the wetting transition should change when vacancies in the boundary layers are suppressed is lacking.
We also note that the crucial ingredient for our study was the availability of a clearcut and useful finite size scaling framework for wetting phenomena (Fig. 1, Eqs. 3-6) as has been developed by us recently [19].However, since the critical slowing down due to very slow interfacial fluctuations still is a problem for the accuracy of our data, no attempt to study the tricritical wetting exponents could be made.In fact, developing numerical transfer-matrix type methods could be a powerful alternative for this problem.Thus, we hope that our study will stimulate corresponding complementary work.

Fig. 2
Fig.2Surface phase diagram corresponding to the groundstate of a semi-infinite Blume-Capel model in the region where the bulk is ordered (D/J < 2) in the plane of the variables reduced surface field (H 1 /J ) versus reduced surface crystal field D s /J .The surface phase transitions are shown as full straight lines (metastable continuations are indicated as broken straight lines).Three phases occur, a wet surface, a nonwet ordered surface (nwo) and a nonwet disordered surface (nwd), as explained in the text

Fig. 3
Fig. 3 Plots of the average absolute value of the magnetization |m| T (a), the magnetization square m 2 (b) and the fourth order cumulant U = 1 − m 4 /(3 m 2 2 ) (c), versus temperature.All data are for L × M lattices in the geometry of L × M with L = 12, 24, and 36, where c = L 2 /M = 9/8 is kept constant, choosing D = 0, D s /J = −19.65.Curves connecting the simulation data (shown by symbols) are only meant as guides to the eye

Fig. 4 a
Fig.4a Temperature variation of the absolute value of the magnetization, for the case H 1 /J = 0.55, D/J = 1.50, D s /J = −19.65,and lattices of size 24 × 512 and 36 × 1152, respectively.Initial conditions were either putting all spins pointing up (data points denoted by circles) or putting spins pointing up in rows from 1 to L/2 − 1, pointing down in rows from L/2 + 1 to L, but S i = 0 if i is in row L/2 (triangles and squares).b Dependence of tha absolute value of the magnetization on the surface field (H 1 /J ), obtained for k B T /J = 0.60 by using lattices of size 18 × 288 and 24 × 512, and two choices of the surface crystal field, as indicated

Fig. 5
Fig. 5 Typical snapshot configurations of the Blume-Capel model in d = 2 dimensions for conditions near a first order wetting transition obtained for a system size of L = 24 and M = 512.The surface field H 1 /J = 0.55 acts always on the right boundary of the strip and while −|H 1 |/J = −0.55acts at the left boundary of the strip.D/J = 1.5 and D s /J = −19.65 are chosen.Spins pointing up are shown in blue (dark-grey when the figure is displayed in black/white tones), spins pointing down in white, and spins S i = 0 in red color (gray when the figure is displayed in black/white tones), respectively.Three temperatures are shown, k B T /J = 0.579 (left), 0.591 (middle), and 0.603 (right).Note that vacancies occur in both "spin up" and "spin down" phases in the bulk, and are strongly enriched at the fluctuating, meandering interface, while they are completely expelled from the surface (Color figure online)

Fig. 6
Fig. 6 Wetting phase transitions in the plane of variables D/J and k B T /J , for the particular choice H 1 /J = 0.55 and two choices of the surface crystal field, namely D s = D (open triangles) and D s /J = −19.65 (full triangles, critical wetting, and diamonds, first order wetting).First order wetting points are evaluated by using lattices of size L = 24, and M = 512.Second order wetting points are evaluated by using the intersection method described in Sect.2.3 with lattices of size L = 12, 24, and 36, by keeping c = L 2 /M = 9/8 constant.The full circles (connected by broken lines) show the bulk second order transition of the Blume-Capel model, while the squares show the bulk first order part of this transition line.Here, the bulk tricritical point at k B T tri b /J 0.609, D tri /J 1.965 is highlighted by a diamond.Note that the wetting transition line for D s /J = −19.65 exhibits a wetting tricritical point (highlighted by an open square)

Fig. 7
Fig. 7 Profiles of the magnetization (a) and density of vacancies (b) for the case of critical wetting with |H 1 |/J = 0.55, D s /J = −19.65,L = 24, M = 512, and the three choices of the bulk crystal field and the temperature, namely D/J = 0.25 (k B T /J = 1.14), 0.50 (k B T /J = 1.05), and 0.75 (k B T /J = 0.96), respectively.c and d same as a and b, but for the case of first order wetting with L = 36, M = 1152, and D/J = 1.5 (k B T /J = 0.585), 1.25 (k B T /J = 0.731), and 1.0 (k B T /J = 0.859)

Fig. 8
Fig. 8 Magnetization profiles (a) and profiles of the vacancy density (b) for a strip choosing L = 36, M = 576, D = 0, D s /J = 19.65,H 1 /J = 0.55 and three choices of the temperature k B T /J , as indicated.c The absolute value of the average magnetization |m| T is plotted versus the reduced temperature k B T /J for three system sizes, L = 18, 24, and 36, choosing a constant aspect ratio L/M = 1/16.d Presents a standard finite size scaling [21,35-37] plot of the data on the left side, i.e. |m| T L β/ν is plotted versus (k B |T − T cb |/J )L 1/ν , using d = 2 Ising model critical exponents (β = 1/8, ν = 1 [26]).The only adjusted value is the critical temperature k B T cb /J = 1.695(5) corresponding to the order disorder transition of the Blume Capel model in the bulk for D = 0.A straight line with exponent β = 1/8 is also included to show the consistency of the data

Fig. 9 a
Fig. 9 a and b Magnetization profiles across a strip of width L = 24 rows (M = 512) for the choice D = 0, k B T /J = 1.393, |H 1 |/J = 0.2, and several choices of D s /J , as indicated.c and d Same as (a-b), but for the profiles of the vacancy density

Fig. 10 a
Fig. 10 a Phase diagram showing the localization of the critical wetting transition curve in the plane of variables |H 1 |/J and D s /J , for D = 0 along an isotherm k B T /J = 1.393 (full dots).The squares indicate the crossover of the boundary layer from nonzero boundary magnetization to (essentially) zero magnetization, due to large vacancy density at the boundary.b as in a but for the choice k B T /J = 0.6.Here full triangles correspond to first order wetting transitions that meet critical wetting transitions (full circles) at a tricritical point (TP) indicated by means of an arrow.First order wetting points are evaluated by using lattices of size L = 24, and M = 512.Second order wetting points are evaluated by using the intersection method described in Sect.2.3 with lattices of size L = 12, 24, and 36, by keeping c = L 2 /M = 9/8 constant.Also, the effective transition points are determined as described in Fig. 11 by using lattices of size L = 24, and M = 512.The acronym B-O & S-D indicates a bulk ordered state with a surface disordered state

Fig. 11 Fig. 12
Fig. 11 Log-linear plots of the fluctuations of the density of vacancies at the surfaces n = 1 and n = L (σ SV ) versus D s /J , with L = 24, M = 512 for the choice D = 0, |H 1 |/J = 0.3, and two temperatures, as indicated.Here, larger fluctuations are observed for each temperature, at the wall where the interface is bound (recall that all data is taken within the non wet phase).The vertical full line placed at D s /J = 1.7 shows the δ−function corresponding to the groundstate

Fig. 13
Fig. 13 Phase diagram showing the location of the wetting transition in the plane of variables k B T /J and D s /J for fixed surface field; |H 1 |/J = 0.55, for D/J = 1.5 (a) and D/J = 0 (b).Triangles show the location of the order disorder transition in the bulk, which is independent of D s , of course (k B T cb /J = 1.153(D/J = 1.5) and k B T cb /J = 1.695(D = 0), respectively).Full squares and circles identify first-and second-order wetting transitions.In (a) the location of the tricritical point (TP) is shown by means of an arrow, while in (b) all wetting transitions are of 2nd order

Fig. 14
Fig. 14 Critical field H 1c (T )/J where the wetting transition in the Blume-Capel model occurs, for D = 0 and three choices of D s as indicated, plotted as a function of temperature k B T /J .Inset shows the same data replotted versus (k B T cb /J − k B T /J ) 1/2 [19,42]r at H 1 = H 1c , if H 1 is varied at constant temperature T .Of course, H 1c (T ) is nothing but the inverse function of T w (H 1 ) in the (T, H 1 ) plane.We have shown[19,42]that the average absolute value of the magnetization |m| of the strip, as well as m 2 and the cumulant U = 1 − m 4 /(3 m 2 2 ) exhibit a simple finite size scaling behavior, when M is varied at constant c = L 2 /M, for large enough L, given by