Critical relaxation and the combined effects of spatial and temporal boundaries

We revisit here the problem of the collective non-equilibrium dynamics of a classical statistical system at a critical point and in the presence of surfaces. The effects of breaking separately space- and time-translational invariance are well understood, hence we focus here on the emergence of a non-trivial interplay between them. For this purpose, we consider a semi-infinite model with $O(n)$-symmetry and purely dissipative dynamics which is prepared in a disordered state and then suddenly quenched to its critical temperature. We determine the short-distance behaviour of its response function within a perturbative approach which does not rely on any a priori assumption on the scaling form of this quantity.


I. INTRODUCTION
It is a trivial observation that any physical system has actually a finite extent; as a consequence, descriptions which assume translational invariance can at most capture its bulk features, with surface effects representing subleading corrections which decay upon moving away from the boundaries. This decay is controlled by the presence of an inherent length scale, which sets the "range" of these surface effects. Rather generically, such a scale corresponds to the correlation length ξ of the system, which encodes the separation beyond which different regions of the extended system are no longer statistically correlated. Accordingly, one can identify three distinct instances: (i) ξ is comparable with the typical linear extent L of the sample and the behaviour of the system is strongly affected by its finiteness, as every point effectively feels the presence of the boundaries. In this case finite-size effects emerge [1] and the thermodynamic quantities explicitly depend on L.
(ii) ξ is significantly smaller than L but one focuses on the behaviour at spatial points located at a distance d ≫ ξ from the boundaries, in such a way that the effects of the boundaries can be neglected and the system can be modeled as being infinitely extended, with translationally invariant properties.
(iii) ξ is significantly smaller than L but one considers the behaviour of the system within a distance d ≪ ξ from the boundaries. If ξ is smaller than the curvature of the latter and there are no wedges or tips, a suitable description is provided by semi-infinite models with a flat boundary and a lingering (approximate) translational invariance in all spatial directions parallel to the surface.
In the following, we focus on case (iii). The differences between the behaviour close to the boundaries and the one far from them, i.e., in the bulk, formally arise due to the explicit breaking of the translational invariance along the direction orthogonal to the surface, which allows an extended freedom in the system: for example, a two-point correlation function C(x, y) is no longer constrained to be a function of the distance |x − y|.
A faithful description of these systems generically requires accounting for all the microscopic features which characterize both the bulk and the surface; in turn, this implies that the corresponding behaviour is highly system-specific. However, when the correlation length ξ becomes large with respect to microscopic scales, collective behaviours emerge which are not determined by the underlying microscopic structure, but by those coarse-grained properties which do not really depend on the considered scale, such as symmetries, range of the interactions and the (effective) spatial dimensionality. These circumstances, which represent a hallmark of systems undergoing a continuous phase transition, lead to universality. In other words, several relevant quantities can be identified which, due to the very fact that the microscopic details become inconsequential for their determination, take the same values in many different systems, which in turn share the same gross features and constitute as a whole the so-called universality class of the transition. Consequently, it is sufficient to study just one representative system in order to gain information on the whole class it belongs to. Continuous phase transitions are typically associated with the spontaneous breaking of some underlying symmetry [2], which is highlighted by the behaviour of the so-called order parameter ϕ (e.g., the local magnetisation for an Ising ferromagnet) upon crossing the critical point.
The emergence of universality is currently understood within the framework of the renormalization group (RG) [2][3][4]. Its transformations act by enlarging the length scale at which a system is described, progressively blurring details at shorter scales. Under the assumption that ξ represents the only inherent non-microscopic length scale, one must conclude that its divergence deprives the system of any typical scale and therefore its physical behaviour must become self-similar under scale dilatations and the RG transformations reach a fixed point. The different physical systems which end up falling into it constitute a certain universality class. As a matter of fact, as long as the interest lies in the study and determination of universal quantities, the resulting mesoscopic description can be formulated in terms of fields on a space-time continuum, which makes it possible to use standard field-theoretical methods in order to calculate many relevant quantities [2].
The discussion below complements the results presented in reference [5] by providing in detail the analytical derivation of the first corrections to the linear response function of a system bounded by a flat surface and subject to a temperature quench. In particular, in section II we briefly discuss the main features emerging in critical systems with spatial and temporal boundaries, while in section III we set up the description of the aforementioned model in terms of a suitable (field) theory on the continuum; by using standard field-theoretical methods and a RG-improved perturbation theory we explicitly calculate the relevant universal quantities and show the emergence of an unexpected edge behavior. Finally, in section V we draw our conclusions.

