Magnetic domain wall creep and depinning: a scalar field model approach

Magnetic domain wall motion is at the heart of new magneto-electronic technologies and hence the need for a deeper understanding of domain wall dynamics in magnetic systems. In this context, numerical simulations using simple models can capture the main ingredients responsible for the complex observed domain wall behavior. We present a scalar-field model for the magnetization dynamics of quasi-two-dimensional systems with a perpendicular easy axis of magnetization which allows a direct comparison with typical experimental protocols, used in polar magneto-optical Kerr effect microscopy experiments. We show that the thermally activated creep and depinning regimes of domain wall motion can be reached, and the effect of different quenched disorder implementations can be assessed with the model. In particular, we show that the depinning field increases with the mean grain size of a Voronoi tessellation model for the disorder.


I. INTRODUCTION
The study of field-induced magnetic domain wall motion in thin ferromagnetic films has received great attention during last decades. Basic research allowed for the promise of new technological developments relying on the motion of domain walls [1][2][3], and received a large impulse in reward. In particular, magnetic thin films with perpendicular anisotropy are good candidates for highdensity magnetic memory devices. One of the advantages in these systems is the narrow domain wall width, of a few tens of nanometers, and the relatively easy control of the domain wall position with external magnetic fields or electric currents [4,5]. Therefore, the prospective development of new technologies based on domain wall motion prompts to deepen the understanding of domain walls dynamics.
How a domain wall in a magnetic material moves is dictated by the interplay between the external drive, thermal fluctuations, ferromagnetic exchange which results in a domain wall elasticity, and the disorder present in the sample. The external force acting over a domain wall can be generically considered to be the result of the application of an external magnetic field favoring the growth of one of the domains separated by the wall. When the magnetic field is small, domain wall motion is strongly hindered by the disorder. The velocity of the domain wall is ruled by activation: where ∆E is a disorder dependent energy barrier, k B T the temperature energy scale (with k B the Boltzmann constant), and V d is a reference velocity corresponding to the vanishing of ∆E. The disorder energy scale depends on the external field as with k B T d a characteristic disorder energy scale, H d the depinning field where the energy barrier goes to zero, and µ the creep exponent (µ = 1/4 for magnetic thin films) [6][7][8]. Equations (1) and (2) imply the so-called creep law, ln V ∼ H −1/4 , which is valid for fields below the depinning field, H < H d . For fields just above the depinning field, H H d , universal power-law behavior for the velocity-field response is due to the underlying zero temperature depinning transition and can be observed in the finite temperature domain wall dynamics [9,10]. Above the depinning field, H > H d , the flow regime is encountered, where the velocity grows linearly with the field with m the mobility. The overall non-linear velocity-field response has been observed in a wide variety of magnetic materials with its universal features characterizing creep and depinning regimes well accounted for by three parameters: the depinning field H d , the depinning temperature T d and the velocity scale V d = V (H d ) [8,10,11]. The use of numerical models assists to account for the full domain wall dynamics. Simple models as the elastic line in disordered media has been useful to unveil universal features of domain wall motion [7,[12][13][14]. The approach of the elastic line has the great advantage of allowing to obtain very precise exponents describing the systems dynamics in the elastic limit, which connects with analytical results. However, the purely elastic description leaves behind several experimentally well known features of domain wall dynamics: topological defects, fingering, overhangs, bubbles, plasticity, multi-valuated interfaces. Even more, nucleation phenomena cannot be assessed with this approach, thus rendering impossible to recreate the vast majority of experimental protocols.
Besides, two-dimensional spin models, as Ising, XY, and Heisenberg, have been adapted for the study of creep and depinning in domain wall motion [15][16][17][18][19][20]. Such spatial models permit indeed to simulate bubble domains and domains with overhangs, but their intrinsic periodic pinning made these models not truly realistic or comparable to the experiments. E.g., most simulations of driven domain walls with these approaches were done for random-field instead of random-bond disorder type.
Moreover, micromagnetic simulations stand as a relevant technique to address material specific properties. They have been intensively used to capture domain walls static and dynamic features, particularly in low dimensions and small systems [21][22][23][24]. However, this approach being detailed and exhaustive, it is not always helpful to distinguish and individualize the essential ingredients ruling the domain walls dynamics. On the computational side, the main disadvantage of this technique is the large amount of resources or time needed for its simulation [25]. Micromagnetic simulations are mainly used to study glassy domain wall dynamics close to the depinning transition, and in most of the studies only the T = 0 K case is considered. However, recently this technique has also been used to address the creep regime of domain wall motion in Pt/Co/Pt thin magnetic films [26], where one needs to simulate extended domain walls, i.e. domain walls whose extent is far much larger than its internal width. Although the creep regime has been reached [26], some features that are not fully compatible with experimental observations have also been observed, as for example, two distinct creep regimes.
When possible, it is desirable for numerical models and methods to mimic experimental protocols. Polar magneto-optical Kerr effect microscopy (PMOKE) is commonly used to measure domain wall velocity [27][28][29][30][31][32][33]. In a typical experimental protocol, one or several nuclei are first created, which usually present a bubble-like configuration. Then, finite time magnetic field pulses are applied, impelling the original domains to grow. The measured domain wall displacement is proportional to the pulse duration, thus giving a measure of the domain wall velocity. The insight that these experimental techniques can provide are naturally limited by several experimental factors: the camera resolution, magnetic field pulse characteristics as maximum amplitude and minimum width, control of the sample temperature and sample characteristics as the defect density and disorder of the sample under study. Therefore, having a model capable of reproducing the experimental conditions is highly desirable and should allow one to reach more quantitative comparisons between experiments and simulations.
Here, we adapt a very well known model in statistical physics, a two-dimensional scalar field model with a double well potential, to describe the phenomenology of domain wall motion in thin ferromagnetic films. The model lays in a mesoscopic scale, between the elastic line and micromagnetic models, allowing to cover large spatial and temporal scales while preserving a fairly detailed control of system parameters. After presenting the model and key considerations to obtain domain wall velocities, we show that simulated velocity field characteristics display the well acknowledge shape in both depinning and creep regimes, including the µ = 1/4 creep exponent value. Furthermore, we investigate the dependence of the domain wall dynamics under different quenched disorders, stressing how the present model can be used to study geometrical properties of magnetic domains.

II. MODEL
We are interested in the study of magnetic domain wall dynamics in thin films with strong perpendicular anisotropy. In this kind of systems, the magnetic moment of the material is given by the time-dependent vector field m( ρ, τ ), where ρ and τ are the two-dimensional space and time coordinates, respectively. m( ρ, τ ) is constrained to point perpendicularly to the sample plane, that we are going to take as the x − y plane. When domains are nucleated in the sample, the magnetization inside domains will still point perpendicularly to the sample plane (zdirection), with the same magnitude as in the rest of the sample, but with a different orientation. In the domain wall region, typically much smaller than the domain region, the magnetization will change smoothly from one value of magnetization to the other. In a system with a strong perpendicular magnetic anisotropy, the magnetization's x and y components will be approximately zero in the whole sample, except for the domain wall region. As the universal domain walls glassy dynamics is independent of the domain walls magnetic structure we will consider the evolution of the magnetization z-direction, neglecting the contribution of the remaining magnetization components.
The scalar field ϕ( ρ, τ ) = m z ( ρ, τ ) will represent the value of the magnetization z-component, taking real values in the interval [−1, 1], at position ρ in the x−y plane. This scalar field is a non-conserved variable: it may alter its value without a corresponding flux. The evolution of such non-conserved scalar field can then be modeled, in the limit of strong perpendicular anisotropy and strong damping [34], through where Γ is a damping parameter, H is the free energy of the system that may contain different terms describing the interactions and disorder present in the system, and ξ( ρ, τ ) represents an uncorrelated thermal bath modeled as a white noise, with ξ( ρ, τ ) = 0 and ξ( ρ, τ )ξ( ρ ′ , τ ′ ) = 2ΓT δ(τ − τ ′ )δ( ρ − ρ ′ ), with T acting as an effective temperature [35]. Equation (4) is the simplest stochastic dynamical model in which a single non-conserved scalar field ϕ( ρ, τ ) is in contact with a constant temperature heat bath. It has been already used in related problems such as the formation of magnetic patterns [36,37] or geometric pinning in magnetic films [38,39]. We model the system free energy Hamiltonian H by following the modified φ 4 model, as discussed by Jagla in Refs. [36,37]. In our implementation the model has three main contributions, H = H loc + H rig + H ext , as described in the following. The local term, H loc , mimics the out-of-plane easy axis magnetization and thus favors the values ϕ = ±1. It is given by with α proportional to the out-of-plane magnetic anisotropy constant. A rigidity term discourages spatial variations of ϕ, with an intensity β proportional to the exchange stiffness constant. Finally, the external magnetic field is incorporated through the term with a positive H favoring the ϕ = +1 state. We introduce two supplementary features to this simple model. First, we consider a prescription from the micromagnetic approach ensuring saturation of the local magnetization, which amounts to adding a saturation term (1 − ϕ 2 ) multiplying the external field H (see Ref. [37] for a discussion). Secondly, we introduce structural quenched disorder by perturbing the value of α in the H loc term. Instead of α we now use (α + εζ( ρ)), with ζ( ρ) a short-range correlated random variable with uniform distribution in [−1, 1] and ε the intensity of the disorder. This implementation of the disorder is compatible with the so-called random-bond disorder. The value of (α + εζ( ρ)) is then a spatially fluctuating quantity giving the height of the two well potential, which controls the strength of the system anisotropy energy, and is a measure of the local field required to revert an isolated magnetic moment.
Using in Eq. (4) the Hamiltonian H = H loc + H rig + H ext with quenched disorder in the local term plus a saturation prescription, the evolution of the field ϕ( ρ, τ ) is given by In a sense, the model in Eq. (8) is a simplification of the phenomenological Landau-Lifshitz-Gilbert equation, which provides a widely acceptable micromagnetic description of the evolution of the local magnetic moment direction of the material. With some variations, it has been proven successful in modeling the magnetization of quasi-two-dimensional systems [34,[36][37][38][39].
For simplicity, under a linear transformation Eq. (8) can be reduced to the form The last equality is imposed in order to ensure the proper correlation of the new effective temperature variable. From now on, all results will be expressed in reduced units, r, t, and h. In order to numerically solve Eq. (9) and obtain φ( r, t), we work with discretized time and space variables. We define a two-dimensional square grid with L × L cells. In each cell, φ has a uniform value updated at each step of the calculus. For the time integration of the equation, we use the first-order numerical Euler method, with a time step of 0.1 and given initial values. In order to implement the semi-implicit method to stabilize the numerical solution, we go through a Fourier transformation on the space variables, evaluating the exchange term at t + ∆t rather than at t. For more details on the numerical solution of Eq. (9) the reader may refer to [36].

III. RESULTS
In this Section, we first describe the adopted protocol and how the velocity of the domain wall is computed. Then we present results within the creep regime of domain wall motion and discuss temperature effects and fitted parameters. Finally, we present results depending on how the quenched disorder is implemented in the model.

A. Domain wall velocity
To measure domain wall velocities we used the following protocol inspired by experiments. As the initial condition for all simulations, the scalar field φ( r, t) is set to the value −1 in all system cells except those cells inside a circle of radius R 0 , centered at the middle of the system, where it takes the value +1. This initial condition is then relaxed by letting the system evolve at zero field (h = 0) for a time ∆t 0 until the circle area reaches a stationary FIG. 1. Evolution of the effective domain radius (circles), when a field square pulse of h = 0.07 is applied (also shown, with dashed lines) in a system at zero temperature and with a uniform disorder. The straight black line is a linear fit of the data during the application of the field pulse, which slope is indistinguishable from the obtained domain wall velocity at this field, as ∆R/∆t (see text). In the inset image, the spatial distribution of φ for a system with L = 4096 cells is shown. Black color indicates the value φ = −1, while gray and white correspond to φ = +1. The gray circle corresponds to the initial domain (before the field pulse) and the white part is the growth of the initial domain after the field pulse.
value. In order to apply an external field promoting domain wall motion, a constant field pulse of intensity h is then applied during a finite time ∆t. Finally, during a time ∆t ′ 0 the system relaxes, evolving at zero field again. In a system of size L = 4096 with ε = 1 in Eq. (9), ∆t 0 = ∆t ′ 0 = 10 3 is enough to ensure that the domain area reaches a stationary value at zero field. These parameters are kept fix at that value throughout the rest of the numerical simulations. Note that this sequence of steps is equivalent to the sequence in which magnetic fields are applied to a sample in a PMOKE microscopy experiment, where first the sample magnetization is saturated in the −z-direction, a nucleation field is applied in order to generate a domain with magnetization in the z-direction and a square pulse is applied in order to accomplish the domain growth [28].
Domain wall velocities are hence computed measuring the increase in domain area during the application of the magnetic field pulse. The area of the domain corresponding to φ = +1, A + , is calculated and registered during the whole simulation. Assuming a circular shape for the domains, effective radius is computed as R = A + /π. The domain velocities are then estimated as v = ∆R/∆t. ∆R = R(∆t 0 +∆t+∆t ′ 0 )−R(∆t 0 ) is the effective domain radius computed as the difference between the effective domain radius before applying the field pulse and after a time ∆t ′ 0 following the field pulse. As an example, the effective domain radius evolution for a square field pulse of intensity h = 0.07 is shown in Fig. 1.
In order to be consistent with PMOKE experiments, it is important that numerical results for the velocity do not depend on domain size nor pulse duration. Therefore, we check that the measured velocities are stationary and independent of the domain size. Figure 2 presents results for two values of the applied field for different initial domain sizes, R 0 , and different durations of time pulses, ∆t. We find that if R 0 or ∆t are too small, velocities may be underestimated or overestimated, respectively, especially for small values of h close to the depinning field (see below). The underestimation of the velocities for small domain radius may be due to the domain curvature since the effective field sensed by the domain wall is corrected with a term proportional to the inverse of the domain radius (h ef f = h − c/R). This effect may not be assessed experimentally with PMOKE microscopy since it occurs at much smaller scales than the camera resolution. For instance, a typical domain wall width is ∼10 nm. The curvature effect according to Fig. 2 is important for R 0 100 simulation cells, that are equivalent to 1 µm by following the transformations of Eq. (11), with the domain wall width estimated as β/α =10 nm. On the other side, the overestimation of the velocities at small durations of the field pulse may be due to a memory effect of the domain walls [14]. Henceforth, to ensure a representative value for the velocity, we use R 0 = 10 3 for all simulations and a carefully chosen value of ∆t for each field, in the range 10 3 to 5 × 10 6 .
When a system at zero temperature and no disorder is considered, a trivial linear behavior for domain wall velocities is found, as shown in Fig. 3 with open squares, which corresponds to a linear flow regime. The mobility m of the domain wall is the proportionality factor between velocity and field and depends on its internal structure. The particular form of the domain wall, i.e. the domain wall profile, needs to be considered in order to estimate the mobility. It is interesting to note that an estimation of domain wall velocities in the flow regime can be extracted from Eq. (9). Lets consider a system of size A with a φ = +1 single domain of area A + ; correspondingly the rest of the system, A − = A − A + , has φ = −1. The total system magnetization M can thus be written as where the integral is taken over the whole system. Taking time derivatives in Eq. (11), and using that A = A + +A − , one obtains To further simplify the problem, we can consider a rectangular portion of the system, of length l, containing one domain wall at a position x 0 (t), and hence A + = lx 0 (t).
Under the action of an applied field h, the domain wall velocity can be obtained as Equations (12) and (13) therefore relate the domain wall velocity with the time evolution of the scalar field φ(t), which is described by Eq. (9). For the case of a system without disorder (ε = 0) at zero temperature (T = 0), as the one considered in Fig. 3, the velocity can be expressed in a simple form: where the integral was solved by using a functional form of the domain wall profile given by the expression For this simple model, the mobility is thus equal to the domain wall width δ. Fig. 3(c) presents three domain wall profiles φ(x) for the direction (x, L/2), taken with a time difference of t = 10, corresponding to the case h = 1 and without disorder at zero temperature. These profiles can be well fitted with Eq. (15), giving a value δ = 1.4 [40]. In Fig. 3 we show with a dashed line the linear relationship of Eq. (14) between v and h, using δ = 1.4 for the mobility, showing a fairly good agreement with the measured velocities in the so-called flow regime. When disorder is considered (at zero temperature) the same linear behavior is observed at large field values, as shown in Fig. 3(a) (circles) for a uniform disorder with ε = 1. However, when the field is decreased the domain wall movement is strongly impeded due to the presence of disorder, resulting in a strong decrease of the velocity below h ≈ 0.06, as shown in Fig. 3(b). A closer inspection of this behavior is shown in Fig. 4. At zero temperature a power-law vanishing of the velocity is expected when the depinning field is approached from above, v ∼ (h − h d ) β , with h d the depinning field and β the depinning exponent (see Ref. [14] and references therein). In order to estimate the depinning field from the numerical results,  Fig. 4), the depinning field corresponds to the point where the theoretical β = 0.245 value [42] is reached, resulting in h d = 0.0598. This value is indicated with a pointed vertical line in the main panel of Fig. 4.

B. Creep and depinning regimes
Domain wall velocities for finite temperature values as a function of the applied field are shown in Fig. 4 for two different non-null temperatures, T = 0.001, 0.01. Stationary velocities values are observed at fields smaller than the depinning field, h < h d (T = 0), since temperature allows the activation over energy barriers, as expected in the creep regime. As indicated by Eq. (1) and the field dependence of the energy barrier, a linear relationship between ln v and h −1/4 should be observed in the creep regime. Such a creep plot is shown in Fig. 5 for the two finite temperature data sets. It shows that the numerical data is compatible with a creep exponent µ = 1/4 for the smaller field values. The inset of Fig. 5 shows the dependence of the velocity with the pulse duration ∆t in a creep plot, showing how the stationary velocity limit is reached at increasing ∆t for low fields. This should be carefully taken into account in numerical simulations.
In order to discern how far one can progress on the comparison between the model and experimental results, we use the same fitting procedure as recently used for experimental data [10,11]. This allows one to extract the three key parameters describing the glassy dynamics within creep and depinning regimes: the depinning field h d , the depinning temperature T d , and the velocity scale v d = v(h d ). The fitting procedure is described in detail in Ref. [11]. In brief, the depinning field and the velocity scale are first estimated using the inflection point of the v(h) curve, which allows one to estimate the depinning temperature from the slope of the creep plot. Then the full model, Eqs. (1) and (2), is fit allowing to adjust the three values. Finally a fine tunning is achieved using that, just above depinning, the velocity presents signals of the zero temperature depinning transition [43], with y 0 = 0.65 a fixed universal constant and ψ = 0.15 the thermal rounding exponent [9,10,13]. Results of the fit using the creep law, Eqs. (1) and (2), and the depinning transition scaling, Eq. (16), to the velocity-field numerical data are plotted in Fig. 6  This feature of the model is due to a large crossover between the creep and the flow regimes, that is also observed in velocity-force curves obtained with the elastic line model [44]. Overall, we have shown that the numerical data can be fit using the same fitting procedure as used to deal with experimental data.

C. Models of disorder
Finally, since the specific model of disorder is, at least partially, responsible of the domain wall dynamics, and in order to stress potential applications of the present model, we show how the velocity-curve depends on the underlying disorder model. We then study the variation of domain wall velocities using three different disorder models. In the first disorder type, already presented, the values assigned to the disorder (ζ( r) in Eq. (9)) were randomly chosen from a uniform distribution over the range [−1, 1], independently for each numerical cell in the system. For the second disorder type, we use a Voronoi tessellation of the system with N V = 1.5 × 10 6 Voronoi grains and give for each grain a constant ζ value between -1 and 1 from a uniform distribution. Finally, as a third disorder model, a filtered disorder, is built by using a standard low pass filter over an independent random uniform distribution. These three disorder types are shown in Fig. 7 for a fraction of 50 × 50 cells of the two-dimensional system.
Numerical results for the velocity-field curve for the different disorder models are shown in Fig. 8(a). As can be observed, velocity scales are visibly dependent on the type of implemented disorder. For a given field, velocity decreases when the considered disorder model passes from the filtered disorder to the uniform disorder and to the Voronoi disorder. In fact, the lowest depinning field is obtained for the filtered disorder, while the greatest depinning field corresponds to the Voronoi tessellation model. One can also observe that although h d changes with the disorder, and T d and v d probably too, the general shape of the velocity-field curve seems to be preserved. This means that universal features, as the critical exponents, are not presumably changing. In fact, the creep plot presented in Fig. 8(b) shows that the creep regime for the three different disorder models can be well described using the universal creep exponent µ = 1/4. The present model can also be used to investigate the effect of different disorder types on domain walls geometrical properties. Figures 8(c-d-e) show the shape of the domain for the three studied disorder types, all obtained at the same velocity within the creep regime, as indicated by full large symbols in Figure 8(b). A simple inspection shows that the roughness of the domain's shape increases with the value of the depinning field, depending on the type of disorder model used.
Until now, we showed qualitatively how different domain geometries and depinning fields may be obtained by changing the disorder implementation, but the comparison was not fair in the sense that the uniform dis- order and the filtered disorder have different correlation lengths and intensities. On the other side, a uniform disorder can be recast as an extreme case of a Voronoi tessellation, where the smallest possible area σ for the Voronoi grain sizes is considered. Thus, the amplitude of the disorder is not changed but the correlation length is. To explore deeper on this point, we tested two other Voronoi mean grain sizes, by generating Voronoi tessellations of N V = 5 × 10 6 and N V = 10 7 cells. The four Voronoi tessellations correspond to a mean area of the grains of σ = 1, σ ∼ 1.6, 3.4, 11.2, when decreasing the number N V , respectively. Velocity-field curves are shown in Fig. 9 for the four elections of N V . The main feature to highlight is that smaller depinning fields are obtained for smaller grain sizes of the Voronoi tessellation. This dependence of the depinning field with the Voronoi mean grain sizes was observed before in micromagnetic simulations [45], although we show here that the geometrical properties of domain walls also depend on the Voronoi mean grain size ( Fig. 9(b-e)).

IV. CONCLUSIONS
In summary, we have presented a study of domain wall dynamics in thin magnetic films using a versatile effective two-dimensional model. The model can be recognized as a generalization of the φ 4 -model of statistical mechanics, commonly used to study phase transitions and critical phenomena [35]. It includes exchange interactions, external field, effective temperature and disorder and can be easily extended to consider dipolar interactions.
With the aim of conciliating numerical simulations with experiments, we treated the numerical system using the same protocol as in experiments. For example, the same sequence of applied magnetic field pulses was considered, and we discussed how to obtain stationary velocity values, independent of the initial domain size and the pulse duration. We showed that the lowest field velocity results are compatible with the thermally activated creep regime. This is an important numerical milestone, it opens the possibility to study creep dynamics at large length and time scales with a simple but realistic and material-parameters-tunable numerical model. We showed that our numerical results are well described by critical exponents commonly used in thin magnetic systems: µ = 1/4, β = 0.245, and ψ = 0.15. This suggests that, in the range of parameters explored here, the full spatiotemporal description of the domain wall is compatible with the quenched Edwards-Wilkinson universality class. A more quantitative comparison, especially viable for the depinning regime, is left as a second step, however. In particular, a question mark is opened to know to which extent elastic depinning scaling will hold in situations where plasticity, bubbles, and overhangs become more dominant. Furthermore, with the same fitting procedure used to analyze experimental results, we obtained values for the key non-universal parameters needed to describe domain wall dynamics in the creep and depinning regimes. Different disorders were finally considered, stressing the versatility of the model. Properly modeling the disorder landscape in thin ferromagnets is key to the understanding of domain walls dynamics and its influence on materials design. The Voronoi tessellation disorder model appears as a tractable model in this direction [45]. In particular, we found that within the Voronoi tessellation disorder model the depinning field and the domain wall roughness both increase with the mean size of the Voronoi grains. We expect that fitting experimental results with the presented model would provide experimental values for the parameters characterizing the disorder, such as the mean grain size and the energy scale of the disorder landscape. Furthermore, an overall systematic exploration of disorder-type effects on phase field and micromagnetic models for domain wall dynamics is somehow missing in the field, and the approach here presented appears as a good starting point on this direction. Prominent features of the studied model are its adaptability to realistic model parameters and versatility to study many different experimentally-inspired protocols that may be difficult to actually perform in the lab. For example, besides the domain wall dynamics, a careful study of domains nucleation for different disorder types with varying intensity can be performed with the same model. Financial support by the France-Argentina Project ECOS-Sud No. A12E03 is acknowledged. This work was partly supported by the Argentinian projects PIP11220120100250CO (CONICET), PICT2016-0069 (MinCyt) and UNCuyo Grants No. 06/C490 and 06/C017.