Nonequilibrium Critical Behavior of Magnetic Thin Films Grown in a Temperature Gradient

We investigate the irreversible growth of $(2+1)-$dimensional magnetic thin films under the influence of a transverse temperature gradient, which is maintained by thermal baths across a direction perpendicular to the direction of growth. Therefore, different longitudinal layers grow at different temperatures between $T_1$ and $T_2$, where $T_1<T_c^{hom}<T_2$ and $T_c^{hom}=0.69(1)$ is the critical temperature of films grown in homogeneous thermal baths. We find a far-from-equilibrium continuous order-disorder phase transition driven by the thermal bath gradient. We characterize this gradient-induced critical behavior by means of standard finite-size scaling procedures, which lead to the critical temperature $T_c=0.84(2)$ and a new universality class consistent with the set of critical exponents $\nu=3/2$, $\gamma=5/2$, and $\beta=1/4$. In order to gain further insight into the effects of the temperature gradient, we also develop a bond model that captures the magnetic film's growth dynamics. Our findings show that the interplay of geometry and thermal bath asymmetries leads to growth bond flux asymmetries and the onset of transverse ordering effects that explain qualitatively the shift observed in the critical temperature. The relevance of these mechanisms is further confirmed by a finite-size scaling analysis of the interface width, which shows that the growing sites of the system define a self-affine interface.


Introduction
The importance of thin film technology has been widely recognized in the realms of experimental and applied science, from the manufacture of electronics (layers of insulators, semiconductors, and conductors from integrated circuits) to optics (reflective and anti-reflective coatings, self-cleaning glasses, etc) and packaging (e.g. aluminium-coated PET films). Indeed, the increasing role of thin films in basic and applied research relies on the development and refinement of nanoscale deposition techniques such as sputtering and molecular beam epitaxy, which allow a single layer of atoms to be deposited at a time [1,2,3,4].
Since the growth temperature is one of the critical parameters in the formation of ordered thin films, several experiments have focused on the influence of a temperature gradient during film growth. In an early experiment by Tanaka et al. [5], magnetic Tb-Fe films were grown between two substrates with a temperature gradient, reporting the observation of perpendicular magnetic anisotropies and other gradient-driven structural features. More recently, Schwickert et al. [6] introduced the "temperature wedge method" where a calibrated temperature gradient of several hundred Kelvin was established across the substrate during co-deposition of Fe and Pt on MgO(001) and MgO(110) substrates. These experiments generated the L1 0 ordered phase of FePt, which is currently the leading candidate material for ultrahigh density heat assisted magnetic recording (HAMR) and bit patterned magnetic recording (BPMR) media ( [7,8] and references therein). Other experiments by Yongxiong et al. [9] have investigated the evolution of Fe oxide nanostructures on GaAs(100) by using a multi-technique experimental setup that included transmission and reflection high energy electron diffraction and scanning electron microscopy. In these studies, nanoscale epitaxial Fe films were grown, oxidized, and annealed using a gradient temperature method, which led to nanostripes with uniaxial magnetic anisotropy. As a result of the experimental advances on this field, many technological applications have been envisioned as well. For instance, magneto-optical recording studies of signal reproduction [10] have suggested that recording media having multiple magnetic layers in a transverse temperature gradient may suppress magnetic noise from tracks adjacent to the target track during information storage and reproduction [11].
From a theoretical perspective, gradients have been studied extensively in the context of diffusion processes and later extended to thermal conductivity and heat conduction problems. The so-called gradient percolation method was originally introduced to study percolation transition models where the density is the control parameter [12] and later applied to a variety of problems, such as fractal diffusion fronts [13,14,15], overlapping disks in a concentration gradient [16], bond percolation for the Kagomé lattice [17], invasion percolation under gravity [18], porous media [19], as well as in the study of vegetation distribution [20]. Very recently, the gradient method has been extended as a powerful tool to study first-and second-order irreversible phase transitions in far-from-equilibrium systems such as the Ziff-Gulari-Barshad model and forest-fire cellular automata [21,22]. In magnetic systems, damage spreading processes in a temperature gradient [23] and studies of several one-dimensional models [24,25,26,27] have been followed by the investigation of the kinetic Ising model in two dimensions under a variety of dynamics [28,29,30,31,32].
Within the broad context of these recent experimental and theoretical investigations, the aim of this paper is to study the irreversible growth of magnetic thin films in a temperature gradient and to provide a full characterization of the gradientinduced critical phase transition. The magnetic thin film growth process under far-from-equilibrium conditions is investigated by using the so-called magnetic Eden model (MEM) [33,34,35], an extension of the classical Eden model [36] in which particles have a two-state spin as an additional degree of freedom. Earlier studies have shown that, growing in (d + 1)-dimensional strip geometries in homogeneous thermal baths, MEM films are noncritical for d = 1 [34]. In contrast, for d = 2 they undergo an order-disorder phase transition that takes place at T hom c = 0.69(1) in the thermodynamic limit. The critical exponents are ν hom = 1.04 (16), γ hom = 2.10(36), and β hom = 0.16 (5), which intriguingly agree within error bars with the exact exponents for the kinetic Ising model [34]. Since the MEM growth process is irreversible and newly deposited particles are not allowed to flip and thermalize once they are added to the growing cluster, the observed correspondence between the MEM and the equilibrium Ising model remains puzzling.
In this work, we focus on the critical case (i.e. d = 2) and show that, when applying a transverse temperature gradient maintained by thermal baths between temperatures T 1 and T 2 , where T 1 < T hom c < T 2 , the system undergoes a continuous phase transition at a higher critical temperature (T c > T hom c ) and with different critical exponents. We also develop a growth bond model and show the existence of bond flux asymmetries caused by the interplay of geometry and thermal bath asymmetries, which shed some light on the growth dynamics and explain qualitatively the shift observed in the critical temperature. The growth bond model analysis is further supported by the fact that the growing interface is self-affine, thus ensuring that the growing sites are correlated at all size scales.
The rest of the paper is laid out as follows. In Section 2, we introduce the model and describe the Monte Carlo algorithm used to simulate MEM thin films in a temperature gradient. In Section 3, we present our results and a discussion. Finally, Section 4 consists of concluding remarks.

