Dynamical depinning of chiral domain walls

The domain wall depinning field represents the minimum magnetic field needed to move a domain wall, typically pinned by samples' disorder or patterned constrictions. Conventionally, such field is considered independent on the Gilbert damping since it is assumed to be the field at which the Zeeman energy equals the pinning energy barrier (both damping independent). Here, we analyse numerically the domain wall depinning field as function of the Gilbert damping in a system with perpendicular magnetic anisotropy and Dzyaloshinskii-Moriya interaction. Contrary to expectations, we find that the depinning field depends on the Gilbert damping and that it strongly decreases for small damping parameters. We explain this dependence with a simple one-dimensional model and we show that the reduction of the depinning field is related to the internal domain wall dynamics, proportional to the Dzyaloshinskii-Moriya interaction, and the finite size of the pinning barriers.


I. INTRODUCTION
Magnetic domain wall (DW) motion along ferromagnetic (FM) nanostructures has been the subject of intense research over the last decade owing to its potential for new promising technological applications 1,2 and for the very rich physics involved. A considerable effort is now focused on DW dynamics in systems with perpendicular magnetic anisotropy (PMA) which present narrower DWs and a better scalability. Typical PMA systems consist of ultrathin multi-layers of heavy metal/FM/metal oxide (or heavy metal), such as Pt/Co/Pt 3,4 or Pt/Co/AlOx [5][6][7] , where the FM layer has a thickness of typically 0.6 − 1 nm. In these systems, PMA arises mainly from interfacial interactions between the FM layer and the neighbouring layers (see Ref. 8 and references therein). Another important interfacial effect is the Dzyaloshinskii-Moriya interaction (DMI) 9,10 , present in systems with broken inversion symmetry such as Pt/Co/AlOx. This effect gives rise to an internal inplane field that fixes the DW chirality (the magnetization rotates always in the same direction when passing from up to down and from down to up domains) and it can lead to a considerably faster domain wall motion 10 and to new magnetic patterns such as skyrmions 11 or helices 12 . Normally, DWs are pinned by samples' intrinsic disorder and a minimum propagation field is needed in order to overcome such pinning energy barrier and move the DW. Such field is the DW depinning field (H dep ) and it represents an important parameter from a technological point of view since a low depinning field implies less energy required to move the DW and, therefore, a energetically cheaper device.
From a theoretical point of view, DW motion can be described by the Landau-Lifshitz-Gilbert (LLG) equation 13 which predicts, for a perfect sample without disorder, the velocity vs field curve depicted in Fig. 1 and labelled as Perfect. In a disordered system, experiments have shown that a DW moves as a general onedimensional (1D) elastic interface in a two-dimensional disordered medium 3,4 and that it follows a theoretical velocity vs driving force curve, predicted for such interfaces 14,15 (also shown in Fig. 1 for T = 0 and T = 300K). Moreover, this behaviour can be reproduced by including disorder in the LLG equation [16][17][18] . At zero temperature (T = 0) the DW does not move as long as the applied field is lower than H dep , while, at T = 0, thermal activation leads to DW motion even if H < H dep (the so called creep regime). For high fields (H >> H dep ) the DW moves as predicted by the LLG equation in a perfect system. Within the creep theory, the DW is considered as a simple elastic interface and all its internal dynamics are neglected. Conventionally, H dep is considered independent of the Gilbert damping because it is assumed to be the field at which the Zeeman energy equals the pinning energy barrier 19,20 (both damping independent). Such assumption, consistently with the creep theory, neglects any effects related to the internal DW dynamics such as DW spins precession or vertical Bloch lines (VBL) formation 21 . The damping parameter, for its part, represents another important parameter, which controls the energy dissipation and affects the DW velocity and Walker Breakdown 22 . It can be modified by doping the sample 23 or by a proper interface choice as a consequence of spin-pumping mechanism 24 . Modifications of the DW depinning field related to changes in the damping parameter were already observed in in-plane systems 23,25 and attributed to a non-rigid DW motion 23,25 . Oscillations of the DW depinning field due to the internal DW dynamics were also experimentally observed in in-plane similar systems 26 . Additional dynamical effects in soft samples, such as DW boosts in current induced motion, were numerically predicted and explained in terms of DW internal dynamics and DW transformations 27,28 .
Here, we numerically analyse the DW depinning field in a system with PMA and DMI as function of the Gilbert damping. We observe a reduction of H dep for low damping and we explain this behaviour by adopting a simple 1D model. We show that the effect is due to the finite size of pinning barriers and to the DW internal dynamics, related to the DMI and shape anisotropy fields. This article is structured as follows: in Section II we present the simulations method, the disorder implementation and the H dep calculations. The main results are outlined and discussed in Section III, where we also present the 1D model. Finally, the main conclusions of our work are summarized in Section IV.
DW velocity vs applied field as predicted by the LLG equation in a perfect system and by the creep law at T = 0 and T = 300K.

