Dynamic depinning phase transition in magnetic thin film with anisotropy

The dynamic pinning effects induced by quenched disorder are significant in manipulating the domain-wall motion in nano-magnetic materials. Through numerical simulations of the nonstationary domain-wall dynamics with the Landau–Lifshitz–Gilbert equation, we confidently detect a dynamic depinning phase transition in a magnetic thin film with anisotropy, which is of second order. The transition field, static and dynamic exponents are accurately determined, based on the dynamic scaling behavior far from stationary.


Introduction
In recent years, the domain-wall dynamics of magnetic materials is widely investigated in both experiments and theories, due to the potential applications in data storage and logic devices [1,2]. How to control the domainwall motion is a major focus of attention. The Landau-Lifshitz-Gilbert (LLG) equation is fundamental in the theoretical study of micromagnetics, and very important in understanding dynamical properties of magnetic materials, such as the domain-wall motion [3,4]. Driven by an external field, for example, the domain-wall propagation undergoes the Walker breakdown, and the underlying mechanism has been deeply investigated based on the LLG equation [5][6][7][8][9].
On the other hand, the dynamic pinning effects induced by quenched defects are significant in manipulating the domain-wall movement in nano-magnetic materials. For the domain wall pinning induced by notches in a strip, for example, the depinning phase transition has been investigated [10,11], which looks like a first-order one or a crossover, presumably due to the relatively simple pinning mechanism. For the domain wall pinning induced by disorder such as random fields, the depinning phase transition is typically of second order, and the experiments clearly indicate the existence of such a dynamic transition in magnetic thin films [12][13][14]. With a constant driving field H, the depinning transition occurs at zero temperature, and the order parameter is the domain-wall velocity. It is well known that when the external field H is above the Walker threshold H W , the domain-wall motion enters an oscillation mode [5,9]. Therefore, one usually considers the case that the transition field H c is below the Walker threshold H W . If the temperature is non-zero, the depinning transition becomes smeared, and the creep motion of the domain wall is thermally activated even below the transition field H c [12].
There have been theoretical efforts with the quenched Edwards-Wilkinson (QEW) equation and the Monte Carlo simulation, to understand the depinning phase transition of the domain-wall motion [15][16][17][18][19][20][21]. The QEW equation is a simple phenomenological model, where a domain wall is considered to be an elastic string, and detailed microscopic structures and interactions of the materials are not concerned. The elastic string is mathematically described by a height function h(y, t), where y is an one-dimensional coordinate, and t is time. Usually h(y, t) is assumed to be a single-valued function. For real materials, however, the domain wall is generally more complicated, such as with overhangs and islands, and a single-valued height function is not sufficient to describe the domain wall. The Monte Carlo method provides also only an artificial dynamics of the domain wall, typically for the Ising-like model. Both the QEW equation and the Ising model with the Monte Carlo method  cannot describe the orientations of the magnetic domains, the vortex structure, the spin-wave propagation, and  the interactions of the spin wave with the domain-wall motion, as well as the Walker breakdown etc. To describe  the domain-wall motion of real magnetic materials, in general, the Heisenberg-like model with a realistic  dynamic equation is required. For a comprehensive understanding of the depinning phase transition in the experiments, the fundamental approach based on the LLG equation looks promising [12][13][14], but up to now the progress is much limited [22][23][24]. The depinning phase transition is a dynamic phase transition, and usually of second order. The numerical simulation in the stationary state close to the transition point is very difficult due to critical slowing down. Recently, the transition point in a magnetic thin film was estimated, with the static exponent β as input [23]. On the other hand, the non-stationary dynamic approach looks efficient in tackling the dynamic phase transitions, since the measurements are carried out in the short-time regime of the dynamic evolution [18,19]. The domain-wall motion in a magnetic nanowire was analyzed based on the LLG equation, and a dynamic scaling behavior was observed [24]. However, the data collapse was performed in a rather wide range of the external field, much beyond the standard critical regime. The power-law behavior of the order parameter, which should be the crucial signal for a second-order phase transition, has not been detected. More importantly, most relevant experiments of depinning phase transitions were performed with two-dimensional thin films [12][13][14]. It is still unclear whether there exists a real depinning phase transition in the magnetic nanowire. Therefore, further numerical study along this direction is important and necessary.
In this paper, numerical simulations based on the LLG equation with the Heisenberg-like model are performed for the nonstationary domain-wall motion in a magnetic thin film with anisotropy. It is demonstrated that the transition point, and both the static and dynamic critical exponents can be accurately determined, with the refined short-time dynamic approach. We choose the model parameters such that the effects of the spin waves are negligible, thus we may focus on the intrinsic domain-wall motion [9, 12-14].