The Model and the Monte Carlo Simulation Method
The MEM in (2 + 1)−dimensions is studied in the square lattice by using a rectangular geometry L x × L y × L z , where L z ≫ L x = L y ≡ L is the growth direction. The location of each spin on the lattice is specified through its coordinates (x, y, z) . The starting seed for the growing cluster is a plane of L×L parallel-oriented spins placed at z = 1 and cluster growth takes place along the positive longitudinal direction (i.e., z ≥ 2). Across one of the transverse directions (the y−axis), a temperature gradient is applied by thermal baths at fixed temperatures linearly varying between T 1 and T 2 . Therefore, in our setup each layer at fixed y is subjected to a constant temperature T (y) = T 1 +(T 2 −T 1 )(y−1)/(L−1) maintained by a thermal bath. We adopt open boundary conditions along the y−direction, while continuous boundary conditions are considered along the x−direction.
Clusters are grown by selectively adding two-state spins (S xyz = ±1) to perimeter sites, which are defined as the nearest-neighbor (NN) empty sites of the already occupied ones. Let us recall that the substrate is a 3D cubic lattice and therefore each lattice site in the bulk has 6 NN sites. Considering a ferromagnetic interaction of strength J > 0 between NN spins, the energy E of a given configuration of spins is given by where the summation xyz, x ′ y ′ z ′ is taken over occupied NN sites. The Boltzmann constant is set equal to unity throughout, and both temperature and energy are measured in units of J. The probability for a perimeter site at (x, y, z) to be occupied by a spin is proportional to the Boltzmann factor exp(−∆E/T ), where ∆E is the change of energy involved in the addition of the spin and T is the temperature at the perimeter site. At each step, all perimeter sites have to be considered and the probabilities of adding a new (either up or down) spin to each site must be evaluated. Using the Monte Carlo simulation method, after all probabilities are computed and normalized, the growth site and the orientation of the new spin are both simultaneously determined by means of a pseudo-random number. Notice that the MEM's growth rules require updating the probabilities at each time step and lead to very slow algorithms compared with analogous equilibrium spin models. Since the observables of interest (e.g. the mean transverse magnetization along the x−direction and its higher moments) require the growth of samples with a large number of transverse planes of size L × L, clusters having up to 10 9 spins have typically been grown for lattice sizes in the range 12 ≤ L ≤ 96. Also, let us point out again that, although Eq. (1) resembles the Ising Hamiltonian, the MEM is a nonequilibrium model in which new spins are continuously added, while older spins remain frozen and are not allowed to flip, detach, nor diffuse.
As in the case of the classical Eden model, the magnetic Eden model leads to a compact bulk and a self-affine growth interface [33] (see Sect. 3.4 for a detailed finitesize scaling analysis of the interface width). The growth front may temporarily create voids within the bulk, usually not far from the rough growth interface. However, since the boundaries of these voids are also perimeter sites, they ultimately become filled at some point during the growth process. Hence, far behind the active growth interface, the system is compact and frozen, and the different quantities of interest can thus be measured on defect-free transverse planes.
Notice that the growth of magnetic Eden aggregates in (2 + 1)-strip geometries is characterized by an initial transient length ℓ T ∼ L 2 (measured along the growth direction, i.e. the z−axis) followed by a nonequilibrium stationary state that is independent of the initial configuration [34]. We considered starting seeds formed by L×L up spins (i.e. S xy1 = 1) but any choice for the seed leads to the same stationary states for z ≫ ℓ T . By disregarding the transient region, all results reported in this paper are obtained under stationary conditions. Also notice that, since the films are effectively semi-infinite and the substrate length along the growth direction plays no role, the only characteristic length is the transverse linear size L.