II. MICROMAGNETIC SIMULATIONS
We consider a sample of dimensions (1024 × 1024 × 0.6) nm 3 with periodic boundary conditions along the y direction, in order to simulate an extended thin film. Magnetization dynamics is analysed by means of the LLG equation 13 : where m(r, t) = M(r, t)/M s is the normalized magnetization vector, with M s being the saturation magnetization. γ 0 is the gyromagnetic ratio and α is the Gilbert damping. H eff = H exch + H DMI + H an + H dmg + H zûz is the effective field, including the exchange, DMI, uniaxial anisotropy, demagnetizing and external field contributions 13 respectively. Typical PMA samples parameters are considered: A = 17 × 10 −12 J/m, M s = 1.03 × 10 6 A/m, K u = 1.3 × 10 6 J/m 3 and D = 0.9 mJ/m 2 , where A is the exchange constant, D is the DMI constant and K u is the uniaxial anisotropy constant. Disorder is taken into account by dividing the sample into grains by Voronoi tessellation 29,30 , as shown in Fig. 2(a). In each grain the micromagnetic parameters {M s , D c , K u } change in a correlated way in order to mimic a normally distributed thickness 31 : where the subscript G stands for grain, t 0 is the average thickness (t 0 = 0.6nm) and δ is the standard deviation of the thickness normal distribution. The sample is discretized in cells of dimensions (2 × 2 × 0.6)nm 3 , smaller than the exchange length l ex ∼ 5nm. Grain size is GS=15 nm, reasonable for these materials, while the thickness fluctuation is δ = 7%. Eq. (1) is solved by the finite difference solver MuMax 3.9.3 29 .
A DW is placed and relaxed at the center of the sample as depicted in Fig. 2(b). H dep is calculated by applying a sequence of fields and running the simulation, for each field, until the DW is expelled from the sample, or until the system has reached an equilibrium state (i.e. the DW remains pinned): τ max < (α). τ max indicates the maximum torque, which rapidly decreases when the system is at equilibrium. It only depends on the system parameters and damping. For each value of α, we choose a specific threshold, (α), in order to be sure that we reached an equilibrium state (see Supplementary Material 32 for more details). The simulations are repeated for 20 different disorder realizations. Within this approach, H dep corresponds to the minimum field needed to let the DW propagate freely through the whole sample. In order to avoid boundaries effects, the threshold for complete depinning is set to m z > 0.8, where m z is averaged over all the realizations, i.e. m z = N i=1 m z i /N , where N = 20 is the number of realizations. We checked that, in our case, this definition of H dep coincides with taking H dep = Max{H i dep }, with H i dep being the depinning field of the single realization. In other words, H dep corresponds to the minimum field needed to depin the DW from any possible pinning site considered in the 20 realizations 33 .
Following this strategy, the DW depinning field is numerically computed with two different approaches: (1) by Static simulations, which neglect any precessional dynamics by solving This is commonly done when one looks for a minimum of the system energy and it corresponds to the picture in which H dep simply depends on the balance between Zeeman and pinning energies. 34 (2) by Dynamic simulations, which include precessional dynamics by solving the full Eq. (1). This latter method corresponds to the most realistic case. Another way to estimate the depinning field is to calculate the DW velocity vs field curve at T = 0 and look for minimum field at which the DW velocity is different from zero. For these simulations we use a moving computational region and we run the simulations for t = 80ns (checking that longer simulations do not change the DW velocity, meaning that we reached a stationary state). This second setup requires more time and the calculations are repeated for only 3 disorder realizations.
Using these methods, the depinning field H dep is calculated for different damping parameters α.

