Fluctuations of the critical Casimir force

The critical Casimir force (CCF) arises from confining fluctuations in a critical fluid and thus it is a fluctuating quantity itself. While the mean CCF is universal, its (static) variance has previously been found to depend on the microscopic details of the system which effectively set a large-momentum cutoff in the underlying field theory, rendering it potentially large. This raises the question how the properties of the force variance are reflected in experimentally observable quantities, such as the thickness of a wetting film or the position of a suspended colloidal particle. Here, based on a rigorous definition of the instantaneous force, we analyze static and dynamic correlations of the CCF for a conserved fluid in film geometry for various boundary conditions within the Gaussian approximation. We find that the dynamic correlation function of the CCF is independent of the momentum cutoff and decays algebraically in time. Within the Gaussian approximation, the associated exponent depends only on the dynamic universality class but not on the boundary conditions. We furthermore consider a fluid film, the thickness of which can fluctuate under the influence of the time-dependent CCF. The latter gives rise to an effective non-Markovian noise in the equation of motion of the film boundary and induces a distinct contribution to the position variance. Within the approximations used here, at short times, this contribution grows algebraically in time whereas, at long times, it saturates and contributes to the steady-state variance of the film thickness.


I. INTRODUCTION
A fluid close to a continuous phase transition exhibits remarkable universal properties, such as long-ranged fluctuations and scale invariance [1]. A core element of criticality is the notion of an order parameter (OP) φ, which, e.g., in the case of a binary liquid mixture is proportional to the deviation of the concentration c A of species A from its critical value c A,c , i.e., φ ∝ c A − c A,c . Confining a critical fluid leads to a critical Casimir force (CCF) acting on the confining boundaries [2][3][4]. Such a situation arises, e.g., when two colloidal particles immersed in a critical solvent come close to another or if one colloidal particle approaches a container wall of the sample. CCFs thus provide means to control the aggregation behavior of colloidal suspensions [5] and are, accordingly, not only of theoretical but also of highly practical interest. Consequently, CCFs and, more generally, Casimir-like fluctuation-induced forces have been extensively studied in as well as out of equilibrium (see, e.g., Refs. [6][7][8][9][10][11] and references therein).
Previous studies mainly focused on the mean value of the thermally averaged CCF, which is finite and universal, i.e., it does not depend on the specific material considered, but only on the bulk and surface universality class [12,13]. However, since the CCF arises from fluctuations of the OP, it is itself a fluctuating quantity. In Ref. [14] the equilibrium fluctuations of the CCF in a generic Gaussian medium subject to Dirichlet BCs have been investigated theoretically, and it has been shown that the static equilibrium variance depends in a non-universal way on the microscopic details of the model (specifically, on the microscopic length scale below which the continuum description breaks down). These findings were confirmed analytically and via Monte-Carlo simulations of lattice models for the case of periodic BCs [15]. The variance of the CCF acting on membrane inclusions has been found to be strongly cutoff dependent as well [16]. In order to understand these results, we note that the (instantaneous) CCF K (per area) is given by the difference between the pressures acting on the two sides of a boundary, K = P f − P b , where the subscripts refer to f ilm and bulk (a notion to be specified below; see also Fig. 1). Concerning the mean value K , all cutoff-dependent quantities present in P f,b cancel in the difference, giving rise to a finite and in fact universal quantity. In contrast, the variance (∆K) 2 = ∆P 2 f + ∆P 2 b (where ∆K ≡ K − K ) is the sum of the individual variances of P f and P b (assuming statistical independence of film and bulk) and thus no cancellation of divergences can occur [14].
Notably, the detailed form of the cutoff dependence of the variance implies a divergence in the strict continuum limit. The corresponding huge variances predicted in this way for the CCF seemingly stand in contrast to the fact that experiments on CCFs in wetting films [17][18][19][20] or in colloidal systems [5,21] did not observe such giant fluctuations. Colloidal particles in critical bulk fluids rather show essentially standard Brownian diffusion with an effective diffusivity modified due to the fluctuations of the field [22][23][24]. Moreover, near a wall, colloids diffusively explore the corresponding effective potential induced by the mean CCF [25,26].
Divergent equilibrium variances of the force have also been obtained for the quantum-electrodynamical Casimir effect [27]. However, it has been shown that the divergence problem is less severe or even absent when considering temporally averaged quantities [14,27]. In fact, such a temporal average is indispensable due to the necessarily finite temporal resolution of a measurement apparatus. This issue can also be appreciated in view of a Langevin description of a Brownian particle [28]: while the random force has correlations proportional to a Dirac δ function and thus is strongly divergent, the resulting observable quantities, such as the position distribution of the particle, are finite [22][23][24]. We finally remark that the static variance of thermal van der Waals forces between dielectric slabs has been found to be finite and independent of a microscopic cutoff, without the necessity for temporal averaging [29].
In the present study, based on the non-equilibrium stress tensor formalism [9,[30][31][32], we investigate static and dynamic equilibrium fluctuations of the CCF in film geometry for various non-symmetry breaking BCs. We focus on the Gaussian approximation. Simulations of lattice models [15] have shown that it constitutes an accurate description of the (static) stress probability distribution. In accordance with Refs. [14,15] we confirm the generic non-universal and divergent character of the static variance of the CCF acting on a fixed boundary. Going beyond previous studies, we show that the dynamic correlations of the CCF are finite and decay algebraically in time with exponents depending only on the spatial dimensionality and the dynamic universality class. We furthermore consider a movable film boundary subject to the fluctuations of the CCF. We find that the mean-squared displacement of the boundary is, at short times, characterized by an algebraic growth (with possible logarithmic corrections), while, at long times, it saturates at a finite value due to the confining potential stemming from the CCF. Accordingly, the strong divergences of the static variance of the CCF as function of the cutoff do not show up in its dynamic correlations or in experimentally observable quantities such as the position of the film boundary.
We present our study as follows: In Sec. II, we introduce the model and briefly review the formalism required to determine time-dependent CCFs. In Sec. III, we analyze static and dynamic fluctuations of the CCF in the case of a film with fixed boundaries. This constraint is released in Sec. IV, where we consider the thermal fluctuations of the film thickness. A summary of the findings of the present study is provided in Sec. V.
with (∇φ) 2 = d α=1 (∂ α φ) 2 , ∂ α = ∂ ∂rα ≡ ∇ α , and the temperature parameter which is directly proportional to the reduced temperature τ ≡ T − T c T c > 0, (2.5) where T c denotes the (bulk) critical temperature. Furthermore, ξ with a bulk critical exponent ν. Within the Gaussian approximation, one has ν = 1/2, such that Eq. (2.6) reduces to ξ = 1/ √ τ . For typical simple fluids one finds ξ (0) + ≈ 1 nm [19]. For the surface Hamiltonian density H s we do not assume a specific form, except that it is strongly localized at a boundary; e.g., for a planar surface at z = L, it takes the form H s (r, L, φ, ∇φ) = δ(z − L)U s (φ, ∇φ) (2.7) with a potential U s . We consider the thin film geometry with macroscopically large lateral directions and, accordingly, we decompose the d-dimensional vector r as r = {r = (r 1 , . . . , r d−1 ), z}. The following BCs of the OP are considered in the transverse direction z: and Neumann : Denoting the physical time byt [to be distinguished from the rescaled one, see Eq. (2.11) below], the OP dynamics is taken to be described by model B [33]: where γ is a mobility coefficient and is the chemical potential. Far away from a boundary, one has µ = −∇ 2 φ + τ φ. Furthermore,ζ is a Gaussian white noise with zero mean and correlations ζ ( The form of Eq. (2.9) and of the noise ensure local conservation of the OP in an unconfined system. In order to guarantee local and global OP conservation in the presence of boundaries, the flux J = −γ∇µ + N of the OP (with noise N such that ∇ · N =ζ) must vanish at the confining surfaces of the system, i.e., J z = 0 for z = 0, L [34]. This is ensured by periodic or Neumann BCs, which we shall henceforth use in the dynamic model [35]. In order to simplify the notation, we follow Ref. [9] and remove the temperature T from the description by introducing a rescaled OP field φ/T 1/2 and, additionally, introduce a rescaled time with a Gaussian noise ζ ≡ζ/T 1/2 γ correlated as The instantaneous CCF alluded to in the introduction is defined in terms of a dynamic stress tensor which can be derived for generic Landau-Ginzburg-type of models [30]. The associated formalism has been developed further in Ref. [31] and has been applied in Refs. [9,32] for the case of a quench of a conserved fluid at criticality. Notably, the expressions for the stress tensor in and out of equilibrium are, in general, different [36]. The concept of a dynamic stress tensor allows one to rigorously define an CCF at any instant during the time evolution of the system. The determination of time-dependent forces on boundaries within model B is detailed in Refs. [9,32]; in the following we provide only the essential equations. The dynamical stress tensor is given by [32] T αβ ≡ T αβ + µφδ αβ , (2.14) FIG. 1. Sketch of the situation considered in the present study. A film of thickness L is enclosed by two boundaries (thick black lines) each of surface area A (spanning the plane perpendicular to the figure). The right boundary is enclosed by a thin imaginary volume Vs. The CCF acting on the right boundary is obtained from an integration of the stress tensor over the surface of Vs, which consists of an inner (∂Vs,i) and an outer part (∂Vs,o) [see Eq. (2.18)]. The film and the bulk consist of the same liquid whereas the gray area corresponds to a solid.
which is expressed here in terms of the standard (grand canonical) stress tensor The last term in Eq. (2.14) involving the chemical potential µ [Eq. (2.10)] generically arises in a non-equilibrium situation for a system with no-flux BCs. Considering the right boundary of the film at position z = L as shown in Fig. 1, the instantaneous generalized force K α (per area A) in direction α of the surface is defined as which can be transformed into [32] [9]). Since the resulting expression for the CCF is identical to omitting the last term in Eq. (2.17), we do not separately deal with this aspect. In general, the instantaneous CCF (per area) K ≡ K z acting on a boundary at position z of the film is obtained as represents the film and the bulk pressure; ∂V s,i/o , correspondingly, denotes the bounding surface of the inner/outer fluid and r is the location of the (d − 1)-dimensional surface element ds z . We emphasize that the quantity K is actually a pressure. However, to stay in line with the convention employed in most of the literature, we use the term Casimir force throughout this study. In equilibrium, we assume the boundaries of the film to be fixed at z ∈ {0, L}. However, in order to facilitate later extensions of the model, in our formalism we shall keep z general. For the Landau-Ginzburg Hamiltonian density of Eq. (2.3) and upon using Eq. (2.10), the dynamical stress tensor takes the formT which, as noted after Eq. (2.17), depends only on bulk terms. The mean CCF (per area) acting on a surface at location z follows from Eq. (2.18) as where "{bulk}" denotes the corresponding expression evaluated in the bulk, which is obtained by integrating over ∂V s,o and calculating statistical averages with the corresponding bulk distribution of φ.

C. Order-parameter correlation functions
Here, we summarize some useful expressions for and properties of the equilibrium two-point two-time OP correlation function C(r, r , t) ≡ φ(r, t)φ(r , 0) = C(r −r , z, z , t), where the latter equality follows from translational invariance in the lateral directions. The time-independent, static limit of C is denoted by C st and follows as C st (r , z, z ) = C(r , z, z, t → 0).

Bulk correlation function
The equilibrium bulk correlation function for model B is given by (see, e.g., Ref. [9]) which, as indicated by the notation, is spatially isotropic [37]. In the static equilibrium case (t = 0) for d > 2 and at bulk criticality (τ = 0), one has In d = 2 dimensions, the static bulk correlation function diverges logarithmically for large r [9].
We now analyze C b [Eq. (2.22)] for d > 2 in certain relevant limits. Due to the exponential in the expression for C b , the integrand in Eq. (2.22) gives significant contributions only if q 4 t + q 2 τ t 1. There are two characteristic asymptotic regimes depending on q 2 τ or q 2 τ , as can be seen by splitting the integral in Eq. (2.22) accordingly: (2.24) If τ t 1/2 1, i.e., for long times t τ −2 , the second term on the r.h.s. of Eq. (2.24) is exponentially suppressed, and we can obtain an asymptotic estimate of C b by extending the integration range in the first term to all q: On the other hand, for τ t 1/2 1, i.e., at short times t τ −2 , the second term in Eq. (2.24) dominates, which can be shown as follows: for τ t 1/2 1, the first term in Eq. (2.24) can be approximated as stands for an approximation of the Dirac δ function. The ratio of the first and the second term follows as (τ 2 t) d/4−1/2 /τ d/2 . Accordingly, for d > 2, the first term is indeed negligible at short times. The bulk correlation function can be asymptotically estimated by extending the integration in the second term in Eq. (2.24) over the whole range ofq, resulting in [9] , where 1 F 3 is a generalized hypergeometric function [38] and, according to Eq. (2.11), the scaling variable ψ is dimensionless.

Film correlation function
In a thin film with periodic BCs, translational invariance implies C (p) (r , z, z , t) = C (p) (r , z−z , t) and the two-time correlation function is given by For Neumann and Dirichlet BCs, respectively, one has 10)], we shall consider these BCs only in the static limit given by setting t = 0 in Eq. (2.30). By means of the Poisson resummation formula, the above correlation functions can be expressed in terms of C b as (see, e.g., Refs. [9,39]) with a scaling function C, the explicit expression of which is not needed in the following but can in principle be obtained from Eqs. (2.26) and (2.32).