Gradient-Driven Continuous Pseudo-Phase Transitions in Finite-Size Films
Let us begin by considering a fixed gradient between temperatures T 1 = 0.5 and T 2 = 1.5. The effect of changing this gradient will be discussed later. Figure 1 shows the snapshot of a longitudinal slice for a fixed value of the transverse coordinate x. The temperature gradient is applied along the transverse y−direction, while the  system grows from left to right along the longitudinal z > 0 direction. Notice that the bulk grows compact, because although voids and holes in the bulk may eventually occur, they ultimately become filled at some point during the growth process. The bottom layers grow in contact with thermal baths at cold temperatures, which favor the formation of well-ordered spin domains. In contrast, the top layers grow in contact with hot thermal baths that promote bulk disorder. As will be shown below, the interplay of model growth dynamics, geometry, and thermal bath asymmetries lead to the onset of gradient-driven order-disorder critical phase transitions that can be quantitatively characterized.
In order to take into account the asymmetries introduced by the temperature gradient, we can quantify the degree of order in the system by considering the magnetization of transverse columns at constant temperature (i.e. along the x-axis): (2) Figure 2 shows the probability distributions of m for a system of linear size L = 64 growing in a temperature gradient between T 1 = 0.5 and T 2 = 1.5. Notice that different plots correspond to different layers and, therefore, to different temperatures within the gradient's range, as indicated. As expected for a continuous order-disorder phase transition, the low-temperature distributions are bimodal and peaked at the spontaneous magnetization m = ±m sp (0 < m sp < 1). As the temperature is increased, the peaks approach each other and merge smoothly, ultimately leading to a Gaussian distribution peaked at m = 0 for high temperatures, which is characteristic of the disordered phase. Indeed, the smooth shift of the distribution maxima across T ≃ T c (L), from the low-temperature nonzero spontaneous magnetization m = ±m sp to the high-temperature Gaussian centered at m = 0, is the signature of true thermally-driven continuous phase transitions [37]. Notice that the distributions in Figure 2 are symmetrical, since there exists a finite probability for fluctuations to grow and switch the magnetization from m ≃ +m sp to m ≃ −m sp and viceversa. Since Monte Carlo simulations are restricted to finite samples, the standard procedure to avoid these shortcomings due to finite-size effects is to average the absolute value of the order parameter [38]. In this context, the appropriate order parameter is the mean absolute magnetization of transverse columns at constant temperature, i.e.
where ... z denotes averages along the growth direction z > 0 within the stationary region. Figure 3 shows plots of the mean absolute magnetization as a function of : Susceptibility as a function of the layer temperature for different system sizes in the range 12 ≤ L ≤ 96. As expected for a critical system, the peaks become sharper and higher as L is increased.
the layer temperature for different system sizes in the range 12 ≤ L ≤ 96. For any given system size, at low temperatures the system grows ordered and the magnetization is close to unity, while at higher temperatures the disorder sets on and the magnetization becomes significantly smaller. However, fluctuations due to the finite system size prevent the magnetization from becoming strictly zero above the critical temperature, so the transition between the low-temperature ordered phase and the high-temperature disordered phase becomes smoothed out and rounded. As expected, larger systems are less affected by finite-size effects and display sharper transitions. Strictly speaking, however, these results just show evidence of pseudo-phase transitions, which might be precursors of true phase transitions taking place in the thermodynamic limit. In the following, we will characterize in more detail this pseudo-critical phenomenon by measuring other observables on finite-size systems. In the next Subsection, we will use standard finite-size scaling procedures to establish the existence of a non-trivial critical temperature in the L → ∞ limit, as well as to calculate critical exponents that describe the behavior of the infinite system at criticality.
Let us now consider the magnetic susceptibility χ, given by For equilibrium systems, the susceptibility is related to order parameter fluctuations by the fluctuation-dissipation theorem. Although the validity of a fluctuationdissipation relation in the case of a nonequilibrium system is less evident, we will assume Eq. (4) to hold also for the MEM. Indeed, as shown in earlier studies of nonequilibrium spin models [39,40], this definition of χ proves very useful for exploring critical phenomena under far-from-equilibrium conditions. In Section 3.2, we will characterize the critical behavior in the thermodynamic limit through critical exponents and finite-size scaling relations by applying the equilibrium theory to our far-from-equilibrium model. Figure 4 shows plots of χ vs T for different system sizes, as indicated. As with the thermal dependence of the order parameter shown in Figure 3, the peaks of the susceptibility become rounded and shifted, indicating the existence of pseudo-phase transitions in finite-size MEM thin films. Increasing the system size, the peaks become sharper and higher, as expected for a critical system.
Since the results presented thus far considered a fixed gradient between the temperatures T 1 = 0.5 and T 2 = 1.5, let us now investigate the effects of changing the gradient span ∆T ≡ T 2 − T 1 . With this aim, we keep the same temperature for the thermal bath at the cold end (T 1 = 0.5) and consider a substantially higher temperature for the thermal bath at the hot end (T 2 = 2.5). Figure 5 compares the mean absolute magnetization for these two different gradient ranges (i.e. ∆T = 1, 2) for systems of size L = 12 and L = 96. We observe that increasing the gradient span shifts the magnetization profiles towards higher temperatures. That is, the temperature of a given layer does not uniquely determine its degree of order, since the mean magnetization of the layer also depends on the overall gradient span under which the film grows. Similar shifts towards higher temperatures are also seen in higher-order moments of the order parameter probability distributions, such as the susceptibility and Binder's fourth-order cumulant (not shown here for the sake of space). Alternatively, we can compare systems of different sizes and gradient spans such that the local gradients δ ≡ ∆T /L are the same. The inset to Figure 5 shows a comparison between a system of size L = 16 and gradient span ∆T = 1 (dashed lines) and another system of size L = 32 and gradient span ∆T = 2 (dotted lines), both of which have the same local gradient δ = 1/16. The solid line corresponds to the mixed case L = 32 and ∆T = 1 (i.e. δ = 1/32). We observe that the systems with the same local gradient have similar magnetization in layers at intermediate temperatures (i.e. approximately in the range 0.8 ≤ T ≤ 1.3). However, the magnetization profiles for equal-δ systems differ noticeably in layers closer to the borders of the sample. The arrows in Figure 5 show that the shifts for smaller systems are significantly larger than the corresponding shifts for larger systems. Defining the finite-size pseudo-critical temperature T c (L) as the temperature corresponding to |m| = 0.5, the shift for L = 12 is ∆T c = 0.18, while the corresponding shift for L = 96 is ∆T c = 0.07. This observation suggests that differences arising from changing the gradient span might just reflect finite-size effects that vanish in the L → ∞ limit. In the next Subsection, we study the critical behavior of MEM films and confirm that, in fact, these differences are merely finite-size effects that become irrelevant in the thermodynamic limit.