A. Granular system
Our first result is shown in Fig. 3(a)-(b), which depicts the final average magnetization m z as function of the applied field for different damping parameters. In the Static simulations ( Fig. 3(a)) H dep does not depend on damping, so that a static depinning field can be defined. Conversely, in the Dynamic simulations ( Fig. 3(b)), H dep decreases for low damping parameters. The depinning field is indicated by a star in each plot and the static depinning field is labelled as H s . The same result is obtained by calculating H dep from the DW velocity vs applied field plot, shown in Fig. 3(c). The stars in Fig. 3(c) correspond to the depinning fields calculated in the previous simulations and they are in good agreement with the values predicted by the velocity vs field curve. The dynamical depinning field µ 0 H d , normalized to the static depinning field µ 0 H s = (87 ± 1)mT, with µ 0 being the vacuum permeability, is shown in Fig. 3(d) as function of the damping parameter α. H d saturates for high damping (in this case α ≥ 0.5) while it decreases for low damping until H d /H s ∼ 0.4 at α = 0.02. This reduction must be related to the precessional term, neglected in the static simulations. The same behaviour is observed with different grain sizes (GS=5 and 30 nm) and with a different disorder model, consisting of a simple variation of the K u module in different grains. This means that the effect is not related to the grains size or to the particular disorder model we used.
Additionally, Fig. 4 represents the DW energy 35 as function of DW position and damping parameter for µ 0 H z = 70 mT. At high damping, the average DW energy density converges to σ ∞ ∼ 10 mJ/m 2 , in good agreement with the analytical value σ 0 = 4 On the contrary, for low damping, the DW energy increases up to σ(0.02) ∼ 14 mJ/m 2 . This increase, related to DW precessional dynamics, reduces the effective energy barrier and helps the DW to over- come the pinning barriers. Fig. 4(c) shows the total energy of the system (including Zeeman). As expected 36 , the energy decreases as the DW moves.
Finally, Fig. 5 shows the DW motion as function of time for α = 0.02 and α = 0.5, along the same grain pattern (and therefore along the same pinning barriers). The applied field is µ 0 H z = 70mT, which satisfies H d (0.02) < H z < H d (0.5). The initial DW configuration is the same but, for α = 0.02, VBL start to nucleate and the DW motion is much more turbulent (see Supplementary Material 32 for a movie of this process). At t = 4 ns the DW has reached an equilibrium position for α = 0.5, while it has passed through the (same) pinning barriers for α = 0.02. Thus, one might think that the reduction of the depinning field could be related to the presence of VBL and their complex dynamics 21 . Further insights about this mechanism are given by analysing the DW depinning at a single energy barrier as described in the next subsection.