Further properties and relations
We close this section by collecting a number of relations useful for later developments in this study. Owing to the linearity of the statistical average, one has in the bulk as well as in a film. Moreover, for α ∈ {1, . . . , d − 1}, one has where we used translational invariance in the lateral directions in order to obtain the last equation.

III. FLUCTUATIONS OF THE CCF
In equilibrium, the two-time correlation function of the instantaneous CCF follows from Eq. (2.18) and it is given by where we used the fact that film and bulk pressure are uncorrelated, since they emerge from physically separated parts of the system. Accordingly, the fluctuations of the CCF (per area) given by are correlated in equilibrium as with ∆P(z, t) ≡ P(z, t) − P(z, t) and where we used Eq. (3.1) as well as Equation (3.4) reflects the fact that the variance of a sum of two uncorrelated random processes (P f and −P b , see Eq. (2.18)) is the sum of the individual variances. For t − t → ∞, P f (z, t) and P f (z , t ) become uncorrelated in equilibrium, i.e., P f (z, t)P f (z , t ) → P f (z) P f (z ) ; this applies analogously to P b . Accordingly, one has ∆K(z, t)∆K(z , t ) Since we consider equilibrium dynamics, we can exploit time-translation invariance and henceforth set t = 0.

A. Pressure correlation function
Equation (3.4) can be evaluated by using Eqs. (2.18)-(2.20) as well as the fact the four-point correlation function for a multivariate zero-mean Gaussian process X n ∼ φ or X n ∼ ∂ β φ is given by the following cumulant expansion: (3.7) Furthermore, we use the fact that A d d−1 r A d d−1 r C(r − r ) = A A d d−1 r C(r ), which follows upon applying the (volume preserving) coordinate transformation (r , r ) → (r − r , r ) and by using translational invariance in the lateral directions. In order to illustrate the calculation of P f (z, t)P f (z , 0) , we consider the term where ds z (r ) denotes an area element and represents a contribution to P f (z) P f (z ) , resulting from the first term on the r.h.s. of Eq. Altogether, we obtain where ∂ γ ∈ {∂ r1 , . . . , ∂ r d−1 , ∂ z } and we define∂ γ ∈ {∂ r1 , . . . , ∂ r d−1 , ∂ z }, such that∂ γ = ∂ γ = ∂ rγ for γ = 1, . . . , d − 1 and∂ d = ∂ z ; for brevity, we suppress the arguments of C(r, z, z , t) on the r.h.s. The same expression results for P b (z, t)P b (z , 0) , but with C replaced by the bulk correlator C b . It is useful to note that the quantity ∆P f (z, t)∆P f (z, 0) = P f (z, t)P f (z , 0) − P f (z) P f (z ) has the length dimension of A −1 L −(d+1) , which can be readily inferred from Eqs. (2.27), (2.30), and (3.10). We now specialize Eq. (3.10) to various BCs in thin films and make use of the available spatial symmetries. For Neumann and Dirichlet BCs we will obtain expressions in terms of the correlation function for periodic BCs. The film pressure variances are evaluated further in Secs. III B and III C.