Characterization of Gradient-Driven Critical Behavior in the Thermodynamic Limit
So far, we have found evidence for the existence of a gradient-driven order-disorder phase transition from the analysis of order parameter probability distributions, the order parameter mean absolute value (magnetization) and its fluctuations (suscep- tibility) in the growth of finite-size magnetic films. In this Subsection, we apply standard finite-size scaling techniques to show the existence of this phase transition in the thermodynamic limit and to determine critical exponents that characterize the system's critical behavior and universality class. The Binder cumulant, defined by is a fourth-order cumulant dependent on the variance and the kurtosis of the order parameter probability distribution. One important property of the Binder cumulant is that, for large system sizes, the low-temperature, ordered region tends to the value 2/3, while the high-temperature, disordered region tends to 0. Thus, in the thermodynamic limit, the function becomes discontinuous exactly at the critical temperature. Moreover, since for second-order phase transitions, the scaling prefactor of the cumulant is independent of the sample size, plots of U 4 versus the control parameter lead to a common (size-independent) intersection point that corresponds to the location of the critical value of the order parameter in the thermodynamic limit [41]. Figure 6 shows the Binder cumulant as a function of the layer temperature for different system sizes in the range 12 ≤ L ≤ 96. The Inset to Figure 6 shows a detailed view of the same data for the largest lattices (32 ≤ L ≤ 96), where the intersection region is indicated by a grey vertical strip. Based on this observation, we determine the critical temperature in the L → ∞ limit as T c = 0.84 (2). Interestingly, this value is significantly higher than the corresponding critical temperature for the MEM growing in an homogeneous thermal bath (i.e. in the absence of a temperature gradient), namely T hom c = 0.69(1) [34]. In the next Subsection, we will explore the growth dynamics and explain qualitatively this shift in the critical temperature as due to ordering effects caused by a net transverse growth bond flux induced by thermal asymmetries.
Notice also that, by fixing the temperature range (i.e. the temperatures T 1 and T 2 ) and increasing L, we are effectively considering different gradients δ = (T 2 −T 1 )/L that become smaller as L is increased. Figure 6, which shows a fixed point in the Binder cumulants as the gradients are changed, provides therefore quantitative evidence for the existence of a gradient-independent phase transition taking place at the temperature T c .
According to the finite-size scaling theory, developed for the treatment of finitesize effects at criticality under equilibrium conditions [42,43], the difference between the true critical temperature, T c , and the effective pseudo-critical one, T c (L), is given by where ν is the exponent that characterizes the divergence of the correlation length at criticality. As mentioned above, we define the finite-size pseudo-critical temperature T c (L) as the temperature corresponding to |m| = 0.5. Let us point out that, given the lack of a comprehensive theory of non-equilibrium phase transitions, concepts and definitions developed in the context of equilibrium phenomena are customarily borrowed and applied to far-from-equilibrium phenomena as well. For a review of standard methods, see e.g. [44,45,46]. Indeed, although this approach is ad-hoc and lacks the theoretical foundations of equilibrium systems, it has been used extensively in the literature and has become a powerful means of advancing our knowledge within the realm of non-equilibrium phenomena. For instance, numerical methods such as Monte Carlo simulations or series expansions are restricted to finite systems and it is therefore important to understand how far finite-size effects influence the properties of the system. As known from equilibrium statistical mechanics, finite-size effects are particularly strong close to the critical point, where the spatial correlation length becomes comparable with the linear dimensions of the system. By introducing the system size as an additional parameter, finite-size scaling laws are used to characterize the steady state of finite far-from-equilibrium systems through appropriate scaling exponents, such as, for instance, the exponent ν in Eq.(6) above. Moreover, this procedure allows us to define universality classes of non-equilibrium systems, as reviewed e.g. in Refs. [47,48,49]. Figure 7 shows log-log plots of the finite-size pseudo-critical temperatures T c (L) as a function of the inverse of the system linear size, L −1 , for different gradients and system sizes, as indicated. By rewriting the finite-size scaling relation as we performed different least-squares fits to the data using the mean, upper-bound and lower-bound values for the critical temperature in the thermodynamic limit. The nonlinear least-squares fitting procedure was implemented using the Levenberg-Marquardt minimization method [50] and the results from each independent fit are reported in Table 1. The errors in the table are determined by the fitting algorithm and take into account the statistical errors for each datapoint. Figure 7 shows that the finite-size scaling relation fits the data very well within error bars for the range of values for the critical temperature that was derived from the intersection of Binder's cumulants. From these fits, we obtain the critical exponent ν = 1.53 (6), where the error bars reported reflect the errors derived from the evaluation of T c as well as the statistical errors. Notice that the data for different gradients tends to converge in the L → ∞ limit, therefore confirming that differences arising from changing the   gradient are finite-size effects. On the other hand, finite-size scaling theory predicts that plots of |m| L β/ν vs |T − T c |L 1/ν for different lattice sizes should collapse near the critical region. The inset to Figure 7 shows the data collapse obtained by using β = 0.26 (that is determined from the hyperscaling relation, as explained below) with two separate branches corresponding to the low-and high-temperature regions. An additional characterization of the critical behavior of this system can be obtained by calculating the critical exponent γ, which describes the divergence of the susceptibility at the critical point. Using again the finite-size scaling theory [42,43], the exponent ratio γ/ν is related to the peak of the susceptibility measured in finite samples of size L by The symbols in Figure 8 correspond to the maxima of χ plotted against the system linear size for different gradients, as indicated, while the solid lines are fits to the data using the scaling relation from Eq. (8). Statistical errors for each datapoint are smaller than the symbol size in the Figure. The fitting procedure (which, as mentioned above, was implemented as a nonlinear least-squares algorithm using the Levenberg-Marquardt minimization method [50]) yields γ/ν = 1.66 (3), where the error bars reflect the statistical errors from the fits. Using this ratio and the value already obtained for ν, we determine γ = 2.54 (11). The insets to Figure 8 display plots of χL −γ/ν vs |T − T c |L 1/ν for ∆T = 1 and different lattice sizes in the range 32 ≤ L ≤ 96. Using the critical temperature as determined by the susceptibility peaks, the data collapse is shown separately for (a) the low-temperature branch and (b) the high-temperature branch. In the former case, data from low-temperature layers near T 1 = 0.5 depart from the collapse and have been removed. However, the collapse near the critical region is remarkable and agrees very well with the expectations from the finite-size scaling theory. By replacing the exponents ν and γ in the hyperscaling relation dν − 2β − γ = 0 with d = 2, we determine the exponent β = 0.26 (8), where the error is determined from standard error propagation applied to the hyperscaling relation. Recall that we anticipated this value of β when we considered the data collapse of the scaled magnetization (see the inset to Figure 7 above). The excellent data collapse near the critical region confirms the consistency and robustness of the obtained results.
As a summary, Binder's cumulant method and finite-size scaling analysis allowed us to characterize quantitatively the critical behavior of nonequilibrium magnetic films growing in a temperature gradient. We found that differences arising from changing the gradient are due to finite-size effects that vanish in the thermodynamic limit. The system's critical temperature is T c = 0.84(2), significantly higher than the critical temperature for films grown in an homogeneous thermal bath, T hom c = 0.69(1) [34]. The critical exponents are ν = 1.53 (6), γ = 2.54 (11), and β = 0.26 (8). Based on our findings, we conjecture that magnetic Eden films growing in a temperature gradient belong to a new universality class characterized by critical exponents ν = 3/2, γ = 5/2, and β = 1/4. In contrast, the critical exponents for magnetic Eden films grown in an homogeneous bath agree within error bars with the exact exponents for the Ising model in d = 2 [34], namely ν = 1, γ = 7/4, and β = 1/8.

