Critical Casimir effect in films for generic non-symmetry-breaking boundary conditions

Systems described by an O(n) symmetrical $\phi^4$ Hamiltonian are considered in a $d$-dimensional film geometry at their bulk critical points. A detailed renormalization-group (RG) study of the critical Casimir forces induced between the film's boundary planes by thermal fluctuations is presented for the case where the O(n) symmetry remains unbroken by the surfaces. The boundary planes are assumed to cause short-ranged disturbances of the interactions that can be modelled by standard surface contributions $\propto \bm{\phi}^2$ corresponding to subcritical or critical enhancement of the surface interactions. This translates into mesoscopic boundary conditions of the generic symmetry-preserving Robin type $\partial_n\bm{\phi}=\mathring{c}_j\bm{\phi}$. RG-improved perturbation theory and Abel-Plana techniques are used to compute the $L$-dependent part $f_{\mathrm{res}}$ of the reduced excess free energy per film area $A\to\infty $ to two-loop order. When $d<4$, it takes the scaling form $f_{\mathrm{res}}\approx D(c_1L^{\Phi/\nu},c_2L^{\Phi/\nu})/L^{d-1}$ as $L\to\infty$, where $c_i$ are scaling fields associated with the surface-enhancement variables $\mathring{c}_i$, while $\Phi$ is a standard surface crossover exponent. The scaling function $D(\mathsf{c}_1,\mathsf{c}_2)$ and its analogue $\mathcal{D}(\mathsf{c}_1,\mathsf{c}_2)$ for the Casimir force are determined via expansion in $\epsilon=4-d$ and extrapolated to $d=3$ dimensions. In the special case $\mathsf{c}_1=\mathsf{c}_2=0$, the expansion becomes fractional. Consistency with the known fractional expansions of D(0,0) and $\mathcal{D}(0,0)$ to order $\epsilon^{3/2}$ is achieved by appropriate reorganisation of RG-improved perturbation theory. For appropriate choices of $c_1$ and $c_2$, the Casimir forces can have either sign. Furthermore, crossovers from attraction to repulsion and vice versa may occur as $L$ increases.


Introduction
When media exhibiting low-energy fluctuations are confined by walls or interfaces, or macroscopic objects are immersed in them, their fluctuation spectrum changes [1,2,3]. This may induce long-range effective forces between such objects and boundaries. A much studied example of such fluctuation-induced forces are the Casimir forces between a pair of parallel and grounded metallic plates produced by their influence on the quantum vacuum fluctuations of the electromagnetic field [4]. These quantum electrodynamics (QED) Casimir forces are known to be very weak unless the separation of the bodies between which they act becomes very small. The latter condition is met in micro-electromechanical devices (MEMS). As has recently become clear, Casimir forces of this kind must be taken into account in the design of such systems. Since they tend to be attractive for simple geometries, they may impair the functioning of MEMS, causing stiction [5,6,7,8]. Owing to the interest in -and technological importance of -such small-scale systems, this insight has triggered considerable new activity in the study of QED Casimir forces. Crucial issues are the knowledge and control of their strength, sign, and geometry dependence [9].
Another example of fluctuation-induced forces that has attracted a great deal of attention recently are the so-called thermodynamic Casimir force caused by thermal fluctuations near critical points [2,10,11]. Predicted decades ago [12], they were verified experimentally first in an indirect way through the thinning of 4 He wetting layers near the lambda transition [13] and their effects on wetting films of binary liquid mixtures [14,15]. Subsequently they were measured directly in binary fluid mixtures near their critical points of demixing [16,17].
In this paper we shall be concerned with thermodynamic Casimir forces in systems with a film geometry bounded by two planar free surfaces at a finite distance L. We will consider systems that can be modelled by an n-component |φ| 4 Hamiltonian on sufficiently long length scales and study the Casimir forces directly at the bulk critical point. The internal O(n) symmetry (or Z 2 symmetry if n = 1) is presumed to be neither explicitly broken by surface contributions to the Hamiltonian nor spontaneously for all temperatures T ≥ T c,∞ , the bulk critical temperature. In other words, no ordered phase can occur for finite L when T ≥ T c,∞ . Choosing a generic kind of such non-symmetrybreaking conditions, we will investigate the critical Casimir forces and demonstrate two important features of them. First, they can have either sign, depending on properties of the two surfaces -an expected result since attractive and repulsive critical Casimir forces were found a long time ago for certain boundary conditions that ought to apply on long length scales [18,19,20]. Second, for appropriate surface properties, crossovers from attraction to repulsion of the Casimir forces and vice versa can occur as a function of film thickness L. A brief account of parts of our work was given in [21]. The purpose of the present paper is to provide details of the calculations and further results.
One motivation for our work is the obvious potential for cross-fertilisation between the fields of thermodynamic and QED Casimir effects. Since both types of effects share a number of characteristic features, mutual benefits may be expected. On the other hand, one must be aware of certain essential differences. Common to both types is that they exhibit universal properties and depend on gross properties of the fluctuating media and confining objects, as well as on their shapes and geometry. Two important differences are: (i) When studying QED Casimir forces, one frequently gets away with the analysis of effective Gaussian theories in which the coupling between the electromagnetic field and matter is accounted for by proper choices of boundary conditions. This usually holds even in the case of polarisable and magnetisable media where material properties also enter via dielectric and permeability functions. By contrast, adequate treatment of critical Casimir forces normally involve the study of non-Gaussian theories such as φ 4 theories or corresponding lattice models (e.g., Ising models) in bounded geometries.
(ii) In studies of QED Casimir forces such as Casimir's original work [4], the electromagnetic fields are taken to have zero averages. Hence the Casimir forces are entirely due to fluctuations -if there were no fluctuations, there would be no Casimir force. However, in the case of thermodynamic Casimir forces, the order-parameter densities φ(x) have nonzero averages in ordered phases, where the order may either result from spontaneous symmetry breaking or else be imposed by symmetry-breaking bulk or boundary fields. If the medium undergoes a phase transition from a disordered high-temperature to an ordered low-temperature phase, then the order-parameter profile φ(x) does not vanish in the latter. Such a nontrivial profile contributes to the size-dependent part of the free energy and hence, to the Casimir force. Consequently, a Casimir force will be obtained already at the level of Landau theory -i.e., in the absence of fluctuationswhenever it yields non-vanishing order-parameter profiles. Beyond Landau theory, the thermodynamic Casimir force will therefore consist of a part coming from a non-fluctuating background and a superimposed fluctuation-induced contribution.
The last statement applies in particular to confined binary fluid mixtures [17,22,23]. These are known to generically involve symmetry-breaking boundary fields. Therefore, the thermodynamic Casimir forces they yield will normally have contributions from non-fluctuating backgrounds on either side of the order-disorder transitions. Since even at the bulk critical temperature, φ(x) is not expected to vanish, this holds there as well.
It is well known that the boundary fields needed to describe binary fluid mixtures in contact with walls can have either sign, depending on which one of the two components gets preferentially adsorbed locally at the wall [10,12,22,23,24]. Furthermore, there is clear evidence based on theoretical and experimental work as well as on Monte Carlo simulations [10,14,15,16,25,26,27,28,29,30] that the critical Casimir force in binary fluid mixtures confined between two planar walls can be attractive or repulsive. Corresponding crossovers must occur and were investigated at the level of mean-field theory [31].
Casimir forces in binary fluid mixtures are particularly well suited for direct experimental measurements. However, from our above remarks it is clear that their analogy with QED Casimir forces is limited: owing to the generic presence of a nonvanishing order parameter profile at, above and below the critical temperature, only a part of them is fluctuation induced. Since we focus here on the case in which the symmetry φ → −φ remains unbroken at T c,∞ , we are dealing with critical Casimir forces that are entirely fluctuation induced and hence are more akin to QED Casimir forces, albeit due to thermal, rather than quantum, fluctuations.
The rest of this paper is organised as follows. In section 2.1 the model is specified and the boundary conditions are recapitulated. Next, the free propagator is determined in section 2.2 for general values of the surface interaction constants, its eigenfunction representation is given and the transcendental equations derived whose solutions yield its eigenvalues. In section 2.3, many-point cumulant functions are introduced, their renormalization recapitulated and their RG equations presented. These considerations are extended to the free energy and Casimir force in section 3. In the remainder of this section details of our approach are explained and our analytical results for the scaling functions of the L-dependent part of the free-energy density and the Casimir force are given. In section 3.5 it is shown how RG-improved perturbation theory can be reorganised to achieve consistency with the fractional series expansions in one encounters in the special case of two critically enhanced surface planes. Section 4 discusses our results, deals with the issue of how to extrapolate them to d = 3 dimensions and presents such extrapolation results for the scaling functions of the L-dependent part of the excess free-energy density and the Casimir force. Section 5 contains a brief summary and concluding remarks. Finally, there are 5 appendices describing various calculational details.