B. Single barrier
In order to understand how the DW precessional dynamics reduces H dep , we micromagnetically analysed the DW depinning from a single barrier as sketched in Fig. 6. We considered a strip of dimensions (1024×256×0.6)nm 3 and we divided the strip into two regions, R 1 and R 2 , which are assumed to have a thickness of t 1 = 0.58 and t 2 = 0.62 nm respectively. Their parameters vary accordingly (see Sec. II), generating the DW energy barrier (δσ) shown in Fig. 6(b). A DW is placed and relaxed just before the barrier. The finite size of the DW (π∆ DW ∼ 15 nm, with ∆ DW being the DW width parameter) smooths the abrupt energy step and, in fact, the energy profile can be successfully fitted by using the Bloch profile 22 where x 0 = 20 nm is the step position, while σ 0 and σ 1 are the DW energies at the left and right side of the barrier as represented in Fig. 6(b). This means that the pinning energy barrier has a spatial extension which is comparable to the DW width. By performing the same static and dynamic simulations, we obtain a static depinning field of µ 0 H s = 120 mT and, when decreasing the damping parameter, we observe the same reduction of the depinning field as in the granular system (see Fig. 6(c)). In this case the DW behaves like a rigid object whose spins precess coherently and no VBL nucleation is observed. Hence, H dep reduction does not depend directly on the presence of VBL but on the more general mechanism of spins' precession already present in this simplified case. Nevertheless, an important characteristic of these single barrier simulations is that the barrier is localized and it has a finite size which is of the order of the DW width. Note that the same holds for the granular system: despite a more complex barrier structure, the dimension of the single barrier between two grains has the size of the DW width.
Thus, in order to understand the interplay between the DW precessional dynamics and the finite size of the barrier, we considered a 1D collective-coordinate model with a localized barrier. The 1D model equations, describing the dynamics of the DW position q and the internal angle φ (sketched in Fig. 2(c)), are given by 16 where H K = M s N x is the shape anisotropy field, favouring Bloch walls, with N x = t 0 log 2/(π∆ DW ) 37 being the DW demagnetizing factor along the x axis. H DMI = D/(µ 0 M s ∆ DW ) is the DMI field. H int (φ) represents the internal DW field, which includes DMI and shape anisotropy. H int favours Bloch (φ = ±π/2) or Néel wall (φ = 0 or φ = π) depending on the relative strength of H K and H DMI . In our system, the DMI dominates over shape anisotropy since µ 0 H DMI ∼ 170 mT while µ 0 H K ∼ 30 mT. Hence, the DW equilibrium angle is  φ = π (φ = 0 or φ = π additionally depends on the sign of the DMI). H p (q) is the DW pinning field, obtained from the DW energy profile (Eq. (4)) as follows: the maximum pinning field is taken from the static simulations while the shape of the barrier is taken as the normalized DW energy gradient (see Supplementary Material 32 for more details), The corresponding pinning field is plotted in Fig. 7(a). 38 The results for the dynamical H dep , obtained with this modified 1D model, are plotted in Fig. 6(c) and they show a remarkable agreement with the single barrier micromagnetic simulations. This indicates that the main factors responsible for the reduction of H dep are already included in this simple 1D model. Therefore, additional insights might come from analysing the DW dynamics within this 1D model. Fig. 7(b) and (c) represents the DW internal angle φ and the DW position q as function of time for different damping. The plots are calculated with µ 0 H z = 55 mT which satisfies H dep (0.02) < H z < H dep (0.1) < H dep (0.5). As shown in Fig. 7(b) and (c), below the depinning field (α = 0.1, α = 0.5), both the internal angle and the DW position oscillate before reaching the same final equilibrium state. However, the am-plitude of these oscillations (the maximum displacement) depends on the damping parameter. Fig. 7(d) shows the final equilibrium position as function of the applied field for different damping. The equilibrium position is the same for all damping and it coincides with the position at which H z = H p (q). Conversely, the maximum displacement, shown in Fig. 7(e), strongly increases for low damping parameters. For applied field slightly smaller than the depinning field, the DW reaches the boundary of the pinning barrier, meaning that a further increase of the field is enough to have a maximum displacement higher than the barrier size and depin the DW. In other words, the decrease of the depinning field, observed in the single barrier simulations, is due to DW oscillations that depend on α and that can be larger than the barrier size, leading to DW depinning for lower field. The DW dynamics and the depinning mechanism are further clarified in Fig. 7(f) and Fig. 7(g). Fig. 7(f) represents the DW coordinates {q, φ} for µ 0 H z = 55 mT and different damping. Before reaching the common equilibrium state, the DW moves in orbits (in the {q, φ} space) whose radius depends on the damping parameter. For α = 0.5 (black line) the DW rapidly collapse into the final equilibrium state. Conversely, for α = 0.1 (red open circles), the DW orbits around the equilibrium state before reaching it. If the radius of the orbit is larger than the barrier size the DW gets depinned, as in the case of α = 0.02 (blue full circles). This mechanism is also represented in Fig. 7(g), where the DW orbits are placed in the energy landscape. The energy is calculated as where σ DW is given by Eq. (4). Fig. 7(g) shows that the equilibrium state corresponds to the new minimum of the energy landscape. Furthermore, it confirms that the applied field is below the static depinning field, at which the pinning barrier would have been completely lifted. Nevertheless, while reaching the equilibrium state, the DW moves inside the energy potential and, if the radius of the orbit is larger than the barrier size, the DW can overcome the pinning barrier, as shown for α = 0.02 in Fig. 7(g). At this point we need to understand why the amplitude of the DW oscillations depends on damping. By solving Eq. (5) and Eq.(6) for the equilibrium state (q = 0,φ = 0) we obtaiṅ since µ 0 H DMI µ 0 H K and, therefore, H int ≈ −(π/2)H DMI sin φ. These equations have a single common solution which corresponds to |H p (q)| = H z and φ = φ 0 = π (at which H int (π) = 0). However, at t = 0, the DW starts precessing under the effect of the applied field and, if φ = π when |H p (q)| = H z , the DW does not stop at the final equilibrium position but it continues its motion, as imposed by Eq. (8) and (9). In other words, the DW oscillations in Fig. 7(b) are given by oscillations of the DW internal angle φ, around its equilibrium value φ 0 = π. These oscillations lead to a modification of the DW equilibrium position due to the DW internal field (H int (φ)), which exerts an additional torque on the DW in order to restore the equilibrium angle. As previously commented, if the amplitude of these oscillations is large enough, the DW gets depinned. From Eq. (8) we see that the new equilibrium position (and therefore the amplitude of the oscillations) depends on the DMI field, the value of the DW angle φ and the damping parameter.
In particular, damping has a twofold influence on this dynamics: one the one hand, it appears directly in Eq. (8), dividing the internal field, meaning that for the same deviation of φ from equilibrium, we have a stronger internal field for smaller damping. On the other hand, the second influence of damping is on the DW internal angle: once the DW angle has deviated from equilibrium, the restoring torque due to DMI is proportional to the damping parameter (see Eq. (9)). Hence, a lower damping leads to lower restoring torque and a larger deviation of φ from equilibrium. The maximum deviation of φ from equilibrium (δφ = φ max − φ 0 ) is plotted in Fig. 8(b) as function of damping for µ 0 H z = 40 mT. As expected, a lower damping leads to a larger deviation δφ.
In this latter section, the DW was set at rest close to the barrier and, therefore, the initial DW velocity is zero. Nevertheless, one might wonder what happens when the DW reaches the barrier with a finite velocity. We simulated this case by placing the DW at an initial distance d 1 = 200 nm from the barrier. The depinning is further reduced in this case (see Supplementary Material 32 for more details). However, in the static simulations, the depinning field remains constant, independently from the velocity at which the DW reaches the barrier, meaning that the reduction of H dep is again related to the DW precession. When the DW starts from d 1 it reaches the barrier precessing, thus with a higher deviation from its equilibrium angle, leading to a higher effect of the internal field.