Growth Bond Model and Bond Flux Asymmetries
In this Subsection, we explore the growth dynamics by means of a simple bond representation. Let us recall that the MEM's growth process adds new spins, which are deposited one by one to the growing cluster. Although voids and holes may form within the bulk, ultimately all sites become filled. Hence, to each pair of neighboring sites, we can assign a directed bond that points from the earlier occupied site to the later occupied one. The components of the bond flux field φ at a site (x, y, z) are defined as:  Figure 9 shows the x−, y−, and z−components of the mean bond flux φ as a function of the gradient span ∆T for L = 32 and T 1 = 0.5. As expected from the symmetry along the transverse x−direction, there is no net bond flux in x: φ x = 0 regardless of ∆T . For ∆T = 0, the system is also symmetric along y, so no net bond flux is observed. When a gradient is applied, however, this symmetry is broken. Since the growth probabilities depend on the Boltzmann factor exp(−∆E/T (y)), where T (y) is the layer's temperature, the thermal asymmetries introduced by the gradient favor spin deposition on the colder layers. This phenomenon is captured by the observed net bond flux φ y > 0. Indeed, as shown in Figure 9, the thermal asymmetries cause φ y to grow steeply up to φ y max ≈ 0.75 followed by a moderate decrease for larger gradients, which is due to the onset of bulk disorder within the hotter layers. Since the net transverse growth bond flux is directed from the ordered (cold) layers towards the disordered (hot) ones, this gradient-induced transverse ordering mechanism causes the system's critical temperature to increase from T hom c = 0.69(1) to T c = 0.84 (2). On the other hand, for ∆T = 0, φ z = 1 due to the longitudinal asymmetries in the substrate (i.e., the semi-infinite strip geometry constrains the system to grow along the z > 0 direction). However, when the transverse gradient is applied, two effects contribute to decrease φ z : (i) the onset of the transverse bond flux, which creates transverse domains in the active perimeter and causes some of the added spins to grow backwards; (ii) the bulk disorder induced in the hotter layers (which also causes φ y to decrease, as discussed above).
The mean fluxes shown in Figure 9 were averaged over sites at different temperatures. In order to gain further insight, let us now investigate the dependence of the bond fluxes on the layer temperature T (y) = T 1 + (y − 1) × (T 2 − T 1 )/(L − 1) (where 1 ≤ y ≤ L). Since we are considering different gradient spans for a fixed system size, it is actually more convenient in this case to plot the bond fluxes as a function of the transverse coordinate y. We already observed that the transverse bond flux along x is null regardless of temperature, so we will focus on the growth bond fluxes along the y− and z−directions. Figure 10(a) shows the upwards bond flux φ y (y) vs y for L = 32, T 1 = 0.5, and different gradient spans ∆T , as indicated. In the absence of a gradient, the flux is directed downwards at the bottom and upwards at the top, yielding zero net bond flux. Indeed, because of the open boundary conditions, empty perimeter sites at the confinement walls experience a missing-neighbor effect and the system grows preferentially along the center of the film as compared to the walls. When a gradient is applied, the flux grows steeply in the upwards direction, as expected. However, for larger gradients, the hotter thermal baths are capable of inducing disorder in the bulk and partially break the upwards bond flux on the upper layers. Indeed, this phenomenon causes the overall bond flux φ y (averaged over all layers) to decrease for large values of ∆T , as discussed above. Similarly, Figure 10(b) shows the forward bond flux φ z (y) vs y for L = 32, T 1 = 0.5, and different values of ∆T , as indicated. For ∆T = 0, there is a slight missing-neighbor effect for the sites near the confinement walls. Since forward growth is mostly driven by the substrate asymmetry, this effect is much less noticeable that in the flux along y. When the gradient is applied, two effects contribute to reduce the longitudinal flux, as discussed above. One of them, which is dominant for small gradients, is due to the formation of transverse domains along y, causing some backwards deposition when the bulk is filled in. The other mechanism, which is dominant for larger gradients, is due to the onset of bulk disorder in the hotter layers. Previously, we pointed out the fact that the gradient-driven order-disorder phase transition occurs at a temperature that is significantly higher than the corresponding critical temperature for the MEM growing in an homogeneous thermal bath (i.e. in the absence of a temperature gradient). The results presented in this Subsection allow a qualitative explanation for the shift in the critical temperature. Indeed, the temperature gradient breaks the transverse symmetry along the y−direction, causing the onset of a net transverse growth bond flux. This flux is directed from the ordered layers (that grow at T < T c ) towards the disordered layers (that grow at T > T c ), thus expanding the low-temperature ordering effects across layers at higher temperatures. This gradient-induced transverse ordering mechanism increases the system's critical temperature.