Model and boundary conditions
We consider a continuum model for a real-valued n-component order-parameter field We write the position vector as x = (x 1 , . . . , x d ) = (y, z), where y = (y 1 , . . . , y d−1 ) denotes the d − 1 Cartesian coordinates y j = x j along the slab and z = x d the remaining one across it. The thickness of the slab, L, is taken to be finite. We choose periodic boundary conditions along all y-directions so that the boundary ∂V of V consists of a pair of (d−1)-dimensional planes B 1 and B 2 at z = 0 and z = L. Assuming the absence of long-range interactions, we can choose the Hamiltonian to have the local form where L b and L j=1,2 are local bulk and surface densities depending on φ(x) and its derivatives.
In accordance with our considerations in the introduction we choose these densities to be O(n) symmetric. For the bulk density we make the standard choice where (∇φ) 2 is the usual short-hand for α (∇φ α ) 2 . It has recently been emphasised that lattice systems lacking cubic symmetry, such as systems involving monoclinic or triclinic lattices, lead to φ 4 bulk densities with a gradient-square term of the generalised form 1 2 g ij (∂ i φ) · ∂ j φ where ∂ i stands for the derivative ∂/∂x i and g ij is a non-diagonal, position-independent matrix [32]. We refrain from working with such generalised models because the matrix g ij is nothing but a constant metric [33]. The geometric effects it has can be absorbed by a proper choice of variables which of course enters the way lattice quantities are related to those of the standard bulk φ 4 model with action density (2.2). For details the reader might want to consult reference [33, p. 14-17].
The surface densities are given by The interaction constantsc 1 andc 2 of the two boundary planes are allowed to differ. They are known to reflect the weakening or enhancement of the local pair interactions at B 1 and B 2 , respectively [22,23]. In semi-infinite systems bounded by a plane B 1 , a threshold valuec sp of the enhancement variablec 1 exists such that a phase with longrange surface order appears (in sufficiently high dimensions) in a temperature regime above T c,∞ when the enhancement −δc 1 ≡c sp −c 1 is positive (is "supercritical"). The transitions that occur at T c,∞ in such semi-infinite systems are called ordinary, special and extraordinary, depending on whether the enhancement δc 1 is subcritical (δc 1 > 0), critical (δc 1 = 0) or supercritical (δc 1 < 0). Analogous statements hold for semi-infinite systems bounded by the plane B 2 with surface enhancementc 2 .
In order to rule out spontaneous symmetry breakdown at T = T c,∞ for large L, we require that both enhancement variables δc j ≡c j −c sp satisfy the condition δc j ≥ 0, i.e. are non-supercritical.
From the boundary contribution to the classical equation of motion we get the boundary conditions of Landau theory where ∂ n is the inner normal derivative on ∂V = B 1 ∪ B 2 . Beyond Landau theory these boundary conditions still hold in an operator sense (inside of averages) [22,23,34,35,36].

Free propagator
A quantity of central importance for the renormalization-group (RG) improved perturbation theory we are going to use below is the free propagator G L atτ = 0. It is given by the operator inverse G L = (− ) −1 subject to the boundary conditions (2.4), where = ∇ 2 is the Laplacian. Let |m be the complete set of eigenfunctions of the operator −∂ 2 z ≡ −∂ 2 /∂z 2 on the interval [0, L] satisfying the boundary conditions That is, we have and the eigenstates fulfil the orthogonality and completeness relations m|m = δ mm (2.8) and ∞ m=1 z|m m|z = δ(z − z ), (2.9) where the label m = 1, 2, . . . , ∞ enumerates the eigenvalues k 2 m by size, starting with the smallest one.
To exploit the translation invariance of the system along the direction parallel to the boundary planes B j , we Fourier transform with respect to the difference of y coordinates, writing where we have introduced the notation We thus arrive at the spectral representation (2.12) Since our assumption that both enhancement variablesc j are non-supercritical implies that the associated renormalized variables c j (whose definition will be recalled in the following section 2.3) are nonnegative, we shall need the eigenvalues k 2 m and eigenfunctions |m only forc j ≥ 0. In this case, there are infinitely many eigenvalues k 2 m for finite L, all of which are nonnegative and non-degenerate (see reference [37, Appendix A] and below). A zero eigenvalue k 2 1 = 0 occurs for arbitrary L ∈ (0, ∞) only whenc 1 =c 2 = 0.
Let us briefly summarise the relevant properties of the eigensystem {|m , k 2 m } we shall need in our subsequent analysis. On dimensional grounds the eigenfunctions must have the form z|m = L −1/2 υ m (z/L|c 1 L,c 2 L). (2.13) It is convenient to introduce the dimensionless variables ζ ≡ z/L, κ m ≡ Lk m , C j ≡c j L. (2.14) Setting L = 1 in equation (2.13) one concludes that the functions υ m can be written as The analogue of the boundary condition (2.5) implies that the phase shift ϑ m is given by We can choose it such that 0 ≤ ϑ m ≤ π/2, which leads to sin ϑ m = 1 The limiting values ϑ m = 0 and ϑ m = π/2 are obtained whenc 1 → 0 andc 1 → ∞, respectively.
Using these results along with equation (2.15), the normalisation constant λ m can be computed in a straightforward fashion. One obtains To determine the spectrum {κ 2 m }, we use the analogue of the boundary condition (2.6), (2.20) This yields the transcendental equation for κ m , with Equations (2.5)-(2.7) specify a regular Sturm-Liouville problem for which the following mathematical properties are known (see e.g. references [38], [39,Kap. IX], [40,Kap. IV.3] and [41, chapters 8, 9]): (i) The eigenvalues κ 2 m (C 1 , C 2 ) are real, nondegenerate, countable and accumulate only at ∞. (ii) They can be ordered such that κ 2 m < κ 2 m for m < m . There is a smallest eigenvalue κ 2 1 but no largest one, i.e. κ m → ∞ as m → ∞. (iii) The eigenfunctions υ m corresponding to different eigenvalues are orthogonal with respect to the standard L 2 scalar product in L 2 (0, 1). ‡ (iv) The normalised eigenfunctions are a complete orthonormal set (basis) in the Hilbert space L 2 (0, 1). (v) When C 1 and C 2 > 0, we can use the Robin boundary conditions (2.5) and (2.6) in conjunction with the fact that the eigenfunctions vanish only for Dirichlet boundary conditions at the boundary planes B j to conclude that υ(dυ/dζ)| ζ=1 ζ=0 = −C 2 υ 2 | ζ=1 − C 1 υ 2 | ζ=0 < 0. Fulfilment of this condition guarantees (by theorem 7 of reference [40, p. 234]) that κ 2 1 and hence all κ 2 m are strictly positive. The non-degeneracy of the eigenvalues κ 2 m can be verified explicitly by showing that the positive zeros κ m of R C 1 ,C 2 (κ) are simple. To this end, one can compute the derivative R C 1 ,C 2 (κ) ≡ ∂ κ R C 1 ,C 2 (κ) at κ = κ m and express the trigonometric functions cos κ m and sin κ m in terms of κ m , C 1 and C 2 using equations (2.21) and (2.17) along with the analogue of the latter implied by the boundary condition at ζ = z/L = 1. This gives , (2.23) which is nonzero for all C 1 ≥ 0, C 2 ≥ 0 and κ m > 0.  [20]. However, for general values of the C j , the zeros κ 2 m of the transcendental equation (2.21) cannot be determined in closed analytical form. This is a familiar difficulty, which makes the evaluation of the single and double mode sums m one encounters in the calculation of one and two-loop Feynman graphs of the free energy a nontrivial problem. We shall turn to this issue in section 3. Note also that in some of the above formulae we have tacitly assumed that C 1 and C 2 are both positive. However, these equations remain valid, firstly, for all m when either one of the C j vanishes while the other remains positive, and secondly, for all m > 1 when C 1 = C 2 = 0, because the respective κ m then all approach positive values in the corresponding limits. Hence one can simply set C 1 or C 2 (or both) to zero in the respective equations. The case of m = 1 with C 1 , C 2 → 0 is special in that κ 1 → 0. It is an easy matter to check that the correct limiting eigenfunction υ(ζ|0, 0) = 1 results, independent of the order in which the limits C 1 → 0 and C 2 → 0 are taken.