LLG equation and scaling analysis
The dynamic evolution of a magnetic system at zero temperature is described by the LLG equation, which is more convenient in numerical simulations. We consider the Heisenberg model with anisotropy, and the dimensionless effective field is given by Here j(i) denotes the nearest neighbors of site i, A, K P and K ⊥ are the exchange, easy-anisotropy and hardanisotropy coefficients, respectively, H is a constant external driving field, and the internal random field h i  represents quenched disorder, which could be a consequence of vacancies, defects, impurities and dislocations in the material. The direction of the random field h i  is randomly set on the unit sphere, while the amplitude of the random field, i.e. h h Quenched disorder in other forms such as the random bond usually yields similar results [19,21].
All parameters in our model are set to be dimensionless. For this purpose, the time, length and energy density are measured in units of , respectively, where γ is the gyromagnetic ratio. The exchange coefficient A and anisotropy coefficients K P and K ⊥ in equation , respectively, where d is the lattice constant. In our simulations, we take the parameters, A=1, α=1, K P =1, and K 10 = . The experiment value of the damping coefficient α varies in a wide regime, from α=0.001 to 4.0 in different material [9,12,23,24]. The smaller α is, the faster the domain wall moves, so a larger lattice must be applied. We set α=1, which allows the simulations with a reasonable lattice size. Extra simulations show that at least in the range of α=0.1 to 1.0 the phase transition remains robust. In addition, α=1 together with K 10 = ensures that the transition field is below the Walker threshold. According to [5,9], for example, H K 2 W a =^. In our simulations, H W = 5, while the transition field H c is about ten times smaller.
Simulations are performed on a rectangular lattice with lattice sizes L y = 512 and L z = 64, up to the maximum time t max =2560. We adopt the periodic and anti-periodic boundary conditions for the y and z axes, respectively. In our simulations, the domain wall does not reach the boundary along the z axis, and there is also no finite size effect along the y axis. Total samples of the quenched disorder for average are about 24 000.
As shown in figure 1(a), the domain wall is initially set at the position z=10, with two domains of opposite orientations on both sides [25]. The simulation starts by switching on a driving field along the positive direction of the easy-z axis. During the dynamic evolution, the domain wall propagates and fluctuates. Now denoting m i  by m y z t , ,  ( ), we define the height function of the domain wall [15,18,26], We then calculate the average velocity of domain wall, where á ñ  represents both the statistical average and the average along the y direction. We should note that the statistical average here is a real ensemble average. This is different from the simulation in the stationary state, where the ensemble average is usually replaced by the time average. In our simulations, different samples for the average are realized by different microscopic configurations of the random fields.
In the stationary state, the domain-wall motion driven by a constant external field exhibits a depinning phase transition. The transition field H c separates the pinning state and the sliding state of the domain wall, with the average velocity υ=0 and υ>0, respectively. Since the depinning phase transition is of second order, there exists critical slowing down, i.e., the correlation time is divergent. Numerical simulations of the depinning phase transition in the stationary state is very difficult.
In the past years, many activities have been devoted to the critical dynamics around a continuous phase transition, and it has been confirmed that there emerges already a universal dynamic scaling form even when the dynamic system is far from the equilibrium or in a non-stationary state. Based on the dynamic scaling form, numerical methods have been developed to detect the transition point and measure both the static and dynamic critical exponents [27,28]. Since the measurements are now performed in the short-time regime, one does not suffer from critical slowing down. Recent numerical simulations show that the non-stationary dynamic approach looks also rather efficient in tackling the depinning phase transitions of the domain-wall motion [18,19].
Assuming the depinning phase transition to be of second order, therefore, there should exist a dynamic scaling form in the macroscopic short-time regime, after a microscopic time scale t mic [27][28][29]. In the critical regime and for a sufficiently large lattice, the dynamic scaling form of the order parameter, i.e., the average velocity υ(t) of the domain wall, is described as is the reduced field, β and ν are static critical exponents, and ξ(t) is the spatial correlation length. At the critical point 0 t = , it leads to a power-law behavior, Therefore, searching for the best power-law behavior of t, u t ( ), one may locate the transition field H c , and then measure the exponent β/ν. The exponent 1/ν can be extracted from the logarithmic derivative of equation (6), Further, the roughness function is defined as , . 9 A more informative quantity is the height correlation function , , which describes the spatial correlation in the y direction. To estimate the spatial correlation length, we introduce where M(t) is the global magnetization, and M (2) (t) is its second moment. At the critical point, when the correlation length ξ(t) = L y , the roughness function ω 2 (t) and correlation function C(r, t) should scale as Here ζ is the global roughness exponent, and ζ loc is the local one. Since ω 2 (t) describes the fluctuation in the z direction, and M t M t includes those in both z and y directions, F(t) is proportional to the correlation length ξ(t) [21], In simple cases, t t z 1 x( ) with z being the dynamic exponent. Practically there may exist corrections to scaling. mic~. Although the LLG equation of the Heisenberg model with anisotropy is much more complicated than the QEW equation and the Monte Carlo simulation of the Ising model, the nonstationary dynamic approach to the depinning phase transition still shows its very efficiency. The transition point, and both the static and dynamic critical exponents are rather accurately determined. We should emphasize that in the present refined dynamic approach, the measurement of the dynamic exponent z does not affect those of the transition point and all other critical exponents [30]. This is particularly crucial when there exists a strong correction to scaling for ξ(t).
The large easy-anisotropy coefficient K P seems restricting the configuration space of the local magnetization M i  to two small regions, similar to that of the Ising spin. However, the depinning phase transition of the domainwall motion is a dynamic one, and its universality class is closely associated with the equation of motion.  (6), respectively. This should lead to a different understanding of the experimental data [13], e.g., where β=1/4 is assumed following the QEW equation. In addition, in a recent numerical simulation with the LLG equation [23], β=1/3 is taken as input to estimate the transition point H c . A larger exponent β should reflect that there exist more and stronger moving mechanisms of the domain wall in the Heisenberg model with the LLG equation. On the other hand, the universality class of a continuous phase transition is spatial dimension dependent. In [24], the depinning transition of the one- . Dashed lines slightly offset from the data points in (b) and in the inset of (a) just show the slopes of the curves without considering corrections to scaling. dimensional nanowire was simulated, and one obtained the critical exponent β=0.30 (1). Compared with this result, β=0.59(1) for the two-dimensional film looks also reasonable. The static exponent ν=1.10(3) of the LLG equation is about 20% smaller than those of the QEW equation and the Monte Carlo simulation. The dynamic exponent z=1.41(2) of the LLG equation happens to be the same as that of the QEW equation within errors. On the other hand, the tilt symmetry characterized by the scaling law (2−ζ)ν=1 is an interesting property of the QEW equation, which relates the roughness length scale and the correlation length scale at the depinning phase transition. This tilt symmetry is violated in the Monte Carlo simulation of the Ising model, possibly due to the overhangs and islands [21]. However, (2−ζ) ν=0.98 (3) indicates that the tilt symmetry holds for the LLG equation of the Heisenberg model with anisotropy. Actually, the dynamic effect of the overhangs and islands looks not very strong for the LLG equation.
The QEW equation is only a phenomenological equation for an elastic string. The Monte Carlo simulation of the Ising model with quenched disorder is an approximate method for a partial understanding of the dynamic properties of magnetic materials. The LLG equation for the Heisenberg-like model comprehensively describes the dynamic properties of magnetic materials.

Conclusion
Based on the LLG equation with the Heisenberg-like model, we numerically simulate the nonstationary dynamic behavior of the domain wall in a magnetic thin film with anisotropy, and confidently detect a depinning phase transition induced by quenched disorder, which is of second order. Power-law behaviors of the domain-wall velocity and relevant observables at the critical point are clearly observed. Thus the transition field, static and dynamic exponents are accurately determined, with the dynamic scaling behavior far from stationary. The results demonstrate the very efficiency of the nonstationary dynamic approach to the LLG equation, and show that the universality class of the LLG equation of the anisotropic Heisenberg model is different from that of the Monte Carlo dynamics of the Ising model, and also the QEW equation. The method may be further applied to understand the complex dynamic structure of the domain wall, such as the vortex domain walls, different orientations of the magnetic domains, and the interactions of the spin wave with the domain-wall motion.