Scaling behavior of the growth interface
According to the analysis presented in Section 3.3, the system's critical temperature is increased due to gradient-induced transverse ordering mechanisms that originate in the cold layers. It could be argued, however, that, since the distance between the cold layer at T 1 and a layer at a fixed temperature T > T 1 becomes larger as L is increased (while keeping fixed the gradient span ∆T ), then the transverse flux effects may become weaker and eventually negligible in the large-L limit. In this Subsection we show that the growth interface is self-affine and its shape is stable and independent of size. Therefore, we confirm that the effects described in Section 3.3 are still relevant in the thermodynamic limit.
In order to track the evolution of the growth interface, we compute the position of the perimeter sites in the active region every time a monolayer of L 2 spins is deposited. The interface width at time t is defined by where the sum is taken over the N p perimeter sites in the active growth region, z i is the longitudinal coordinate of the i−th perimeter site, and z c = (1/N p ) i z i is the center of the interface. Figure 11 shows the interface width as a function of time for different lattice sizes in the range 12 ≤ L ≤ 96, where the time unit corresponds to the deposition of L 2 spins. After a short transient period, the interface reaches a stable saturation width, w sat , analogously to other surface growth phenomena [51,52]. The Inset to Figure 11 shows the dependence of w sat on the system size L, where the statistical errors are smaller than the size of the symbols. The fit to the data (solid line) shows that w sat ∝ L α , where the roughness exponent is α = 1.01 ± 0.01. That is, the saturation width scales linearly with the system size. By subtracting the interface center, z c , we can compute the average interface profile in the stationary regime, as shown in Figure 12. We find that, when scaled by the mean interface width, interface profiles for different system sizes collapse into a universal shape for the growth interface. This nearly-linear universal shape is qualitatively consistent with the roughness exponent α = 1, as it has quantitatively been determined. In passing, notice also that this universal profile shows a good agreement with the instantaneous snapshot displayed in Figure 1. Thus, we conclude that the active growth interface is self-affine and has universal features: the detailed analysis of bond flux asymmetries presented in Section 3.3 for a fixed lattice size (L = 32) remains valid for larger systems. In particular, our analysis focused on the influence exerted by the low-temperature layers into higher-temperature layers, which therefore is a gradient-induced growth mechanism still relevant in the thermodynamic limit. In addition to the roughness exponent α, self-affine interfaces are also characterized by a fractal dimension d f . In fact, within short length-scales such that ∆z ≪ ℓ, where ∆z is the longitudinal interface distance along the growth direction between two points separated by a transverse distance ℓ, the fractal dimension is d f = 2 − α [52]. In the long length-scale limit, moreover, the fractal dimension of the self-affine interface is d f = 1, irrespective of its roughness [52]. Therefore, we conclude that the growth interface of the MEM in a thermal gradient is d f = 1 at both short and long length-scales. Indeed, these conclusions are in full agreement with the nearly-linear shape of the growing interface that is consistent with a unitary fractal dimension at all length-scales, as seen in Figure 12.
Here, it is useful to compare the obtained results with those from related, wellstudied growth models. For the standard MEM growing in a homogeneous temperature bath at high temperatures, the attachment of spins is a stochastic (random) process that becomes independent of the interaction energy and the temperature [33]. Thus, in this limit, the growth interface of the standard MEM [33,34,35] agrees with that of the classic Eden model [36], which is well known to belong to the Kardar-Parisi-Zhang (KPZ) universality class [52]. The most accurate simulation results for the KPZ model in (2 + 1)−dimensions yield α = 0.393 ± 0.003 [53], which agrees well with some of the formerly reported values for KPZ [54] and the Eden model [55]. Using this value for the roughness, the fractal dimension of the growth interface of the MEM in an homogeneous thermal bath at high temperatures crosses over from d f ≃ 1.6 (at short length-scales) to d f = 1 (at long length-scales). Moreover, we can safely expect that the self-affine properties of the growth interface of the standard MEM (in an homogeneous bath) be independent of the temperature within a wide range around and below the critical temperature, at least insofar the occurrence of a layering/roughening transition, in the sense of that observed in the 3−dimensional Ising model [56,57], can be neglected. Therefore, we conclude that the fractal and self-affine characteristics of the growth interface of the MEM in a constant temperature bath are quite different than those of the same model growing in a temperature gradient.
On the other hand, due to the fact that the MEM grown under a temperaturegradient constraint exhibits a second-order transition, one may also consider the self-affine properties of the interface between the ordered and the disordered phases. Although for systems under equilibrium conditions a useful (alternative) approach is the evaluation of the damaged interface [23], the damage spreading technique can not straightforwardly be applied for the evaluation of an interface in an irreversible growth model. Furthermore, the implementation and application of a cluster counting algorithm [12,15,16,18,19,21,22,58,59] to our model would require formidable computational task that is beyond the aim of this paper. Additional shortcomings for this kind of calculations are the definition of the suitable cluster (e.g. Swendsen-Wang vs physical clusters) and the occurrence of noticeable corrections to scaling that one needs to evaluate in order to obtain reliable exponents [15], which is also a task that lies beyond our computational capabilities. However, from heuristic arguments based on the standard scaling relationship α ord−dis = ν/(1 + ν) [12,15,23], where α ord−dis is the roughness exponent of the order-disorder interface, we can conjecture that α ord−dis = 3/5, which yields a self-affine order-disorder interface with short length-scale fractal dimension d ord−dis f = 7/5. For comparison, by applying the same scaling relationship to the standard MEM growing in an homogeneous thermal bath, we obtain α ord−dis = 0.51 ±0.09 and d ord−dis f = 1.49 ±0.09. Thus, although the growth interface is very significantly affected by the temperature gradient compared to the thermally homogeneous system, the geometry of the order-disorder interface is not so markedly affected by the gradient constraint and the results for both systems agree within error bars.