Many-point cumulants, their renormalization and RG equations
We proceed by providing some of the necessary background on the renormalization of the theory defined above. § Let us begin by recalling the bulk and boundary counter-terms required to absorb the ultraviolet (UV) singularities of many-point cumulant functions. To this end we consider the (N +M 1 +M 2 )-point cumulant functions (2.25) involving N fields φ α j (x j ) located at points in the interior of V and M 1 +M 2 fields φ β k on the boundary planes B 1 and B 2 . As indicated, we take the first M 1 boundary points to lie on B 1 and the remaining ones on B 2 .
Equations (2.12)-(2.17) imply that the p transform of the free propagator for general nonnegative values ofc 1 andc 2 can be written aŝ . (2.26) Its explicit form is known [33,42]. It is given by the expression into which the right-hand side of equation (B5) of reference [33] turns upon replacement of itsκ ω by i.e. (2.28) Both expressions (2.26) and (2.28) are exceedingly difficult to work with. To understand the nature of the UV singularities of the theory, it is better to treat the boundary terms L j as interactions and work with the free propagator forc 1 =c 2 = 0. The latter is nothing but the free propagatorĜ NN is the free bulk propagator forτ = 0 in the pz-representation. As expounded in reference [22], the j = 0 contribution from the first term in square brackets in equation (2.29) is the origin of the familiar primitive bulk UV singularities. The j = 0 and j = −1 contributions from the second term in square brackets are singular at coinciding points on B 1 and B 2 , respectively. The additional UV singularities they produce can be absorbed by counter-terms with support on these boundary planes, i.e. the surface counter-terms required for the corresponding semi-infinite systems [22,35,36]. All other contributions in the square brackets do not diverge at coinciding points and hence do not give rise to additional UV singularities. The upshot is that the usual bulk reparametrizations of the surface enhancement variablesc j and the boundary operators φ B 1 (y) ≡ φ(z, 0) and φ B 2 (y) ≡ φ(z, L) suffice to absorb the UV singularities of the functions G (N ;M 1 ,M 2 ) .
Following reference [43], we choose the factor that is absorbed in the renormalized interaction constant as where = 4 − d and γ E is the Euler-Mascheroni constant. This choice ensures that the bulk renormalization factors Z φ , Z τ and Z u as well as the surface renormalization factors Z 1 and Z c , when determined by minimal subtraction of poles in , up to second order in u reduce to the two-loop results given in references [22], [34], [35], [36], [44] and used in Krech and Dietrich's work [19,20]. The quantityτ c,∞ is the critical bulk value ofτ . Further,c sp is the critical enhancement value associated with the special surface transition. In a theory whose UV singularities are regularised by means of a large-momentum cut-off Λ, these quantities diverge ∼ Λ 2 and ∼ Λ, respectively. In our calculations below we shall utilise dimensional regularisation. Note also that the surface renormalization factors Z 1 and Z c depend exclusively on u and but not on c j (nor on L) when fixed by the requirement that the UV poles in be minimally subtracted. In our calculations of the L-dependent part of the free energy at the bulk critical point and the critical Casimir force to be described in section 3 we shall need Z c merely to first order in u. We therefore quote the result for convenience. Upon introducing the renormalized cumulants we can now derive RG equations for these functions by varying µ at fixed bare parameters u,τ ,c 1 ,c 2 , and L. Let us define the exponent functions (where ∂ µ | 0 indicates a µ-derivative at fixed bare parametersů,τ ,c 1 ,c 2 , L) and the operator Then the RG equations can be written as They are completely analogous to those for the renormalized cumulants of the respective semi-infinite geometries (L = ∞ with either M 2 = 0 or M 1 = 0). We have given them here for general τ , even though we will restrict ourselves to the critical case τ = 0 in our calculations in section 3.

Definitions
We next turn to the free energy and the Casimir force. Introducing the partition function Z and the total free energy F via we define the reduced free-energy density per hyper-surface area A, This quantity can be decomposed as where f b and f s are the bulk free-energy density per d-dimensional volume LA → ∞ and the excess free-energy density per hyper-surface area A → ∞, both of which are defined by appropriate thermodynamic limits (see e.g. references [22] and [45]). The former depends only on the bulk interaction constants, the latter additionally on the boundary interaction constants. It is a sum of the reduced surface excess free-energy densities f s,j of the corresponding semi-infinite systems bounded by either B 1 or B 2 : The remaining term in the decomposition (3.3), the reduced residual free-energy density f res , contains the full L-dependence of f − Lf b . The effective ("Casimir") force per hyper-surface area to which it gives rise is where the ellipsis stands for all other (bulk and surface) variables on which it depends.