II. EQUILIBRIUM TRANSITIONS AT SURFACES AND NON-EQUILIBRIUM CRITICAL DYNAMICS AFTER A QUENCH
A. Spatial boundaries: Equilibrium critical behaviour at surfaces Being originally devised for describing the behaviour of unbounded and uniform systems, which provide a good approximation of regime (ii) above, the RG has been subsequently generalized in order to describe finite-size effects (regime (i)) [1] and the presence of boundaries [6][7][8][9] (regime (iii)). In the latter cases, the breaking of translational invariance plays a fundamental role, leading to the appearance of novel scaling behaviours which can be observed, for example, in the correlation functions of the order parameter ϕ upon approaching the boundaries [6][7][8][9][10] and are characterized by new universal exponents which cannot be inferred from the bulk ones. Within the RG approach, the additional parameters which describe the gross features of a boundary (e.g., a different interaction strength at the surface with respect to the bulk) typically give rise to a splitting of the original bulk universality class into a set of surface subclasses characterized by these novel (surface) exponents. A number of analytical [8,11,12], numerical [13][14][15] and experimental (see, e.g., reference [16]) studies have been carried out to investigate semi-infinite and film geometries, whereas wedges, edges [17,18], as well as curved and irregular surfaces [7,9] have been studied to a lesser extent.
In order to exemplify some of the features mentioned above, we focus now on the Ising universality class, which is effectively described by an effective Ginzburg-Landau free-energy where r ∝ T − T c controls the distance (in temperature) from the critical point and g > 0 sets the strength of the interaction. In the presence of a surface, the spatial integration above is restricted to the half-space x ⊥ > 0, x ⊥ being the coordinate in the direction orthogonal to the surface.
Furthermore, a term F 1 which accounts for the surface properties has to be added to equation (1). In view of the eventual application of RG, one can conveniently focus on the most relevant surface term which is allowed by dimensional analysis: it turns out to have the form where the integration only runs on the coordinates parallel to the surface, and effectively encodes the boundary condition [7] The surface enhancement c 0 accounts for the differences in the interaction strength between the bulk and the surface. Its presence modifies the characteristic phase diagram of the Ising model, as shown in figure 1. In addition to the usual ferromagnetic (SO/BO in the figure) and paramagnetic (SD/BD) phases, one may find a third one (SO/BD) in which the surface (S) is ordered (O) whereas the bulk (B) is disordered (D). Accordingly, the number of different transitions (and therefore universality classes) increases to four and they are referred to as ordinary (the surface orders with the bulk), surface (only the surface becomes ferromagnetic), extraordinary (the bulk orders in the presence of an already magnetized surface) and special (the point at which the critical lines coalesce) [6,7]. Note that every transition except for the ordinary one inherently requires the possibility for the surface to order independently from the bulk, i.e., that its dimensionality d − 1 is strictly larger than the lower critical dimension d lc = 1 of this system. Thus, semi-infinite ferromagnets in d = 2 can only undergo the ordinary transition. We recall that for continuous symmetries (e.g., O(n) with n ≥ 2) one has instead d lc ≥ 2 [2]. The ordinary and special transitions take place with a vanishing order parameter both in the bulk and at the surface; thanks to this, they admit a unified description which differs only by the effective boundary conditions (see equation (3)) cast onto the order parameter, which within the Gaussian approximation are of Dirichlet (i.e., c 0,ord = +∞) and Neumann (c 0,sp = 0) type, respectively. As stated above, the magnetization m s at the surface and the one m b in the bulk show different behaviours when varying the temperature T in the critical regime, i.e., which require the introduction of a new critical exponent β 1 ; this exponent is reflected also in the generic dependence of m on x ⊥ upon approaching the surface, which is determined by a so-called short-distance expansion (SDE) [7] m( Thus, one can characterize the critical surface behaviour by means of the algebraic dependence of certain quantities on the distance from the surface.

B. Temporal boundaries: short-time critical dynamics after a temperature quench
Rather surprisingly, the same formalism can be applied in a dynamic framework as well, after performing a sharp variation of one of the thermodynamic control parameters of the system (e.g., the temperature) [19]. Intuitively, the instant t 0 at which such a quench is performed separates the equilibrium regime (t < t 0 ) from the non-equilibrium one (t > t 0 ) and therefore acts qualitatively exactly in the same way as a spatial surface which stands between the outside and the inside of a system. In this case the distance x ⊥ from the surface is given by the time t − t 0 elapsed from the quench.
The emerging scaling properties of the dynamics of a statistical system near criticality are more easily discussed in terms of the evolution of the associated order parameter field ϕ on the continuum. For a broad class of systems (the so-called model A universality class in the notion of reference [23], see further below) the features of such an evolution can be effectively captured via a Langevin equation of the form where F is the effective free-energy, Ω a kinetic coefficient, while η is a Gaussian white noise with which accounts for the fluctuations due to a thermal bath at temperature T and represents an external source of dissipation (for simplicity, units are chosen such that k B T = 1). These kind of equations can be mapped onto a field-theoretical description [20][21][22] via the introduction of the so-called response field ϕ. The corresponding action reads The newly-introduced variable ϕ actually encodes the response properties of the system: in fact, to an external perturbation applied at time s which couples linearly to ϕ. In fact, introducing a source Analogously to the spatial case, the quench limits the integration to times t ≥ t 0 . Moreover, the determination of the dynamics needs an initial condition ϕ 0 for the field ϕ to be specified, for example via its probability distribution P(ϕ 0 ), which can be conveniently written in the exponential form The associated boundary action S 0 is integrated only over spatial coordinates and is thus akin to F 1 in the spatial case illustrated in section II A. By an RG argument, one can analogously account only for the most relevant terms in S 0 ; considering again the Ising universality class as an example, one would then have With the additional assumption that the initial state is very far from criticality, one can also neglect the second addend and thus obtain a Gaussian initial condition with variance τ −1 0 . Although the structure of S 0 looks the same as the one of F 1 in equation (2), one has to take into account the fact that P(ϕ 0 ) is a probability and thus τ 0 can be neither vanishing nor negative. Hence, without explicitly breaking the Ising symmetry ϕ 0 ↔ −ϕ 0 of the action, the only fixed (i.e., critical) point one can identify is the equivalent of the ordinary one with τ 0 → +∞ [19], corresponding to a state with vanishing correlations, typically related to a very high temperature.
Depending on the gross features of the dynamics, such as conservation laws etc., a splitting of the equilibrium universality classes is found [23], associated to the appearance of new universal quantities, such as the dynamical exponent z which encodes the difference in scaling dimension between space and time coordinates [23] and thus describes how the typical linear relaxation time t R ∼ ξ z grows upon approaching the critical point. For example, a classical Ising model on a lattice which evolves in time via thermally-activated independent flips of the spin (Glauber dynamics) belongs to the so-called model A universality class (in the notion of reference [23]) and is characterized by having, within the Gaussian approximation, z = 2. Conversely, if the evolution conserves the total magnetization (Kawasaki dynamics), i.e., it makes domains diffuse or split, or combine, the universality class changes to model B, with Gaussian dynamical exponent z = 4, whilst all the universal features of time-independent quantities remain the same. The universal behaviour emerging at a "temporal" boundary takes the form of an initial slip and emerges as discussed above for spatial boundaries, i.e., because of the different scaling behaviour of observables in the initial and final stages of the non-equilibrium dynamics. In fact, analogously to equation (5), this behaviour can be extracted from the short-distance expansion observed for t → t 0 [19].

III. QUENCHING A DYNAMICAL MODEL WITH A SURFACE: EMERGENCE OF EDGE EFFECTS
We consider here a purely dissipative dynamics in the presence of a spatial surface, which can be described by model A reported in equation (6) with the effective Ginzburg-Landau free energy F replaced byF = F + F 1 , which includes the semi-infinite bulk and the surface contributions F (equation (1)) and F 1 (equation (2)), respectively. For the purposes of the present analysis, the kinetic coefficient Ω can be generically set to 1 by rescaling time and noise as t → t/Ω and η → Ωη, respectively. In order to introduce the temporal boundary we consider the system to be prepared in a completely disordered state for t < t 0 = 0, which corresponds to vanishing ϕ and, in turn, to Dirichlet boundary conditions for the dynamics [10,19]. For the spatial boundary, instead, we will consider both the ordinary and special points.
The separate effects on the system of either boundary are well-understood [7,8,11,12,19,[24][25][26]; on the other hand, their interplay has been less extensively studied [5,10,27,28]. Moreover, in reference [27] a power-counting argument has been used for arguing that no new algebraic behaviours arise which are specific to the intersection of the two boundaries (hereafter referred to as the edge) and that all the observed effects are therefore a combination of those independently generated by the surface and the quench; almost all the subsequent studies formulated scaling ansatzes based on this assumption. The analysis of reference [5], instead, is not based on it but proceeds to a direct calculation of the effects of the edge. As we detail below, this leads to the emergence of new (fieldtheoretical) divergences which are sharply localised at the edge and therefore highlight non-trivial modifications to the scaling laws of observables in its proximity. In particular, we present here the calculations of the first-order corrections to the two-point functions of the theory, generalizing it to a O(n)-symmetric model, i.e., one in which the order parameter ϕ is a n-vector whose components ϕ i satisfy Langevin equations of the form (6). For the case n = 1, the resulting predictions for the emerging scaling behaviour were in fact confirmed by Monte Carlo simulations, as briefly reported in reference [5]. Since we have assumed -motivated by the central limit theorem -the noise η to be Gaussian, we can integrate over it in equation (8), thus obtaining an effective action which depends explicitly only on ϕ and ϕ. In the case at hand (see equations (1) and (2)), this corresponds to for the bulk action and for the surface one.

A. Dynamical response and correlation functions
Exploiting the translational invariance along the surface, one can study the spatial Fourier transform of the two-point correlation and response functions C( k ; x, t; y, s) = ϕ( k ; x, t)ϕ(− k ; y, s) and R( k ; x, t; y, s) = ϕ( k ; x, t) ϕ(− k ; y, s) , respectively, which depend on the wavevector k (in the following denoted just by k for simplicity) and on the distances x and y of the two points from the surface, in addition to the times t and s elapsed since the quench. The remaining two-point function ϕ ϕ vanishes identically [29]. For every value of k , the ordinary and special transitions correspond to the boundary conditions {ϕ(x ⊥ = 0, t) = 0, ϕ(x ⊥ = 0, t) = 0} and {∂ x ⊥ ϕ(x ⊥ = 0, t) = 0, ∂ x ⊥ ϕ(x ⊥ = 0, t) = 0} (the latter being valid only within the Gaussian approximation) which are supplemented by the initial condition ϕ(x ⊥ , t = 0) = 0. The resulting correlation and response function within the Gaussian approximation [ (0) ] turn out to be [10][11][12]19]: where the upper and lower signs refer to the special and ordinary phase transitions, respectively, and f ± (α) = (e α ±e −α )/2; C (b,eq) are the corresponding functions in the bulk at equilibrium, i.e., in the absence of boundaries, which are given by In order to highlight the effect of the boundaries on the collective dynamics, hereafter we focus on the case in which the system is quenched right at its critical point, i.e., we fix r = r c for t > 0, where r c is the critical value of the parameter r; beyond the Gaussian approximation r c still vanishes if the analysis is done by using dimensional regularization for the calculation of the relevant integrals. Note that C (0) (b,eq) in equation (17) is constructed from R (0) (b,eq) via the classical fluctuation-dissipation theorem, which holds in equilibrium [30].

B. First-order corrections
The presence of the interaction ∝ g ϕϕ 3 in equation (12) can be accounted for in perturbation theory, which gives rise to corrections to R (0) and C (0) , represented by Feynman diagrams (see, e.g., reference [2]). In the following, we focus on the first non-vanishing correction R (1) to the response function, since the one for the correlation function yields the same critical exponents; this term is represented by the diagram in figure 2, where directed lines correspond to R (0) in equation (15) and undirected ones to C (0) in (16). The corresponding expression is where k q k y, s x, t z, τ with 2z, 2τ ). The upper and lower signs distinguish the special from the ordinary transition. Within dimensional regularisation [2,7] one finds B 0 = 0 and where γ (α, w) = w 0 dz z α−1 e −z is the incomplete gamma function. Each B i (z, τ ), once integrated as indicated on the r.h.s. of equation (18) yields a contribution R i , in terms of which In the following, we set for simplicity k = 0 and t > s; furthermore, we focus on the asymptotic behaviour of the R i s in the proximity of the spatial (y = 0) and temporal (s = 0) boundaries. Note that within the present perturbative expansion an algebraic behaviour ∼ x α with power α = α 0 + α 1 g + O g 2 is signaled by a logarithmic term, because
Recalling that d = 4 − ǫ and that we are employing dimensional regularisation [31], the integral over l in the previous equation has actually to be interpreted as with A =x ϑ/(1 − ϑ) and B =ỹ (1 − ϑ)/ϑ. Due to the fact that the integral is regular for ǫ → 0, we can conveniently fix d = 4. We now use the identity 2f ± (Al)f ± (Bl) = cosh ((A + B)l) ± cosh ((A − B)l) and express the hyperbolic cosines as cosh x = ∞ m=0 x 2m /(2m)!; by recalling that 2 ∞ 0 dl(l 2m−2 e −l 2 − δ m,0 ) = Γ(m − 1/2) we are able to calculate the integrals over l and find With an additional change of variables (1 − ϑ)/ϑ = √xỹ β the integral over ϑ above becomes It can be proved that the resulting series is equation (31) is pointwise convergent forxỹ > 0 and that the only singular term forxỹ = 0 is the one with m = 0, i.e., where K ν is the modified Bessel function of the second kind, whose asymptotic behaviour is K 0 (z → 0) ∼ − ln z. For the special transition this implies that, in the limitxỹ → 0, equation (32) is dominated by the logarithmic singularity which, after an allowed multiplication by cosh (xỹ/2) ≈ 1, can be cast in the form As in the case of R 1 in equation (25), this behaviour can be traced back to the emergence of an algebraic behaviour at the Wilson-Fisher fixed point; this algebraic behaviour reproduces the one predicted for the surface scaling [7,11]. Accordingly, this term correctly and completely captures the surface divergence in the special case. As mentioned in section V (see, e.g., equation (53)), b can be actually expressed [7] in terms of the bulk and surface exponents β and β 1 , respectively, introduced in equation (4) which correctly reproduces the previously-known results for the ordinary transition [12].