Bulk and periodic BCs
In the following, we jointly analyze the pressure both for the bulk system and for thin films with periodic BCs in the transverse direction. We note that translational invariance [see Eq.
with the analogous expression applying in the bulk. Upon using Eqs. (2.36) and (3.11), the mean equilibrium CCF for a film, with one surface at z = 0 and with periodic BCs, follows from Eq. (2.21) [9]: The film pressure correlator in Eq. (3.10) reduces to where we again suppressed the arguments of C(r, z − z , t) on the r.h.s. In the special case z = z = 0, Eq. (3.13) can be simplified using Eq. (2.35): . (3.14) For bulk pressure correlations P b (z, t)P b (z , 0) , the same expressions as in Eqs. (3.13) and (3.14) apply, but with C (p) replaced by the bulk correlation function C b .

Neumann BCs
The pressure correlation function for a film with Neumann BCs is given by Eq. (3.10) with C replaced by C (N) [see Eq. (2.32b)]. In the special case z = z = 0, the resulting expression can be simplified by noting that Eq. (2.32b) implies Using these relations in Eq. (3.10) renders , (3.16) where, as indicated on the r.h.s., the correlation function C (p) for periodic BCs has to be evaluated for a film thickness of 2L.

Dirichlet BCs
We begin by noting that Eq. (2.32c) implies Using these relations to evaluate the film pressure correlator in Eq. (3.10) for z = z = 0 renders The equilibrium variance of the CCF at a fixed boundary (located at z = z = 0) follows by setting t = t = 0 in Eq. (3.4), resulting in This expression is evaluated in the following by replacing in Eq. (3.10) C by the expression of C st for the respective BCs. Below we discuss the final results at bulk criticality (τ = 0).

Bulk system
The static equilibrium variance of the bulk pressure ∆P 2 b ≡ P 2 b − P b 2 is obtained by inserting the equilibrium bulk correlation function C b,st [Eq. (2.23)] for C into Eq. (3.10). In order to simplify this calculation, we note that, due to the translational invariance of C b , ∆P 2 b takes the same form as the variance for periodic BCs reported in Eq. (3.14) (with C (p) replaced by C b ). We furthermore recall the Schwinger-Dyson relation (see, e.g., Ref. [9]) which implies the following identity for the equilibrium bulk correlation function: Using this relation in Eq. (3.14) renders this end, we evaluate their Fourier space representation by taking into account a wavenumber cutoff q max 1/ε and, accordingly, regularize the δ function as [14] δ(z) = Analogously, based on the Fourier representation in Eq. (2.22), we obtain the regularized form of the bulk correlator (at τ = 0) The term −n · [. . .] ∂A δ(z) in Eq. (3.22) vanishes because the boundary ∂A of the surface is assumed to be located basically at infinity, such that r ,α = 0. Altogether, the variance of the bulk pressure at bulk criticality (τ = 0) turns out to be This expression is singular at d = 2 and positive for all other integer dimensions d.

Periodic BCs
According to Eq. (2.32a), the static OP correlator for a thin film with periodic BCs follows as st . We note that in Eq. (3.26) the bulk correlation function corresponds to the term with m = 0.
Using Eq. (2.32a) in Eq. (3.21) yields the following identity for the film correlation function: which can be evaluated in terms of the Riemann zeta function ζ. For m 1 = 0 or m 2 = 0, the integral over A does not converge at its lower limit (r → 0) and is thus evaluated using a small-distance cutoff ε for the radial coordinate r . The integral over A generally converges at its upper (large distance) limit. Altogether, we obtain the following film pressure variance for periodic BCs with τ = 0: The power laws for L and ε emerge by using the scaling form for C b,st given in Eq. (3.26). It turns out that (∆P , which is only a function of d and the dimensionless parameter L/ε, remains positive for all d > 2 and all L/ε > 0. In the relevant regime L/ε 1, the variance (∆P

Neumann BCs
Equation (3.16) renders, after some algebra analogous to the one leading to Eq. (3.28), the static variance where ε denotes a small length scale required for regularizing the integral over the surface area. The term ∝ 1/(Aε d+1 ) scales like the bulk pressure variance [see Eq. (3.25)], but its detailed form differs from the one in Eq. (3.28) due to the different structures of Eqs. (3.14) and (3.16).

Dirichlet BCs
Upon introducing, as above, a small-distance cutoff ε at the lower limit of the areal integral in Eq. (3.18), the static pressure variance follows as [40] (∆P The universal subdominant term (first term in the square brackets) coincides with Eq. (14) in Ref. [14], as does the dominant (bulk-like) scaling behavior for small ε [given by the last term in Eq. (3.30) and by Eq. (13) in Ref. [14]].

Discussion
According to Eq. (3.19), the variance of the CCF (per area A and thermal energy k B T ) is given by the sum of the bulk and film pressure variances. The film and bulk pressure variances determined above [see Eqs. (3.25) and (3.28)-(3.30)] pertain to a thin film geometry (A = L d−1 → ∞) and thus represent the leading terms in an expansion in terms of 1/A (for ε nonzero). We thus expect the above results to represent a reasonable estimate for systems with a sufficiently large aspect ratio L /L, such that the lateral OP modes approximately form a continuum spectrum. Thus, in the limit of small ε → 0, the static equilibrium variance of the CCF turns out to be dominated by bulk-like contributions, which induce the scaling behavior We emphasize that Eq. (3.31) quantifies the typical fluctuations of the Casimir pressure K. Those of the actual Casimir force acting on a surface of area A are characterized by A (∆K) 2 1/2 and thus exhibit the expected thermodynamic scaling behavior [29]. The essential scaling behavior expressed in Eq. (3.31) has been confirmed by Monte Carlo simulations of various lattice models (see Refs. [15] and [41]). Compared to the mean values of the CCF [42], the dispersion (∆K) 2 1/2 of the CCF is orders of magnitude larger and non-universal, as noted previously in Ref. [14]. In order to address this issue, it has been argued in Ref. [14] that the static variance is unobservable because any measurement device experiences a force averaged over a finite time interval. The variance of the average forceK is thus reduced by a factor N = t res /t corr , where t res is the temporal resolution of the measurement device and t corr is the correlation time of short-wavelength fluctuations (which provide the dominant contribution to the variance). For typical experimental setups, N is estimated to be of O(10 4 ), which reduces the dispersion ∆K 2 1/2 = N −1/2 ∆K 2 1/2 to a value comparable to the mean force K [14].
In the remaining part of this study, we shall analyze temporal pressure correlations as well as experimentally observable quantities influenced by the fluctuations of the CCF. We find that in these cases the cutoff dependence is mitigated or even disappears, without the need to invoke the temporal resolution of the measurement device and to formulate specific assumptions about it.

C. Dynamic pressure correlations
We first focus on a film with periodic BCs. The associated two-time correlation function of the CCF at a fixed surface (z = z = 0) results from Eq. (3.4) as In the following, we evaluate the dynamic bulk and film pressure correlations based on the formalism developed in Secs. III A 1 and III A 2.

Bulk pressure correlations
In the short-time limit t τ −2 (which includes the critical case τ = 0), terms proportional to τ or τ 2 in Eq. (3.10) are subdominant and can be asymptotically neglected [43]. Accordingly, in this regime, the dynamic bulk pressure correlations at a surface located at z = 0 are obtained by inserting Eq. (2.26) into Eq. (3.14): where is the surface area of the d-dimensional unit sphere and F d (ψ) is a time-independent function which is constant for ψ → 0 and d > 2, while it vanishes exponentially for ψ → ∞. Its explicit form is rather lengthy and thus it is not reported here. A numerical evaluation of H d , which is finite in spatial dimensions d > 2, gives H 3 5.299 × 10 −3 and H 4 5.924 × 10 −4 . Conversely, at long times t τ −2 , using the scaling form given in Eq. (2.25) for Eq. (3.14), the bulk pressure correlations are dominated by the term (1/2)τ 2 C(r , 0, t) 2 . Evaluating the remaining integral over r renders

Periodic BCs
In order to determine the dynamic correlations of the film pressure at a fixed surface (z = 0), we insert Eq. (2.27) into Eq. (3.14). A typical term is given by, e.g., where S(p, k, t) is reported in Eq. (2.29). In summary, the dynamic film pressure correlation at a fixed surface is given by The integrand I(p, k m , k n ) in Eq. (3.38) has the symmetry I(p, k m , k n ) = I(p, −k m , −k n ). While it is not symmetric upon exchanging k m and k n , it can be written in such a form by using the fact that m and n are summed over. (In the following discussion, k m and k n refer to the second and third argument of I, respectively.) For small p and τ = 0, the integrand behaves as The difference between Eq. (3.39) and the corresponding scaling behavior in the bulk [Eq. (3.34)] stems from the continuum spectrum of modes with small k associated with the infinite transverse (z) direction of the bulk system. We finally remark that, in a completely finite volume one expects a long-time relaxation behavior different from Eq. (3.39) due to the isolated zero mode, which is absent for conserved dynamics [9]. At short times (t L 4 , still keeping t τ −2 ), instead, modes with k j t −1/4 contribute significantly to the integral in Eq. (3.38). We first focus on dimensions d ≥ 4, in which case the integrand in Eq. (3.38) is typically finite for p → 0. In order to estimate ∆P Note that the notions of a short-and long-time limit are different in the film and in the bulk (compare Sec. III C 1). For dimensions d < 4, the contribution to the integrand in Eq. (3.38) pertaining to k m = 0, k n = 0 is diverging for p → 0 (as discussed above) and thus has to be analyzed separately in order to properly determine the short-time behavior. The scaling behavior of this contribution is found to be where σ 1 is an arbitrary dimensionless parameter introduced to enable the small-P expansion of the integrand, and R (σ) is a remainder, which is exponentially suppressed and will be omitted henceforth. Since the first equation is independent of σ, the sum in the second equation must be independent of σ, too.
Once this problematic contribution has been removed from Eq. (3.38), the continuum approximation leading to Eq. (3.40) can again be applied, such that, in total, the short-time scaling behavior is given by the sum of Eqs. (3.40) and (3.41). Specifically in d = 3, this renders where the quantity H d is defined in Eq. (3.34). Case t τ −2 . From t τ −2 , it follows that t −1/4 (τ t) −1/2 . The form of the exponential terms in Eq. (3.38) then implies that (τ t) −1/2 provides an upper bound for q below which the integrand contributes. Thus, in this long-time regime, one has q 2 /τ 1/(t 1/2 τ ) 1 [44], such that the square bracket in Eq. (3.38) essentially reduces to τ 2 . We apply this approximation analogously to the exponentials in Eq. (3.38) and distinguish the two cases τ t L 2 and τ t L 2 . In the first case we have τ 2 t max(1, τ L 2 ), which allows us to set k m = k n = 0. Under these assumptions, Eq. (3.38) renders the following long-time behavior of the film pressure correlations: For τ −1 τ t L 2 , on the other hand, we can evaluate Eq. (3.38) by replacing the sum over k j by an integral. As before, this is equivalent to inserting Eq. (2.25) into Eq. (3.14), keeping only the dominant term (1/2)(τ C (p) ) 2 . This gives the following intermediate asymptotic behavior:  As it was the case for periodic BCs, we find that (∆P f ) 2 (t) has a weak logarithmic divergence for d = 3, but is finite for d > 3 and t > 0. Proceeding analogously to the previous subsection, one obtains the following asymptotic  . For the evaluation of ∆P f (t)∆P f (0) at bulk criticality (τ = 0) we used a low momentum cutoff λ 10 −9 /L (this specific value is chosen for illustrative purposes). For times shorter than those covered by the plot, the curves corresponding to L 2 τ = 10 5 (green) and 10 (red) approach the asymptotic critical (τ = 0) short-time behavior given in Eq. (3.43). Due to the smaller value of L 2 τ , the red curve lies closer to this asymptote and crosses over to the off-critical asymptote [Eqs. (3.44) and (3.47b)] later than the green one. Except for the case of the power law ∝ t −3/4 , the asymptotic laws reported in these plots do not involve a fitting parameter. We note that ∆P f (t)∆P f (0) has the same dimension as A where the function F  Note that these asymptotic predictions do not involve a fitting parameter. In order to regularize the logarithmic infrared divergence of ∆P f (t)∆P f (0) in d = 3 at bulk criticality (τ = 0) [see Eqs. (3.43) and (3.47e)], we use a value λ = 10 −9 /L for the infrared cutoff. The value of the prefactor is sufficiently small in order to satisfy the condition of considering a thin film, but otherwise arbitrary and chosen for illustrative purposes. For the parameter σ, which is rather technical and the origin of which is explained in the context of Eq. (3.41), we use a value σ 0.1/L motivated by numerical considerations. We remark that, due to the logarithmic dependence, using values of these parameters which differ even by an order of magnitude do not noticeably affect the quality of the asymptotic approximations.

IV. FLUCTUATIONS OF THE BOUNDARY
In the preceding section, we have studied static and dynamic pressure correlations at a spatially fixed boundary of the film. We now relax this assumption and consider a responsive, i.e., movable film boundary.

A. Dynamical model
We start with a description in terms of the physical timet [see Eq. (2.9)] and denote the position of the film boundary by R(t), which initially is at R 0 ≡ R(t = 0) = L [see Fig. 1]. We consider the motion of the boundary to be overdamped and subject to the instantaneous force K z [Eq. (2.16)] as well as to Gaussian white noise η with covariance η(t)η(t ) = 2δ(t −t ). Accordingly, we propose the following Langevin equation for the time evolution of R: where Γ is a mobility coefficient (which has the same dimension as L 2 /(Tt)) [45]. The prefactor T A in front of K arises because we have defined K as a force per area A of the film at temperature T [see Eq. (2.18)]. A reflecting boundary, representing a wall impenetrable to the moving boundary, is taken to be positioned at R = 0. In the second equation in Eq. (4.1), the CCF K is split according to Eq. (3.3) into a mean, K φ , and a fluctuating part, ∆K, (with respect to the OP field). Since we consider thermal equilibrium for the OP field, the mean force K does not explicitly depend on time. The white noise η accounts for the molecular momentum exchange between the solvent and the surface, which would be present even without a coupling to the OP. Upon introducing the rescaled time t = γt as defined in Eq. (2.11), Eq. (4.1) turns into where the parameter D is the bare diffusivity with dimension [D] = [L −2 ]. The fact, that the stochastic process ∆K(R,t)) has a nonlinear dependence on R, raises the issue of the proper stochastic calculus (e.g., Ito or Stratonovich) to be used in Eq. (4.1) [28,46]. In principle, this can be addressed by deriving Eq. (4.1) from a more basic model, e.g., by adiabatically eliminating the OP degrees of freedom from a coupled system of Markovian Langevin equations for the boundary and the fluid medium [22-24, 28, 47-49]. While this is beyond the scope of the present study, it turns out that, in the Markovian limit, the variance of ∆K is approximately independent of R, rendering the choice of the stochastic calculus to be immaterial for the present analysis (see Sec. IV C 2). Equation (4.2) describes the random motion of R in the Casimir potential associated with K φ , subject to the Brownian noise η [see Eq. (4.6) below] and an additional non-Markovian noise ∆K. While Eq. (4.2) resembles previously proposed Langevin equations for inclusions in critical media [22][23][24]50], here we assume that the presence of a movable boundary does not modify the equation of motion of the OP [Eq. (2.12)]. In fact, the only effect of the boundary on the OP is to impose BCs [see Eq. (2.8)]. Accordingly, the noise provided by ∆K is not balanced by a corresponding friction term, implying that detailed balance is not satisfied by the model defined by Eqs. (2.12) and (4.2). This is also reflected by the fact that the distribution of the boundary position does not approach the expected equilibrium form at long times (see Appendix A). We thus refer to our model as having a "passive boundary". For a more adequate description of the dynamics of a colloidal particle in a critical fluid in thermal equilibrium [5,8,21,26], detailed balance has to be re-established [51]. This can be accounted for by either adding appropriate coupling terms to the OP dynamics in Eq. (2.12) (see Ref. [50]), or by suitably modifying Eq. (4.2) in order to satisfy the fluctuation-dissipation theorem. We will return to the latter approach in Sec. IV C 2 below. However, compared with the models in Refs. [22][23][24]50], a distinctive feature of the present model is that it allows one to enforce strict no-flux BCs at the location of a boundary.
The mean equilibrium CCF (per area A and temperature T ) has the form where Ξ is a universal scaling function depending on the boundary conditions. Explicit expressions within the Gaussian approximation are reported in, e.g., Refs. [42,52]. At criticality (i.e., ξ → ∞), one has Ξ(0) < 0 for periodic and Neumann BCs [see Eq. (3.32)]. At short distances, the electrostatic repulsion between the fixed wall and the mobile boundary must be taken into account [21,53], thereby regularizing also the singularity of the CCF [Eq. (4. 3)] appearing within the continuum description in the limit R → 0. Phenomenologically, this electrostatic repulsion can be described by the potential (defined, analogously to the CCF, per area A and temperature T ) where α is an effective parameter (of the same dimension as L) and D denotes the Debye screening length. Introducing formally the Casimir potential V C (z) associated with the mean CCF via K(z) φ = −dV C (z)/dz, we define an effective potential which gives rise to a force −U (z) replacing K(z) φ in Eq. (4.2). Accordingly, Eq. (4.2) is replaced by the Langevin equation The effective potential U(z) has a minimum at a finite distance z 0 from the wall. In order to facilitate an analytical treatment, we consider its quadratic approximation around its minimum: An explicit expression for the (temperature dependent) parameter κ (which has the same dimension as L −d−1 ) can in principle be obtained from Eqs. (4.3) and (4.4), but is not necessary in the following. As a further simplification, we focus on the limit in which the dynamics of the field φ is faster than that of the position R (adiabatic approximation).
Since within the Gaussian model B [33] the relaxation time of a critical OP fluctuation grows ∝ L 4 , where L is the characteristic system size, the adiabatic approximation requires that both the film and the "bulk" part of the system are finite (see Fig. 1). Still, we take the bulk to be sufficiently large so that a continuum approximation is appropriate. Requiring L to be a measure of the largest length scale in the system (which implies also L A 1/(d−1) ) allows us to introduce a dimensionless "adiabaticity" parameter In terms of χ, Eq. (4.2) can be expressed as (4.9) Thus the adiabatic limit requires χ 1 and facilitates a perturbative solution of Eq. (4.9). For typical experimental conditions of a colloid immersed in a critical solvent, the condition χ 1 is fulfilled [26].

B. Boundary located in the bulk
Here, we assume here that R(t) is sufficiently far from the wall such that the effective force U in Eq. (4.6) can be neglected. Accordingly, we have a bulk-like system on both sides of the boundary. Integrating Eq. (4.14) The factor 2 in Eq. (4.12) arises because the bulk pressure acts on both sides of the boundary, while in Eq. (4.13) the prefactor A is introduced in order to render M independent of A. We note that, due to the presence of ∆K, the process described by Eq. (4.2) is highly nonlinear. In order to evaluate Eq. (4.13), we note that in Eq. (4.9) the OP-induced noise ∆K is of higher order in χ than the white noise η. Accordingly, solving Eq. (4.9) (assuming, as before, U = 0) in the regime χ 1 to first order in χ yields a free Brownian motion for R(t): which is to be used in Eq. (4.14) [see Sec. IV D below]. The required two-time joint probability distribution for a free Brownian process is given by [28] p(R s , R s , R 0 , s, s ) = 1 (4.16) which renders the following variants of the characteristic function: e iQ(R(s)+R(s )) = e 2iQR0−DQ 2 (s+s +2 min(s,s )) . (4.17b)

C. Boundary located near a wall
Here, we retain the effective force U in Eqs. (4.6) and (4.9). Within the regime χ 1, the fluctuating contribution of the CCF ∆K is subdominant compared to the white noise η, while the deterministic term related to U has to be regarded as formally to be of the same order in χ as η. (This assignment is also applied, e.g., upon solving a standard Langevin equation of a Brownian particle [28].) Accordingly, for χ 1, the leading order solution R(t) of Eq. (4.9) [based on the assumption formulated in Eq. (4.7)] is an Ornstein-Uhlenbeck process [28]. However, the calculation of an average as in Eq. (4.13) over this process leads to intractable expressions. We therefore solve Eq. (4.9) instead by applying two complementary approaches: in one (see Sec. IV C 1), we formally integrate Eq. (4.9) in time and focus on short times, in which case the Ornstein-Uhlenbeck process reduces to a free Brownian motion. In a second approach (see Sec. IV C 2), we determine a Markovian approximation of ∆K in the adiabatic limit χ 1. Although the second approach holds at all times, we use it specifically to study the long-time behavior. where M b is given by Eq. (4.14) and Note that, in order to to simplify the calculation, we consider the film boundary to be initially located at the potential minimum z 0 . Upon using Eqs.

Markovian approximation of ∆K
Here, we seek an alternative solution of Eq. (4.9) which is valid for all times. In order to identify a suitable approximation, we define a new time variablet via t =t/χ (4.24) and introduce the quantitiesR(t ) = R(t/χ),η(t ) = 1 √ χ η(t/χ) [with η(t )η(t ) = 2δ(t −t )]. In terms of them, Eq. (4.9) becomes  In the limit χ 1, one expects from Eqs. (4.25) and (4.26) that ∆K decorrelates rapidly in time, justifying a Markovian approximation forN . Furthermore, we assume ∆K to be Gaussian, such that cumulants higher than the second one reported in Eq. (4.26) are zero [54]. From standard adiabatic elimination techniques [28,55] it follows that ∆K can be approximated by a Markovian noise ∆K M with the variance The resulting Stratonovich-type Langevin equation [28,55] is given by and corresponds to the following Fokker-Planck equation for the probability density P (R, t) of the position R of the movable boundary: where the term involving N (R) is a consequence of the Stratonovich stochastic calculus, also referred to as "spurious" drift [28]. The expression for the noise amplitude N (R) (which has the same dimension as L 4−2d ) is obtained by using Eqs.
where the wavenumbers used in the sum are defined in Eq. (2.28) [see also the comment after Eq. (3.46)]; K d is defined after Eq. (3.38), and n + 2p 4 ] k 2 m 2p 2 1 + e 2iRkm 1 + e 2iRkn + k 2 n 3e 2iR(km+kn) + e 2iRkm + e 2iRkn + 3 − 2k m k n 2k 2 n + p 2 1 + e 2iR(km+kn) + 1 + e 2iRkm 2p 2 k 2 n + 2k 4 n + p 4 1 + e 2iRkn , (4.31b) for τ = 0 [56]. The expressions of these quantities for τ = 0 are rather lengthy and thus they are not stated here. In general, the integral over p diverges both at its lower (0) and upper (∞) limits and thus has to be regularized, which will be discussed below. Symmetry arguments show that m,n Q (N) is real valued. Note that Q (N) is also a function of R, while Q b and Q (p) are not. However, a numerical analysis reveals that Q (N) is only weakly depending on R (provided Lλ 1), which allows us to make the approximation Q (N) (R) Q (N) (R = 0) so that N (N) (R) N (N) (R = 0). This renders N to be independent of R for all BCs considered here. Accordingly, Eq. (4.29) can be approximated by Next, we determine the Fokker-Planck equation (FPE) of a boundary in a system with detailed balance. To this end, we assume that the corresponding FPE takes the same form as in Eq. (4.32), up to a different (but still spatially constant) mobility Γ eff , i.e., (4.33) We now require that the fluctuations of the boundary are described by the equilibrium distribution (see Appendix A) In order that the steady-state solution of Eq. (4.33) agrees with P eq , the following relation must hold [28]: This can be regarded as a consequence of the fluctuation-dissipation theorem [57]. We shall refer to a boundary described by the model in Eqs. (4.33) and (4.35) as "thermalized". Importantly, asymptotically at short times, the behaviors of the MSD resulting from Eqs. (4.32) and (4.33) are identical and independent of the mobility (see Sec. IV D). The derivation of the effective mobility of a boundary from a more rigorous approach (following, e.g., Refs. [50,58,59]) is left for future studies.
Denoting the second term in the square bracket in Eq. (4.30) by m b (where b stands for bulk), power counting (taking into account the integrals over k m and k n ) implies the asymptotic behaviors (4.36c) The first term in the square brackets in Eq.
independently of τ . In order to regularize the integral over p in Eq. (4.30), we replace its lower and upper cutoff by wavenumbers λ and Λ, respectively. Accordingly, within the Markovian approximation in d = 3, with the characteristic (lateral) system size L if λ 1/L is taken for the lower integration limit.