Conclusions
In this work, we studied magnetic thin films growing under far-from-equilibrium conditions in (2 + 1)-dimensional strip geometries, where a temperature gradient is applied across one of the transverse directions. We modeled the thin film growth process by means of extensive Monte Carlo simulations performed on the magnetic Eden model (MEM), in which spins are deposited on a growing cluster with probabilities dependent on a ferromagnetic, Ising-like configuration energy.
Firstly, we studied the thermal dependence of order parameter probability distributions, the order parameter mean absolute value (magnetization), the order parameter fluctuations (susceptibility) and its higher moments (Binder cumulant) on finite-size magnetic films, which showed the existence of gradient-driven pseudophase transitions. Secondly, we applied Binder's cumulant method and finite-size scaling analysis in order to characterize quantitatively the critical behavior of MEM films growing in a temperature gradient. The system's critical temperature is T c = 0.84(2), significantly higher than the MEM's critical temperature when growing in an homogeneous thermal bath, namely T hom c = 0.69(1) [34]. The critical exponents are ν = 1.53 (6), γ = 2.54 (11), and β = 0.26 (8), which also differ from the MEM's exponents in the absence of a temperature gradient [34]. By changing the gradient span, we observed finite-size effects that vanish in the thermodynamic limit. Hence, the critical temperature and exponents are universal for MEM films growing in a temperature gradient. We also investigated the system's growth dynamics by means of a bond model. We found that the interplay of geometry and thermal bath asymmetries leads to growth bond flux asymmetries and the onset of transverse ordering effects that explain qualitatively the shift observed in the critical temperature. Finally, we analyzed the self-affine growth interface and obtained a collapse of the scaled average growth profile in the stationary regime for different lattice sizes, which shows that growth bond flux asymmetries play a relevant role in the model's growth dynamics even in the thermodynamic limit.
In the context of a great experimental and theoretical interest in magnetic systems growing in temperature gradients, as well as a wide variety of technological applications that benefit from these efforts, we hope that this work will contribute to the progress of this research field and stimulate further work.