C. Different DMI and pinning barriers
Finally, by using the 1D model it is possible to explore the dependence of H dep on the pinning potential amplitude H s (related to the disorder strength) and on the DMI constant D. The depinning field as function of damping for different values of H s is plotted in Fig. 9(a). The reduction of H dep is enhanced for larger values of H s (strong disorder). This is consistent with our explanation, since for strong disorder we need to apply larger fields that lead to larger oscillations of φ. Fig. 9(b) represents the dynamical H dep as function of  damping for µ 0 H s = 120 mT and different DMI constants (expressed in term of the critical DMI constant D c = 4 √ AK 0 /π = 3.9 mJ/m 2 ) 39 . In this case, the reduction of H dep is enhanced for low DMI, until D = 0.05D c , but a negligible reduction is observed for D = 0. This non-monotonic behaviour can be explained by looking at the dependence of δφ and H int on the DMI constant. Fig. 10(a) shows the maximum fluctuation δφ as function of DMI for µ 0 H z = 30 mT. δφ increases for low DMI and it has a maximum at πH DMI = H K , which in our case corresponds to D = 0.014D c . The increase of δφ for small values of D is due to the smaller restoring torque in Eq. (9). This holds until πH DMI = H K , where shape anisotropy and DMI are comparable and they both affect the DW equilibrium configuration. As a consequence, the reduction of H dep is enhanced by de- Another contribution is given by the amplitude of the internal field, H int . Fig. 10(b) depicts µ 0 H int as function of δφ and D. The maximum δφ, obtained at µ 0 H z = 30 mT, is additionally marked in the plot. The internal field decreases with the DMI but this reduction is compensated by an increase in δφ, which leads to an overall increase of µ 0 H int , as discussed in the previous part. However, at very low DMI, the internal field is dominated by shape anisotropy and, independently on the DW angle displacement, it is too small to have an effect on the depinning mechanism. Note, however, that the amplitude of H int should be compared with the amplitude of the pinning barrier H s . Fig. 9(b) is calculated with µ 0 H s = 120 mT and the internal field, given by shape anisotropy (H K /2 ∼ 15 mT), has indeed a negligible effect. However, larger effects are observed, in the case D = 0, for smaller H s , with reduction of H dep up to H d /H s ∼ 0.6, as shown in Fig. 9(c), which is calculated with µ 0 H s = 30 mT. In other words, the reduction of the depinning field depends on the ratio between the pinning barrier and the internal DW field.
Finally, it is interesting to see what happens for weaker disorder and different DMI in the system with grains. Fig. 11 shows the dynamical H dep , for different pinning potential and different DMI, obtained in the granular system. The results are in good agreement with what predicted by the 1D model for different disorder strengths. However, we observe a smaller dependence on the DMI parameter. This is due to two reasons: (1) in the system with grains the static pinning barrier is µ 0 H s = 87 mT and the dependence of the depinning field with DMI is smaller for smaller barriers, as shown in Fig. 9(c). (2) The DW motion in the granular system presents the formation of VBL which might also contribute to the reduction of the depinning field. The mechanism is the same: a VBL is a non-equilibrium configuration for the DW (as a deviation of φ from equilibrium) that generates additional torques on the DW, which contribute to the DW depinning.