New singularities and edge corrections to scaling
Finally, we consider the third term R 3 ; its expression is the same as in equation (28) with Γ (d/2 − 1) replaced by the incomplete gamma function γ(d/2 − 1, z 2 /(2τ )); by introducing the same change of variables as in equation (28) one arrives at Since γ(α, w → 0) vanishes as ∼ w α , it constitutes a sufficient regularisation to make the integral in l convergent. One can therefore set d = 4 from the outset, noticing that, correspondingly, γ(1, w) = 1 − e −w . Thus, this integral becomes with the same A and B as in equation (30), C = 2ϑ(1 − ϑ)/(ϑ +s) ands = s/∆t. In order to extract the possible logarithmic contributions localised at the boundaries x = 0, y = 0, s = 0, we rewrite I as and note that I 1 represents a more regular version of I (see equation (30)); thereby, for the special (+) case this term is always regular and one can entirely focus on I 2 . On the other hand, in the ordinary case f − (0) = 0 and I 2 vanishes. Again, by applying to I 1 what we have found above for I, one can restrict to considering the first order in the expansion of the hyperbolic functions, i.e., f − (Al)f − (Bl) = sinh (Al) sinh (Bl) ≃ ABl 2 . Remarkably, the analysis reported below for the special case can be similarly repeated for the ordinary transition. Making use of the identity and introducing µ = (1 − ϑ)/ϑ one finds The argument of the square brackets in the expression above behaves as −µ 2 for µ → 0, whereas for µ → ∞ it vanishes as µ −2 /s for s > 0 and approaches 1 − √ 3 for s = 0. Therefore, even in the absence of the exponential (i.e., for x = y = 0) the integral is convergent for every s > 0. We also notice that the integral is still finite for s = x = 0, y > 0, as the exponential regularises the behaviour at µ → ∞. Thus, a singularity can be obtained only for y = s = 0, independently of x. Because of this, we introduce the representatioñ in terms of u = s 2 + (ỹ 2 ) 2 , which acts as an effective "radial" coordinate in the space of directions orthogonal to the spatial and temporal boundaries, and the corresponding "polar angle" α. Note that we have also accounted for the different scaling of time and space, described by the dynamical exponent z = 2 + O(ǫ 2 ) [19]. Accordingly, one has Q(x, √ u cos α, u sin α) ≡ Q(u, α), where we do not indicate explicitly the dependence onx since it plays no significant role. As in the previous cases, we expect a logarithmic behaviour Q(u, α) = f (α) ln u+O (1) to emerge for u → 0. In order to highlight it and calculate the coefficient f (α), we derive this function with respect to u and introduce µ = γ/ √ u, which yields In these terms, f (α) is given by J 1 (0, α) + J 2 (0, α), provided it is finite. In order to calculate it, the former addend J 1 (0, α) can be rewritten as where J (α) denotes the contribution of the second term within the brackets (i.e., the square root). We then rewrite J 2 (0, α) as where the last equality follows from an integration by parts. This confirms that the divergence of Q(u, α) for u → 0 is indeed logarithmic in nature. Moreover, it proves that the coefficient f (α) = ( √ 3 − 1)/2 is actually independent of the choice of α, which means that the divergence is the same when approaching the edge from any "direction" in the y z − s plane. Thus, the divergent part (43) can be rewritten as which entails the emergence of a novel, algebraic law R(x, t, u) ∼ u −θ E in the vicinity of the edge (i.e., for u → 0) with θ E = θ E,0 +θ E,1 g+O g 2 . The actual exponent can be inferred by multiplying the expression above by (n + 2)g/6 which yields θ In our notation, they read detail proves indeed that this is the correct way of approaching the problem; as expected on these grounds, we have in fact identified a new logarithmic singularity in our perturbative expansion which is strictly localised at the edge and therefore can only modify the scaling behaviour of observables in its proximity (see, e.g., equation (48) with the definitions (44)). Further confirmation of the presence of edge effects has been sought in reference [5], which reports a numerical Monte Carlo study of the classical Ising model in three spatial dimensions (i.e., n = 1 and d = 3) with Glauber dynamics. The comparison between the theoretical predictions presented here and these numerical data requires some care: in fact, the predictions for θ E reported here are limited to the first-order corrections in a dimensional expansion, which are generally known not to be quantitatively accurate. However, their qualitative features are typically more robust: for example, their sign typically dictates the one of the entire expansion they belong to. In the present case, it is therefore convenient to design the simulations in such a way to probe the correlation and response functions in a time regime in which s is varied within the window t ≫ s ≫ y z because this analysis would highlight an algebraic decay C(. . . , s) ∼ s 1−θ (see equation (56)) in the absence of edge effects, or C(. . . , s) ∼ s 1−θ−θ E (see equation (55)) in their presence. (In passing we mention that A turns out to be or order 1 in the numerical simulations of reference [5].) Noticing that for the ordinary and special cases the exponent θ E in equations (49) and (51) takes opposite signs, one can look for differences between the two transitions, which are expected to display faster/slower power laws with respect to the regime s ≪ y z ≪ t, which is dominated instead by the initial-slip physics and invariably leads to observing C(. . . , s) ∼ s 1−θ . The results of reference [5] are in agreement with such predictions and therefore constitute evidence of the existence of the edge modifications to the scaling laws discussed above. As a further confirmation, a crossover from the former to the latter regime is also observed upon moving away from the surface and eventually leads to a collapse of the data on a master curve which is indeed independent of the surface transition undergone by the model.