D. Short-time behavior of the MSD
Based on the expressions derived above, in the following we discuss the evolution of the mean-squared displacement (MSD) of the movable boundary at short times and elucidate how the fluctuations of the CCF affect it.

Preliminaries
The expressions of M f [Eq. (4.21)] for periodic BCs and M b [Eq. (4.14)] for the bulk are structurally similar and can thus be analyzed together. M (p) f is obtained by inserting the Fourier transform of the correlation function given in Eq. (2.27) into Eq. (3.13). In contrast to the calculation of the pressure variance at a fixed position discussed in Sec. III C 2, we now have to retain the dependences on the z-coordinate (as well as on the second time argument t , which was previously set to 0) and subsequently average over the distribution of the boundary position z = R. Moreover, we cannot use the simplifications stemming from Eqs. (2.35) and (3.15). In order to keep the expressions tractable, we focus on the bulk critical point (τ = 0).
We illustrate the calculation of M f by considering a typical contribution to the integrand in Eq.  where the function S (which depends only on even powers of p) is given in Eq. (2.29). In order to proceed, we approximate R(s) as a Brownian motion as formulated in Eqs. (4.15) and (4.22). We recall that, in the bulk, this approximation holds at all times, whereas it requires t 1/(DAκ) in the case of a film. Next, we replace in Eq. (4.39), as well as in all other contributions stemming from Eq. (4.21), the characteristic function by the expression in Eq. (4.17a) (setting Q = k m + k n ). After performing the time integrals, one obtains in total for a film with periodic BCs , (4.40b) and where K d is defined in the context of Eq. (3.38). The corresponding expression for M b is given by replacing in Eq. (4.40) the sums over k m,n by integrals according to which holds for any function f (k). For a film with Neumann BCs, M (N) f follows from using Eq. (2.32b) in Eq. (3.10) and then by proceeding analogously to the calculation which leads to Eq. (4.40). We note that here both variants of the characteristic function in Eq. (4.17) are needed and that the resulting expression has to be evaluated with 2L instead of L. Since the result is rather lengthy it is not reported here. With the exception of the short-and long-time limits (see below), M f,b have to be determined numerically.