IV. CONCLUSIONS
To summarize, we have analysed the DW depinning field in a PMA sample with DMI and we found that H dep decreases with the damping parameter with reductions up to 50%. This decrease is related to the DW internal dynamics and the finite size of the barrier: due to DW precession, the DW internal angle (φ) deviates from equilibrium and triggers the internal DW field (DMI and shape anisotropy) which tries to restore its original value. At the same time, the internal field pushes the DW above its equilibrium position within the energy barrier. This mechanism leads to DW oscillations and, if the amplitude of the oscillations is higher than the barrier size, the DW gets depinned for a lower field. Deviations of φ from equilibrium and DW oscillations are both damping dependent and they are enhanced at low damping.
In the system with grains the mechanism is the same but deviations from the internal DW equilibrium include the formation of VBL with more complex dynamics. The effect is enhanced for low DMI (providing that πH DMI > H K ) and for stronger disorder since we need to apply larger external fields, which lead to larger DW oscillations. These results are relevant both from a technological and theoretical point of view, since they firstly suggest that a low damping parameter can lead to a lower H dep . Furthermore, they show that micromagnetic calculations of the depinning field, neglecting the DW precessional dynamics can provide only an upper limit for H dep , which could actually be lower due to the DW precessional dynamics. creasing the sample dimension along the x direction, we increase the probability of finding the highest possible hj in the single realization and the average of H i dep will tend to the maximum. 34 This is solved by the Relax solver of MuMax with the assumption α/(1 + α 2 ) = 1. 35 The DW energy is calculated as the energy of the system with the DW minus the energy of the system without the DW (uniform state). The profile is obtained by moving the DW with an external applied field and then subtracting the Zeeman energy. (1998). 38 The same results are obtained with a Gaussian barrier, meaning that the key point is the finite size of the barrier rather than its shape. 39 For D > Dc, DW have negative energies and the systems spontaneously breaks into non-uniform spin textures.

Appendix A: Maximum torque and equilibrium state
In this section we show in more detail how the maximum torque represents an indicator of the equilibrium state. Maximu torque is defined as over all cells with label i = {1, ..., N = N x · N y }. MuMax3.9.3 29 can provide this output automatically if selected. We perform the same simulations as indicated in the main text, without any stopping condition, but simply running for t = 20 ns. Fig. 12(a) shows the average m z component for α = 0.2 and B z = 10 mT, while Fig. 12(b) depicts the corresponding maximum torque. We can see that, once the system has reached equilibrium, the maximum torque has dropped to a minimum value. The same results is obtained for different damping but the final maximum torque is different. Numerically this value is never zero since it is limited by the code numerical precision and by the system parameters, in particular by damping. Fig. 12(c) represents the maximum torque as function of applied field for different damping. The maximum torque is clearly independent on the applied field but depends on the damping value. Finally, Fig. 12(d) shows the max torque as function of damping. The maximum torque decreases with damping and it saturates for α ≥ 0.5 since we have reached the minimum numerical precision of the code 29 . For higher damping the maximum torque oscillates around this minimum sensibility value, as shown in the inset of Fig. 12(d). The value obtained with these preliminary simulations is used to set a threshold (α) for the depinning field simulations in order to identify when the system has reached an equilibrium. Furthermore, additional tests were performed, without putting any max torque condition, but simply running the simulations for a longer time (t = 80, 160 ns) and calculating the depinning field in order to ensure that the results obtained with these two method were consistent, i.e., that we have actually reached an equilibrium state with the maximum torque condition.