Renormalization of the residual free energy and Casimir force
Inspection of the perturbation series of the reduced free energy F/k B T reveals that the bulk reparametrizations (2.31) and surface reparametrizations (2.32) are not sufficient to absorb all UV singularities. This is because both the bulk free-energy density f b as well as the surface free-energy densities f s,j have additional primitive UV singularities [22]. To cure these, additive bulk and surface counter-terms are required. Since these counter-terms must merely absorb primitive UV singularities of f b and f s,j , respectively, they can be chosen independent of L. One convenient way to fix them is by subtracting from f b the Taylor expansion in δτ , and from f s,j that in δτ and δc j about nonvanishing reference values to the appropriate orders (see e.g. [ Hence this quantity satisfies a homogeneous RG equation, namely The latter can be solved in a standard fashion by means of characteristics. Let us define the running variablesū( ) andc j ( ) as solutions to the initial value problems Let u * be the infrared-stable zero of the beta function β u for d = 4 − < 4. We shall need its expansion below only to O( ). Solving the flow equations (3.9) forc j yields where Φ is the surface crossover exponent whose expansion is known to O( 2 ) [35,36]. Since we shall need the expansion of the RG eigenexponent y c to first order in , we recall the result In the large-length-scale limit → 0, the running variablesū andc j behave as respectively, where E * c (u) is a non-universal amplitude. We now set τ = 0, solve the RG equation (3.7), choose the flow parameter as µL = 1 and use the above limiting expressions forū andc j . This gives the asymptotic behaviour with The analogous scaling form of the reduced critical Casimir force per hyper-surface area follows by differentiation with respect to L. It reads To appreciate these results, recall that the large-L form of the residual free-energy density at the bulk critical point T = T c,∞ , for given asymptotic large-scale boundary conditions BC, is conventionally written as where ∆ BC is an L-independent universal, but BC dependent, number, called "Casimir amplitude". According to our result (3.15), a scaling function D appears for general values of c 1 and c 2 in place of the Casimir amplitude -that is, the Casimir amplitude becomes scale-dependent.
On the other hand, various Casimir amplitudes ∆ BC can be recovered from the knowledge of D(c 1 , c 2 ). Since we assumed both surface enhancement variables c j to be non-supercritical, and also ruled out the breaking of the O(n) symmetry via boundary terms, the case of Robin boundary conditions we consider includes the four cases The first two of these amplitudes, ∆ (O,O) and ∆ (O,sp) , are known to have expansions in integer powers of . The leading two terms of these series expansions were determined in reference [20]. The corresponding results can be written as and and By contrast, the Casimir amplitude ∆ (sp,sp) does not have a power-series expansion in . As shown in references [43] and [46], the presence of the k 1 = 0 mode leads to a breakdown of the expansion of ∆ (sp,sp) and produces additional half-integer powers j/2 with j = 3, 5, . . ., which are modulated by powers of ln when j ≥ 5. The expansion is known to order 3/2 ; it reads [43,46] The -expansion results for the scaling function D(c 1 , c 2 ) we are going to derive below must be consistent with the series expansions (3.21)-(3.25) and equation (3.20). This requirement provides nontrivial checks for these results.

Perturbation theory
We next turn to the loop expansion of the reduced free-energy density and its computation forτ = 0 to two-loop order. The zero-loop contribution f [0] vanishes in the disordered phase. The next two terms can be written as and for generalτ ≥ 0. Here we have introduced the functions For the parameter values for which they are needed in equation (3.29), a lengthy but straightforward calculation yields with the polynomials (3.32) The p-integrals in the above equations (3.28) and (3.29) can be handled in a standard way using dimensional regularisation. The main difficulty one then is faced with is the calculation of the resulting single and double sums over the not explicitly known eigenvalues κ 2 m . We have done this by means of complex integration, modifying the Abel-Plana techniques described in reference [37] for our purposes. ¶ The technical details are given in Appendix A-Appendix B. Here we just present our results.
As before, we use the notation Q [l] to specify the l-loop term of a quantity Q. We choose the additive constant in the free-energy density such that the dimensionally regularised bulk free-energy density f b vanishes forτ = 0. Hence Our results for the surface free-energy densities read and f [2] s,j (τ = 0,ů,c j ) = where we have introduced the familiar quantity In order to write the one-and two-loop residual free-energy densities in a compact fashion, it is helpful to define the functions (3.39) ¶ As we show in Appendix A.2, it can also be derived from the result (2.26) for the free propagator. and Z (d) In terms of these, our results can be written as and (3.42)

Renormalized residual free-energy density and scaling functions
The functions X (d,1) c 1 L,c 2 L have simple poles at d = 4, caused by the behaviour of the respective pre-factor of the integral in equation (3.38). These UV singularities imply that the two-loop term (3.42) is not regular at = 0. It has a simple pole originating from the terms proportional to X (d,1) . It is an easy matter to check that this pole gets cancelled by the contribution ∝ u/ one obtains from f [1] res (L; 0,c 1 ,c 2 ) upon expressing the bare variablesc j in terms of their renormalized analogues c j via equations (2.32) and (2.34). Substitution of the one-and two-loop terms (3.41) and (3.42) into equation (3.6) therefore yields indeed a UV-finite O(u) result for the renormalized residual free-energy density f res,R (L; 0, u, c 1 , c 2 , µ), namely and J (σ) .

(3.46)
To obtain the expansion of the scaling function D, we set u = u * . Using (3.47) one sees that the result is consistent with the predicted scaling form (3.14) and yields the expansion For the special values c j = 0 and ∞, the functions D 0 and D 1 can be analytically computed in a straightforward manner. One finds where a 1 (n) is the coefficient defined in equation (3.24). These results confirm the validity of the relations (3.20) to first order in . However, the present form (3.48) of our result does not yield the contribution ∼ 3/2 to D(0, 0) = ∆ (O,sp) /n. The reason is that we assumed that c 1 + c 2 > 0 so that the lowest eigenvalue k 2 there is a zero mode in the free propagator. As already mentioned, this causes a breakdown of the conventional RG-improved perturbation theory at T c,∞ [43,46,47]. In the following we will improve our results in such a manner that they fully comply with all of the small-expansions of Casimir amplitudes given in equations (3.21)-(3.25), including the fractional expansion (3.25) of ∆ (sp,sp) to O( 3/2 ).

3.5.
Modified RG-improved perturbation theory 3.5.1. Formulation and results It is known from references [43], [46] and [47] how one can cope with the mentioned zero-mode problem one encounters at τ = 0 when c 1 = c 2 = 0: one can reorganise RG-improved perturbation theory such that it becomes well-defined at τ = 0. This suggests an obvious strategy to ensure consistency with the small-expansions of ∆ (sp,sp) to 3/2 . We should reorganise RG-improved perturbation theory for general c 1 > 0 and c 2 > 0 in a way that reduces to the one used in previous work [43,46] for the case c 1 = c 2 = 0.
We now setτ = 0 and define an effective (d − 1)-dimensional field theory with Hamiltonian H eff [ϕ] by integrating out ψ: (3.55) To determine H eff [ϕ], we use perturbation theory. The contribution to the two-point vertex caused by the coupling between ϕ and ψ, to first order inů, originates from the graph . (3.57) It represents a local interaction corresponding to a shift δτ L ofτ . The Hamiltonian To obtain the free-energy density f (L; 0,ů,c 1 ,c 2 ) in the present modified perturbation scheme, we must add to f ψ (L; 0,ů,c 1 ,c 2 ) the contribution associated with the ϕ-field. We denote it as f ϕ and define it via ϕ (L; 0,ů,c 1 ,c 2 ) + f [2] ϕ (L; 0,ů,c 1 ,c 2 ) + . . . .

(3.61)
Here f [1] ϕ and f [2] ϕ are the one-and two-loop terms of a (d − 1)-dimensional φ 4 theory with quadratic and ϕ 4 interaction constants δτ L +k 2 1 andů∆ 1,1,1,1 /L, respectively. Their explicit expressions may be inferred from equation (4.28) of reference [43]. They read f [1] ϕ (L; 0,ů,c 1 , Building on the reorganisation of perturbation theory just described, we now wish to compute the residual free-energy density f res (L; 0,ů,c 1 ,c 2 ) and express it in terms of renormalized variables to obtain improved results for its scaling function D(c 1 , c 2 ) and that of the Casimir force. Clearly, from these results all Taylor expansions in reported in sections 3.  [46] and [43] when c 1 = c 2 = τ = 0.
What makes this approach technically difficult to implement is that the eigenvalues k 2 m , whose behaviours for smallc 1 +c 2 we gave in equation (2.24), are not explicitly known for general values ofc 1 andc 2 . Hence we shall have to resort again to numerical means.
Let us start by taking a look at the effective two-point vertex whose one-loop approximation appears in the expressions (3.62) and (3.63) for f [1] ϕ and f [2] ϕ . Its renormalized counterpart is given by γ (2) . Thus its Fourier p-transform γ where c j denote the scaling variables defined in equation (3.16). Further, E * h is a non-universal amplitude whose representation as a trajectory integral can be found in equation (3.85c) of reference [22] but will not be needed in the remainder. Since this amplitude drops out of the Casimir force, we need not keep track of it and may therefore set it to 1 henceforth.
In Appendix C we determineγ R to one-loop order, show that it is UV finite and compute the scaling function Ω(c 1 , c 2 ) to O( ). The result is where κ 1 and λ 1 represent the eigenvalue κ 1 (c 1 , c 2 ) and the normalization factor λ 1 (c 1 , c 2 ), respectively. The extrapolation one obtains for Ω(c 1 , c 2 ; d = 3, n = 1) by evaluating the above O( ) result at = 1 is depicted in figure 2. Since the eigenvalues κ 1 (c 1 , c 2 ) are analytically known for all combinations of (c 1 , c 2 ) with c j = 0, ∞, j = 1, 2, the expansions of the corresponding limiting values of Ω can be determined exactly. The and Ω(c 1 , c 2 ) = c 1 + c 2 + O(c 2 1 , c 2 2 , c 1 c 1 ) + n + 2 n + 8 (3.70) The first term on the right-hand side reflects the small-c j behaviour κ 2 1 (c 1 , c 2 ) ≈ c 1 + c 2 implied by equation (2.24). The remaining terms originate from the expansion of the shift δτ to linear order, whose zeroth-order term is according to equation (14) of reference [46] (where it was denoted as δτ Just as in our calculation of γ R , the pole term gets cancelled by the counter-term provided by the contribution linear inc 1 +c 2 = µZ c (c 1 + c 2 ) to k 2 1 . If we substitute equation (3.71) into k 2 1 + δτ L along with equations (3.72) and (3.73), express the result in terms of renormalized variables and set u = u * , we recover the expansion of Ω(c 1 , c 2 ) to linear order in c j given on the right-hand-side of equation (3.70).
We now insert the scaling form (3.66) of the effective two-point vertex into the expressions (3.62) and (3.63) for f [1] ϕ and f [2] ϕ . Expressing the sum in terms of renormalized variables yields

(3.75)
Here the ellipsis stands for contributions that are O(u * , 2 ) as long as Ω(c 1 , c 2 ) = O( 0 ). This condition is fulfilled whenever c 1 + c 2 > 0. Depending on whether this is the case or not, the function D ϕ (c 1 , c 2 ) has an expansion in integer powers of , namely  (1). The origin of these UV singularities is obvious. The free-energy density f ϕ of the effective (d − 1)dimensional theory at d = 3 involves UV singularities of the form Λ 2 and Ω ln Λ. The bulk and surface counter-terms included so far do not cure these; additional subtractions (of the kind needed for a super-renormalizable effective bulk theory in d − 1 = 2 dimensions, cf. [48]) would be needed. Evaluated at d = 3, the first term on the righthand side of equation (3.75) would then yield a contribution (Ω/8π)[1 − γ E − ln(Ω/4π)] to f ϕ,res,R and hence to D ϕ . Such logarithmic terms are encountered at d = 3 also in Gaussian film models; see e.g. reference [49].
Let us postpone any further discussion of the behaviour of D ϕ at d = 3 for the moment and first compute the free-energy contributions f [1] ψ and f [2] ψ . The former differs from f [1] through the lowest-mode contribution. This is given by equation (3.62) with the shift δτ L set to zero. Hence we have To determine the two-loop graph of f [2] ψ , we must subtract from the two-loop graph of f [2]  and .

(3.82)
To visualise the result, we have plotted in figure 3 the difference between the O( ) expression (3.48) for D(c 1 , c 2 ) and its analogue (3.80) for D ψ (c 1 , c 2 ) one obtains in the scalar case n = 1 upon evaluation at = 1. This difference corresponds to the contribution from the lowest mode m = 1 to the series expansion of D(c 1 , c 2 ) to O( ). It vanishes as (c 1 , c 2 ) approaches the origin and reaches its minimum at the fixed point (c 1 , c 2 ) = (∞, ∞) corresponding to large-scale Dirichlet boundary conditions at both boundary layers B 1 and B 2 .
Combining equations (3.75) and (3.80), we can write the result of the modified RG-improved perturbation theory used in this subsection as (3.83) Figure 3. Contribution from the lowest mode m = 1 to the expansion of the scaling function D(c 1 , c 2 ) for n = 1, evaluated at = 1. This quantity is given by where D ϕ stands for the expression defined through the right-hand side of equation (3.75) in conjunction with equation (3.67), while D ψ represents the O( ) series expansion in equation (3.80) with the expansion coefficients given by equations (3.81) and (3.82), respectively. Note that the result must be utilised with care. It is not suitable for direct evaluation at d = 3 ( = 1). We defer a discussion of this and related issues to section 4.

3.5.2.
Consistency with fractional expansion for c 1 = c 2 = 0 and Ginzburg-Levanyuk criterion Our motivation for working out the modified RG-improved perturbation theory described in section 3.5 was the goal to achieve consistency of the theory for general nonnegative values of the enhancement variables c 1 and c 2 with the approach used in references [46] and [43] to study the zero-mode case c 1 = c 2 = 0 of critical surface enhancements. Since the former approach reduces to the latter when c 1 = c 2 , it is clear that consistency is ensured and the fractional expansion (3.25) must be recovered. To see this explicitly from the results of the previous subsection, recall that κ 1 vanishes as c 1 + c 2 → 0. Therefore, D ψ,0 (c 1 , c 2 ) and D ψ,1 (c 1 , c 2 ) approach the limiting values D 0 (0, 0) and D 1 (0, 0), respectively. Furthermore, the first term on the right-hand side of equation (3.75) yields the 3/2 term according to equation (3.77). That the modified RG-improved perturbation theory is controlled for small , irrespective of whether c 1 + c 2 is positive or zero, can also be seen from a criterion of the Ginzburg type [50,51]. Letγ (4) of the effective four-point vertex γ (4) (y 1 , . . . , y 4 ). The analogue of equation (3.66) tells us thatγ (4) (0, 0, 0) ∼ u * µ 2η L d−5+2η , where the proportionality factor is a function of the scaling variables c 1 and c 2 . Following a standard reasoning (used also by Sachdev in his study of finite-temperature crossovers near quantum critical points [48]), we can construct fromγ (2) Sinceγ (2) (0) is proportional to Ω(c 1 , c 2 ), this measure of the strength of the nonlinearities in H eff behaves for small as Hence the Ginzburg-Levanyuk-type criterion U 1 is satisfied when 1. once the non-universal amplitude E * c (u) has been fixed. As is borne out by figure 4, choices of (c 1 , c 2 ) exist for which these paths start in the region D < 0, enter the region of positive values as L increases and finally return back to the region D < 0. Hence one expects that crossovers of the Casimir force from attractive to repulsive and back to attractive behaviours should occur as a function of L.

Discussion of results and extrapolation to d = 3 dimensions
The result for the scaling function D n=1 (c 1 , c 2 ) of the Casimir force that follows from this extrapolation for D n=1 (c 1 , c 2 ) via equation (3.18) is displayed as a three-dimensional plot in figure 5 and as contour plot in figure 6.
To obtain these plots, we proceeded as follows. We substituted the extrapolated expression [D 0 + D 1 ] =1 for the scaling function D in equation (3.18). For the prefactor d − 1 + y c we used the value 2.718 corresponding to Hasenbusch's recent Monte   Carlo estimate y c (n=1, d=3) 0.718(2) [52]. Let us note that some earlier Monte Carlo calculations led to significantly larger values for this exponent [53,54,55], but others and more recent ones [56,57,58,59] yielded similar numbers. Field-theory estimates based on the expansion to second order and the massive field-theory approach to two-loop order [60,61] gave y c (1, 3) 0.75 and y c (1, 3) 0.85, respectively. For a detailed list of Monte Carlo and field-theory estimates for y c (1, 3) the reader is referred to reference [52]. Choosing a somewhat larger value of y c (1, 3) such as the quoted field-theory estimates would lead to small quantitative, but no qualitative changes in figures 5 and 6.
The behaviour of the plotted function D is qualitatively similar to that of D. As is obvious from figures 5 and 6, the critical Casimir force does exhibit the anticipated crossovers from attraction to repulsion and back to attraction for appropriate choices of (c 1 , c 2 ). According to a theorem for systems satisfying reflection positivity [62], the Casimir force cannot become repulsive for equal enhancements c 1 = c 2 . This conforms with the fact that the Casimir force, by continuity, is attractive (i.e. D < 0) in the vicinity of the 11 diagonal. By choosing (c 1 , c 2 ) sufficiently away from the 11 diagonal, one can ensure that a double sign change of F occurs as L increases. A fascinating consequence is that for such choices of (c 1 , c 2 ), critical values L = L 0 (c 1 , c 2 ) of the film thickness exist for which the critical Casimir force vanishes.
In the extrapolations for D n=1 (c 1 , c 2 ) and D n=1 (c 1 , c 2 ) presented above, we tacitly assumed that c 1 + c 2 > 0. Knowing that the expansion breaks down when c 1 = c 2 = 0 (becoming fractional), we may expect to see indications of this breakdown in the behaviour of these extrapolations near the origin. This is indeed the case, though only very close to it. To demonstrate this, we depict in figure 7   This function [D 0 + D 1 ] n=1 (c, c) has an infinite slope at c = 0. Its asymptotic behavior near the origin can be determined analytically in a straightforward fashion. One finds (4.1) The in-set in figure 7 shows a comparison of [D 0 + D 1 ] n=1 (c, c) and its asymptotic form (4.1) for small c. The term ∝ c 1/2 in equation (4.1) results from that part of the two-loop contribution f [2] res whose mode summation m 1 ,m 2 is restricted to m 1 = 1 in either one of the two loops and to m 2 > 1 in the respective other. The easiest way to recover it is to expand D ϕ (c, c)| n=1 to O( ). The coefficient of the O( ) term is precisely the contribution ∝ √ c in equation (4.1). It is a consequence of the (spurious) infrared singularity one encounters in this expansion due to the m = 1 mode with κ 1 (c, c) → 0, and would imply an infrared singular derivative ∂ c of the function (4.1) at c = 0. This is a spurious infrared singularity since the system is not expected to exhibit critical behaviour at T − T c,∞ = c 1 = c 2 = 0 when L < ∞. The extrapolated scaling function [D 0 + D 1 ] n=1, =1 of the Casimir force behaves analogously near the origin, as should be clear from equation (3.18).
To illustrate that the spurious behaviour of [D 0 + D 1 ] n=1 (c, c) near the origin is due to the zero-mode contribution, we also display the function [D ψ,0 + D ψ,1 ] n=1 (c, c) in figure 7. The latter function is regular at c = 0 and has a finite positive slope there, since it does not involve the zero mode m 1 = 1. Its value at c = 0 agrees with that of The remaining brown curve in figure 7 represents an extrapolation based on equation (3.83), namely the function This means that we have discarded the two-loop contribution to D ϕ and replaced the pre-factor A d−1 /(d − 1) of the one-loop term by its value at d = 4. To understand the rationale for this approximation, one should note the following. If c 1 = c 2 = 0, then this term reduces precisely to the a 3/2 (n) 3/2 contribution to D(0, 0). If we expand, on the other hand, this choice of D ϕ to linear order in when c 1 + c 2 > 0 and substitute it into D ψ,0 + D ψ,1 + D ϕ , we recover the O( ) result D 0 + D 1 . Hence, in the limit c → 0, D app (c, c) approaches the value one obtains for ∆ (sp,sp) | n=1, =1 by evaluating the small-expansion (3.25) at = 1. Furthermore, its derivative ∂ c D app is finite at c = 0. Despite these appealing features, we refrain from showing a plot of this function for all c ∈ (0, ∞) because away from the origin this function deviates considerably from the O( ) extrapolation [D 0 + D 1 ] n=1 . The reason is that Ω(c 1 , c 2 ; = 1) becomes fairly large for large c 1 + c 2 (see figure 2) and the difference between (A 3 /3)Ω(c, c; ) and its extrapolated O( ) expansion (1+ ∂ )(A 3 /3)Ω(c, c; 0)| =1 is not small. Away from the origin, we therefore do not consider D app to be superior to the naïve extrapolation D 0 + D 1 .
As we mentioned earlier, the pre-factor A d−1 /(d − 1), has a UV pole at d = 3. Being guided by the aim to achieve consistency with the small-expansion for c 1 = c 2 = 0, we expanded it in and hence could replace it by its value at d = 4 at the order of our calculation. In an extrapolation based on equation (3.83) one might be tempted to use the value of the UV finite one-loop term of D ϕ one obtains in the manner described below equation (3.77) by subtracting the pole term [Ω(c 1 , c 2 , = 1)] (d−1)/2 and evaluating the difference at d = 3. However, such an approach would mix the RG approach based on the expansion we used for f ψ with elements of a fixed-dimension RG approach. The corresponding extrapolation would again lead to substantial differences from the O( ) extrapolation D 0 + D 1 for large values of c 1 + c 2 . We feel that such attempts to evaluate D ϕ directly at d = 3 should be based on a RG approach in fixed dimension d = 3. Such an approach is known to involve the study of the theory away from the bulk critical point [63,64]. Furthermore, the surface enhancement variables c 1 and c 2 provide two mass parameters in addition to the bulk correlation length one has to deal with. A corresponding massive RG approach in fixed dimension d = 3 based on references [60] and [61] -or combined with elements of the strategy followed in reference [65] -is technically very demanding and beyond the scope of the present work, albeit conceivable.
It is instructive to compare the extrapolated scaling function D n=1 of the critical Casimir force displayed in figures 5 and 6 with their exact analogues res (1; 0, c 1 , c 2 )/n (4.3) of the Gaussian theory in d = 3 and d = 4 dimension. Here, the one-loop term (3.41) must be substituted for f [1] res . A plot of these functions on the diagonal c 1 = c 2 ≡ c is shown in figure 8. As one sees, the extrapolation D n=1 (c, c) lies for most values of c below and above the Gaussian scaling functions (4.3) with d = 4 and d = 3, respectively. Having already discussed the spurious small-c behaviour of the extrapolation, let us add a comment about the limiting behaviour at large c. It has been known for long that deviations of the surface enhancement variables c j from their fixed point values c j = ∞ associated with the ordinary surface transition correspond to irrelevant surface scaling fields ("extrapolation lengths") whose RG eigenvalues are given by the momentum dimension −1 [66,22]. Taking into account the corrections to scaling ∝ L −1 implied by them is crucial for a proper analysis of Monte Carlo simulation results for Casimir forces [26,27,28,29,30]. Our small-results enable us to explicitly verify the presence of such corrections to scaling. In order to comply with a correction to scaling to the residual free energy that is down by a factor 1/L, the scaling function D(c 1 , c 2 ) ought to exhibit the asymptotic behaviours respectively. Let us recall that the behaviour of the extrapolation D n=1 plotted in figure 8 was obtained by substituting [D 0 + D 1 ] =1 into equation (3.18). Therefore, the plotted extrapolated function does not approach its limiting value at c = ∞ as c −1/yc . This is because its asymptotic form contains the logarithmic term ∼ ln(c)/c mentioned above, evaluated at = 1. Hence, the deviation of the plotted extrapolation D n=1 from its value at c = ∞ varies (incorrectly) as ∼ ln(c)/c. For conciseness, we have refrained from designing alternative extrapolations into which the correct asymptotic large-c asymptotics is incorporated via appropriate exponentiation ansatzes.

Summary and Conclusions
Considering an O(n)-symmetric φ 4 model in film geometry, we investigated the fluctuation-induced forces produced by thermal fluctuations between the boundary planes of a film that is held at its bulk critical point. We restricted our attention to the case of generic non-symmetry-breaking boundary conditions. The interaction constants c j of the O(n) symmetric boundary termsc j φ 2 /2 included in the Hamiltonian were chosen to correspond to non-supercritical surface enhancements but otherwise allowed to take arbitrary values. Together with the restriction to the bulk critical temperature T c,∞ (or, more generally, to temperatures T ≥ T c,∞ ), this requirement guarantees that the O(n) symmetry is neither explicitly nor spontaneously broken. Whenever the symmetry φ → −φ is broken, the order-parameter profile φ does not vanish. As discussed in the introduction, this implies that a nonzero Casimir force is obtained even in mean-field theory -i.e., in the absence of fluctuations. Thus in these cases the Casimir force is not fully fluctuation induced but contains a non-fluctuating component. This applies, in particular, to binary liquid mixtures, for which beautiful direct measurements of the Casimir forces were accomplished recently [16,29,67]. A key motivation for excluding such symmetry breakdowns in the present investigation was the intention to keep the thermodynamic Casimir forces entirely fluctuation induced, making them share this crucial feature with the original QED Casimir force [1,4].
Our analysis was based on the field-theoretic RG approach in 4 − dimensions. The use of the expansion in studies of critical and near-critical Casimir forces has shortcomings that go beyond the familiar ones known from its application to critical behaviour in bulk and semi-infinite systems. The series expansions it yields for bulk and surface critical exponents as well as for other properties at bulk critical points -Taylor expansions in -are asymptotic. Naive extrapolations of such expansions to low orders in are not normally quantitatively reliable. However, quantitatively accurate results can be obtained by extending the expansions to sufficiently high orders, combining them with information from appropriate large-order results and using appropriate Borel-Padé resummation techniques [68,64,69,70]. The additional complication one encounters in the case of films is that = 4 − d ceases to be an appropriate expansion parameter when dealing with dimensional crossovers. (The adequate analogue of this parameter for the study of critical behaviour at the transition temperature T c,L of a film of finite thickness L and infinite extension in d−1 Cartesian directions would be 5 ≡ 5−d.) The limitations of the approach become particularly apparent when zero modes appear at zero-loop order at the bulk critical point, as happens for critical enhancement c 1 = c 2 = 0 of both boundary planes. Whenever this occurs, the expansion breaks down at the bulk critical temperature T c,∞ and becomes fractional, involving also powers of ln [43,46,47,48]. By contrast, the expansion remains (formally) intact at T c,∞ provided at least one of the (non-supercritical) enhancements c j is subcritical, although a similar breakdown clearly must occur at a temperature below T c,∞ .
Despite these deficiencies, the field-theoretic RG approach in 4 − dimensions has a number of appealing and valuable features. First of all, it lends itself to nontrivial checks of the scaling forms the residual free energy and the Casimir force are expected to have according to finite-size scaling theory and the RG approach. Since we performed a two-loop calculation which gave f res to order u, we were able to identify the ln(µL) contributions implied by the nontrivial O( ) term of the surface crossover exponent Φ in the scaling arguments (3.13) and thus verify the scaling forms (3.13) and (3.17) to first order in . Second, the fact that general properties such as the exponential decay in the disordered phase (see references [20,Sec. VI] and [43, Sec. IV.C]) and analyticity properties of the free-energy density f (see references [20,Sec. VII] and [43, Sec. IV.C]) must hold at any order of the small-expansion provides nontrivial checks of them. For the sake of clarity, it will be helpful to recall these analyticity properties. For given n and sufficiently large d, the film undergoes a sharp transition to an ordered phase at a temperature T c,L < T c,∞ . In the case of non-supercritical surface enhancements considered here, there should also be no surface transitions of the film at temperatures T ≥ T c,∞ (i.e. its surface layers should not exhibit long-range order). Consequently, the total free-energy density f must be analytic in temperature at T c,∞ . This requirement imposes conditions on the behaviour of the temperature dependent analogue of the scaling function of f res . We refer the reader to references [20,Sec. VII] and [43, Sec. IV.C] for details.
These analyticity requirements were found to be fulfilled by the two-loopexpansion results for anti-periodic, (O, O), and (O, sp) boundary conditions, but violated at O( ) for periodic and (sp, sp) boundary conditions [20]. Extension of the smallexpansion to O( 3/2 ) in the latter two cases via the effective-action approach outlined above increased the order at which a violation of analyticity requirements occurred to 3/2 [43]. Third, our analyses in the latter reference and above yield evidence of a shift of the point c 1 = c 2 = 0 at which the zero-mode field ϕ becomes critical at T c,∞ to an L-dependent location with c 1 + c 2 < 0. By analogy with the shift T c,∞ → T c,L of the transition temperature, such a shift is expected and in conformity with the required analyticity of f at bulk criticality.
One crucial problem one is faced with in studies of critical behaviour and dimensional crossover in films is the difficulty of obtaining accurate results for such shifts in d = 3 bulk dimensions by analytical methods. The reason is that fluctuations on all length scales between the microscopic one (lattice constant) and the bulk correlation length may contribute, as a result of which these shifts may acquire contributions that are non-analytic in the interaction constant [71,72,60,73]. Since we used RG-improved perturbation theory in 4− dimension at T c,∞ to determine the effective action H eff , our approach is not capable to capture such non-analytic contributions. It gave us a shift along the c 1 +c 2 axis of the point where the ϕ component of the order parameter becomes critical that was proportional to u * . This in turn implied a contribution to the critical Casimir force for critical surface enhancements c 1 = c 2 = 0 of the form (u * ) (3− )/2 . To make our results for the scaling functions D(c 1 , c 2 ) and D(c 1 , c 2 ) consistent with the fractional expansions at (c 1 , c 2 ) = (0, 0), we extended the effective-action approach of references [46] and [43] to nonzero values of c 1 + c 2 . We succeeded in achieving consistency with the results of these references for the case c 1 = c 2 = 0 and the fractional expansions. Furthermore, the effective-action approach turned out to be well-behaved for small . Nevertheless, its results must be taken with a grain of salt: they did not lend themselves easily to simple extrapolations to d = 3 dimensions that could be judged as clearly superior to those based on the conventional expansion.
Let us emphasise that our results exhibit a number of interesting and important qualitative features which may be expected to persist in quantitatively more accurate investigations, irrespective of the accuracy of our extrapolation to d = 3.
(i) Generalising previous work [19,20,37,43,47], they clearly demonstrate that critical Casimir forces can be attractive or repulsive depending on the values of the surface enhancement variables c j , even when the internal symmetry (Z 2 for n = 1 and O(n) for n > 1) is neither spontaneously nor explicitly broken.
(ii) They are in full accordance with -and hence corroborate -the scaling forms (3.14) and (3.17) of the residual free energy and Casimir force, respectively. This means that the Casimir amplitude becomes scale dependent.
(iii) As a further important observation we found that crossovers from attractive to repulsive forces and vice versa can occur if the surface interaction constants take appropriate values. In conjunction with (ii) this means that the critical Casimir force goes through zero at a certain thickness L 0 .
(iv) For the case of primary interest -the Ising bulk universality class with short-range interactions and d = 3 -we saw that for properly chosen surface enhancement values c 1 and c 2 , crossovers from attractive to repulsive and back to attractive behaviours should occur as L increases. It would certainly be worthwhile to check this prediction along with (iii) by Monte Carlo calculations.
(v) Finally, let us mention that crossover behaviours analogous to (iii) and (iv) should also be possible for three-dimensional O(n) systems with easy-axis spin anisotropies. Semi-infinite systems of this kind are known to have anisotropic analogues of the isotropic special transition [74,75], characterised by the coincidence of the transition temperature at which the easy-axis component of the order parameter at the surface becomes critical with the transition temperature T c,∞ of the n > 1 bulk system. Hence for finite thicknesses L, the surface interaction constants associated with a given easy axis on one or both boundary planes can be critically enhanced. This suggests that the mentioned analogues of (iii) and (iv) should occur. (Films involving different directions of easy axes on the two boundary planes would require separate analyses.) Appendix A. Calculation of the one-loop free-energy density f [1] Appendix A.1. Calculation based on the Abel-Plana summation formula To compute f [1] , we start from equation (3.28), differentiate this expression with respect toτ and compute the integral over p. This gives The expression in the second line of equation (A.1) varies ∼ L ∞ 0 dk (τ + k 2 ) (d−3)/2 in the large-L limit. We could subtract this bulk term to make the difference well defined. However, we are ultimately interested in the free-energy density f [1] (L;τ ,c 1 ,c 2 ). Term-by-term integration of the series for ∂τ f [1] increases the exponent (d − 3)/2 of the series coefficients by one. Hence, subtracting the limiting bulk contribution Lf [1] b would not suffice to render UV finite expressions for the free-energy density f [1] (L;τ ,c 1 ,c 2 ) in dimensions d 4. To obtain well-defined results for this quantity, we prefer to use analytical continuation in d, rather than making additional subtractions. We start by integrating the series in equation (A.1) term-wise, choosing the integration constant such that That the integration constant has been chosen correctly may not be obvious at this stage. However, once we will have analytically continued the above series (A.2) in d, one can confirm this choice by verifying that the resulting free-energy expressions reduce to known dimensionally regularised results atτ = 0+ in the limitsc 1 ,c 2 → 0 andc 1 ,c 2 → ∞. Furthermore, there are two other ways to arrive at expression (A.2), both of which indicate its consistency with general rules of dimensional regularisation. The first is to use the representation for the logarithm in the second line of equation (3.28), perform the integrations and take into account that integrals of pure powers such as (d−1) p 1 vanish in dimensional regularisation. Another way is to insert the partial-p identity (d − 1) −1 ∇ p · p into the integrand of the integral (d−1) p in the last line of equation (3.28) and integrate by parts, dropping the boundary terms at p = ∞.
According to the result (A.2) we must calculate -and analytically continue -a series of the form where κ m are the positive solutions to equation (2.21). Such series are generalisations of Epstein-Hurwitz ζ functions (see e.g. reference [76]). Let us assume that both C j ≥ 0 and C 1 + C 2 = 0. Then all κ m , m = 1, 2, . . . , ∞, are positive, and as we showed in section 2, non-degenerate. Hence the function ∂ κ ln R C 1 ,C 2 (κ) has a simple pole with residue 1 at each of these κ m . Since (κ 2 + b 2 ) a is regular at κ = κ m , the value of this function at κ m is given by the residue Res . But at κ = κ m R C 1 ,C 2 can be factorized as (A.7) Thus we have Since this function and the one appearing in S C 1 ,C 2 (a; b) are even in κ, we can extend the summation in equation (A.4) to all integer-valued m = 0, using κ −|m| = −κ |m| , and divide by 2. Since (κ 2 + b 2 ) a is regular at κ m , each term of this series can be expressed as (2πi) −1 times an integral dκ (κ 2 + b 2 ) a Υ C 1 ,C 2 (κ) along a contour that passes once around the pole at the respective κ m in a counter-clockwise fashion and contains no other singularities. If a < −1/2, the series S C 1 ,C 2 (a; b) converges and the integrand of the contour integral decays sufficiently fast as κ → ±∞ ± i0, so that the union of all these contour integrals can be deformed into a path γ 1 encircling all poles κ m with 0 = m ∈ Z (see figure A1). We now add and subtract integrals along the paths γ 2 and γ 3 depicted in this figure. Since the integrand is regular at all nonzero κ inside the region bounded by the closed path γ 1 ∪ γ 2 ∪ γ 3 , the integral along it is −2πi times the residue of the integrand at κ = 0. Furthermore, the integral along γ 2 approaches zero as the radius of the circle on which γ 2 is located tends to infinity. Hence we have (A.10) Figure A1. Paths in the complex κ-plane.
The factor (κ 2 + b 2 ) a implies that the integrand has branch cuts from ib to i∞ and from −ib to −i∞. In addition, the function g C 1 ,C 2 (κ) and hence Υ C 1 ,C 2 (κ) have poles at ±iC 1 and ±iC 2 , marked by yellow pentagons in figure A1. We move the path γ 3 infinitesimally close to the imaginary axes, passing around the poles ±C j along semicircles of radius δ → 0. The portions of γ 3 with Im κ > 0 and Im κ < 0 give identical contributions. Taking into account that the limiting values of the power on the right and left rim of the branch cut are lim Re κ→0± Im κ>0 one arrives at (A.12) Here P ∞ b dt means the principal value lim δ→0+ C>+δ dt, where C < and C > are the smaller and larger one of C 1 and C 2 , respectively, and b is assumed to be smaller than C < . The contribution in the second line results from the integrals along the semicircles.
In order that the integral P ∞ b dt exists, we must require a > −1 (to ensure convergence at the lower limit of integration) in addition to the original condition a < −1/2. In view of the behaviour of the integrand near the upper integration limit, we must have a < −1/2 (to guarantee convergence at the upper limit of integration). However, we can split off the limiting value lim κ→±i∞ N C 1 ,C 2 (κ)/R C 1 ,C 2 (κ) = ∓i, obtaining Since the first term vanishes at κ = iC j , the contribution it yields to the term in the second line of equation (A.12) vanishes, i.e. Res contributions to S C 1 ,C 2 (a; b) implied by the term −i in equation (A.13) can be represented as an integral along the original path γ 1 , namely where sgn(x) means the sign function. For a < −1/2, the latter integral is well defined and can be analytically computed. This gives where 2F1 , the regularised hypergeometric function, is an entire function related to the standard hypergeometric function 2 F 1 via 2F1 (α, β; γ; z) = 2 F 1 (α, β; γ; z)/Γ(γ). where I C 1 ,C 2 (a; b) represents the meromorphic function defined by the right-hand side of equation (A.15).

Equation (
A.17) provides the appropriate generalisations of the results of reference [20] for the special cases (C 1 , C 2 ) = (0, 0), (∞, ∞) and (0, ∞) to general nonnegative values of C 1 and C 2 . To check the consistency with Krech and Dietrich's results, let us work out the asymptotic behaviours in these limits. Straightforward analysis yields  19) and (A.20) approach infinity in the considered limits C j → ∞. These do not appear if one sets C 2 = ∞ and C 1 = C 2 = ∞ from the outset because the operations of analytic continuation and of taking these limits do not commute. (As long as a < 0, the subtracted terms approach zero as C j → ∞.) Note also that these terms contribute only to the excess surface densities f [1] s,j , but not to f [1] res . We now insert the above results (A.15) and (A.17) into equation (A.2), set a = (d − 1)/2, and take the limit b → 0 (i.e.τ → 0). The t-integral of equation (A.17) becomes ∞ 0 dt π g C 1 ,C 2 (it) t d−1  for such values of d.