Bulk system
In order to calculate Eq. (4.40) for a bulk system, in the integrals we rescale all momenta by t −1/4 and regularize the integral over p by replacing its lower and upper integration limits by λ and Λ, respectively. The resulting expression is independent of D in the limit t D −2 , which we identify as the short-time regime for the bulk. In this regime one finds the scaling behavior (up to possible logarithmic corrections, see below)

Film with periodic BCs
We now consider a film boundary fluctuating near a fixed wall and note that, for t L 4 , the exponential in Eq. (4.40b) varies mildly over a wide range of momenta. Accordingly, one can approximate the sums in M (p) f in Eq. (4.40) by integrals and, as above, rescale all momenta by t −1/4 . As before, we additionally consider the regime t 1/(DAκ) and note that, consequently, t D −2 is automatically fulfilled owing to the adiabatic approximation

Film with Neumann BCs
A scaling analysis of the expression M

E. Long-time behavior of the MSD
We now turn to a discussion of the MSD of a movable boundary and of the influence of the CCF on it at long times.

Bulk system
The long-time behavior of a boundary in a bulk system is provided by Eq. (4.12) together with Eq. (4.40) [after making the appropriate replacement indicated in Eq. (4.41)]. According to the analysis leading to Eq. (4.42), a longtime regime emerges for times t D −2 . In this limit, the exponential term in Eq. (4.40b) turns out to be negligible, and the remaining expression renders a linear time dependence: dp v(p). (4.48) In order to assess the relevance of the cutoffs λ and Λ, we note that the integrand behaves as v(p → ∞) ∝ p d−4 , while the behavior at small p depends in general on D. For D → 0, one has v ∝ p d−4 , while at nonzero D, a numerical analysis indicates a behavior close to v ∝ p d−3 . In spatial dimension d = 3 with nonzero D, one may thus set λ = 0, while retaining a logarithmic dependence on the microscopic length scale ε ∝ 1/Λ via the dimensionless combination Dε 2 . Since the function v depends on the diffusivity D, in general the calculation of Eq. (4.48) requires a numerical approach (see Sec. IV F for further discussions).

Film with periodic BCs
In order to determine the long-time behavior of the MSD of a fluctuating boundary in proximity to another wall, we invoke the Markovian approximation discussed in Sec. IV C 2, which leads to an Ornstein-Uhlenbeck process for the boundary position [see Eqs. (4.28), (4.32), and (4.33)]. Equations (4.32), (4.33), and (4.35) imply at long times a static Gaussian distribution with variance [28] [ We remark that, off criticality, the term proportional to A/L 2 in the square brackets disappears owing to the much weaker divergences of the expressions in Eqs. (4.36) and (4.37) for τ > 0. Equation (4.49b) provides a first approximation for the variance of the position of a colloidal particle close to a wall and immersed in a nearly critical solvent in thermal equilibrium, consistent with previous studies [5,8,26]. Experimentally, the typical magnitude of the fluctuations is of the order of [∆R] 2 1/2 ∞ ≈ 10 nm [5,60]. In order to evaluate Eq. (4.51), which is experimentally less relevant due to the underlying passive boundary approximation, we take A to be the cross sectional area of a typical colloidal particle, A (1 µm) 2 , and estimate the molecular length scale of the solvent as ε 1 nm, in which case Eq. (4.51) reduces to ∆R 2 pass ∞ 1 κA (1 + 0.75 × DA). Because Eq. (4.8) (with L A 1/2 ) implies DA 1 in the adiabatic limit, the contribution to [∆R] 2 pass ∞ stemming from the fluctuations of the CCF is subdominant relative to those stemming from the thermal white noise in Eq. (4.28).

Film with Neumann BCs
If the fluctuating film boundary is coupled to the OP via Neumann BCs, one obtains at long times the same expression for the variance ∆R 2 ∞ as in Eq. (4.49), but with N based on Eq. (4.31b) (with R = 0). Specifically, for d = 3 and τ = 0 we obtain, by using the same estimates as in Eq. (4.51), the long-time variance  At short times, the MSD behaves identically (up to constant prefactors) in the bulk and in proximity to another wall. This is physically expected, because the influence of the wall requires a certain amount of time to reach the moving boundary. The data in Fig. 3(a) also cover the long-time regime of the bulk contribution M b , which, according to Eq. The MSD depends only weakly on these differences, because in both cases it is dominated, within the adiabatic approximation, by the Brownian diffusivity D.

V. SUMMARY, CONCLUSIONS, AND OUTLOOK
In this study, we have analyzed the fluctuations of the critical Casimir force (CCF) acting on the boundary of a film with periodic, Neumann, or Dirichlet BCs (Fig. 1). Our predictions are obtained within a time-dependent Ginzburg-Landau description of an equilibrium fluid system with a conserved order parameter (model B [33]) and based on a rigorous definition of the instantaneous CCF [9,[30][31][32]. We have considered a film consisting of a spatially fixed boundary and a parallel boundary which is either fixed (see Secs. III B and III C) or movable (see Sec. IV). The movable boundary is coupled "passively" to the OP (see Sec. IV A) and is subject to the action of thermal white noise and the fluctuating CCF. Due to the nature of the passive coupling, the boundary acts on the OP only by imposing BCs, but not by exerting additional friction. By utilizing the fluctuation-dissipation relation, we have shown how this approximation can be improved in order to describe a boundary in thermal equilibrium (see Sec. IV C 2). We have distinguished between the motion in a bulk-like environment and the motion close to that boundary of the film which is kept fixed. In this latter case, an effective confinement potential arises due to the attractive mean part of the CCF and a short-ranged electrostatic repulsion. This situation is analogous to a CCF-induced trapping of a colloid close to a solid wall, i.e., the inner boundary [5,25,26].
Our results can be summarized as follows: the variance acquires contributions from the fluctuations of the CCF [see Eq. (4.49a)]. By contrast, these contributions are absent in the case of a boundary in thermal equilibrium [see Eq. (4.49b)]; the variance of its position is solely determined by the critical Casimir potential [5,21,26].
The formally divergent equilibrium variance of the CCF [see Eq. (3.31)] can be considered as an artifact of the corresponding continuum field theory, which cannot be observed directly. In fact, in previous studies of (critical) Casimir forces [14,27], it has been shown that this divergence problem is mitigated or even absent if time-averaged quantities are considered. In the present study, we have instead directly linked the statistical properties of the CCF to experimentally observable quantities, such as the position of a film boundary. Our analysis renders a position variance which is finite and reduces to the variance of a Brownian particle in an effective potential set by the mean CCF. This finding is consistent with previous experiments on CCFs in wetting films [17][18][19][20] or in colloidal systems [5,21]. Heuristically, this situation is analogous to the fact that observable quantities obtained from a standard Langevin equation are usually finite, despite the random force being correlated ∝ δ(t − t ), which formally diverges in the static limit t → t . In the present case, the force correlations diverge algebraically (see Fig. 2) but are still integrable in the sense of the MSD (see the discussion in Sec. IV D 5). Moreover, due to the limited temporal resolution of a measurement device, experimental measurements are typically not sensitive to fluctuations of the CCF at frequency scales pertaining to the relaxation of the short-wavelength modes [14].
Our analysis of the dynamics of the boundary is based on two central assumptions: first, the moving boundary acts on the OP only via the imposed BCs, not via coupling terms in the equation of motion (passive boundary approximation) [22,23,50]. While, within this approximation, detailed balance is broken and the system is inherently out of equilibrium, we have also shown how to re-establish detailed balance at the level of the Fokker-Planck equation for the position of the boundary [see Eq. (4. 33)]. Second, the relaxation dynamics of the boundary is much slower than the one of the order parameter; this is fulfilled in typical experiments [26] and also ensured in our study because we have assumed a finite, albeit macroscopically large, system size. Within the passive boundary approximation, the actual Casimir potential-which is related only to the mean CCF-cannot be simply inferred by assuming a standard Boltzmann distribution of the position [21], because OP fluctuations render additional contributions to the variance which can vary in space [see Eq. (4.30)]. Asymptotically, in the adiabatic limit these contributions are negligible.
Concerning future research, a systematic derivation of the dynamical equation from a more fundamental model, beyond the passive boundary approximation, would be desirable. It would also be interesting to extend the present analysis towards the dynamics of a spherical colloidal particle which has a finite radius, is immersed in a critical solvent, and floats near a wall [26,61]. It is justified to expect a wealth of additional effects due to hydrodynamic interactions, the non-planar geometry, and the strong preferential adsorption of the OP at the surfaces [7,62,63].

ACKNOWLEDGMENTS
AG acknowledges financial support from MIUR PRIN project "Coarse-grained description for non-equilibrium systems and transport phenomena (CO-NEST)" n. 201798CZL.

Appendix A: Equilibrium distribution of a movable boundary
Here, we discuss the probability distribution of the position of a movable boundary in thermal equilibrium with a critical fluid in a half-space [30,32]. The Hamiltonian of this system is given by F(R; The equilibrium probability distribution of R follows by integrating over φ: The mean equilibrium CCF K (per area A and temperature T ) acting on the boundary [which in Eq. (2.21) is expressed in terms of the averaged stress tensor] can be equivalently obtained from Z(R) [30]: