Barkhausen noise from formation of 360 ◦ domain walls in disordered permalloy thin ﬁlms

Barkhausen noise in disordered ferromagnets is typically understood to originate primarily from jerky ﬁeld-driven motion of domain walls. We study the magnetization reversal process in disordered permalloy thin ﬁlms using micromagnetic simulations, and ﬁnd that the magnetization reversal process consists of gradual formation of immobile 360 ◦ domain walls via a sequence of localized magnetization rotation events. The density of 360 ◦ domain walls formed within the sample as well as the statistical properties of the Barkhausen jumps are controlled by the disorder strength.

In this work, we employ full micromagnetic simulations of disordered Py thin films to properly capture the details of the field-driven magnetization reversal process and the related statistical properties, including the emergence of BN.Unlike simple Ising-type models, micromagnetic simulations provide a proper description of the full vectorial nature of the order parameter (magnetization) and include all the relevant energy terms (such as the demagnetizing energy).Contrary to the paradigm of BN due to jerky motion of DWs, we find that field-driven BN in Py thin films is due to the gradual formation of largely immobile 360 • DWs [19][20][21][22][23][24][25][26][27][28][29][30][31] via a sequence of localized and intermittent magnetization rotation events.Once formed, the 360 • DWs get progressively narrower due to an increasing driving field induced effective anisotropy, before disappearing as the oppositely saturated state is reached.Both the emerging configuration of 360 • DWs (e.g., their density) and the statistics of BN are found to be dependent on disorder strength, suggesting that the magnetization reversal process in disordered Py thin films is governed by disorder-induced criticality [10], but with an important difference compared to Ising models: no infinite, spanning avalanche (i.e., an avalanche spanning the finite system along at least one of its dimen- sions, resulting in a jump of magnetization in an infinite system below the critical disorder) is observed for weak, subcritical disorder.The present work therefore focuses on investigating how the Barkhausen effect and the related jump size distributions emerge from the previously unidentified mechanism of gradual formation of immobile 360 • domain walls in thin films with negligible magnetocrystalline anisotropy, as opposed to the well-known mechanism related to the jerky motion of domain walls.

II. COMPUTATIONAL METHODS
Our micromagnetic simulations are performed using the GPU-accelerated simulation code mumax3 [32].Typically, micromagnetic simulations rely on the dynamical approach where the magnetization dynamics is governed by the Landau-Lifshitz-Gilbert equation, written as where M = M s m is the magnetization vector, M s is the saturation magnetization, m = m x x + m y ŷ + m z ẑ is the unit vector pointing along the magnetization, γ is the gyromagnetic ratio of electron, α is the Gilbert damping constant, and H eff is the effective field.The total energy of the system (and hence H eff ) consists of contributions from exchange, demagnetization, and Zeeman terms as arXiv:2209.10847v2[cond-mat.mtrl-sci]6 Mar 2023 described in detail in the Supplemental Material [33].Due to the vanishing magnetocrystalline anisotropy of Py, the anisotropy energy is neglected here, that is the magnetocrystalline anisotropy constants are set to zero.
An example of a simulated system is shown in Fig. 1(a).We discretize the square thin film on a finite difference grid with spacing of 4 nm in the x and y in-plane directions where the linear system size is 4096 nm and consider a single 16 nm long cell in the perpendicular z direction (the film thickness).Periodic boundary conditions are employed in the in-plane x and y directions.We use the common literature values for the exchange stiffness A = 13 pJ/m and the average saturation magnetization M s = 800 kA/m.Structural disorder is introduced via a spatially varying M s by performing a 2D Voronoi tessellation on each sample with the grain size of 30 nm (the disorder correlation length), and assigning a random value of M s in each grain from a normal distribution with mean of 800 kA/m and standard deviations (the disorder strengths) between 5 % and 35 % of the mean value (negative M s values are avoided by re-drawing until a positive value is obtained) [34,35].This way, we establish a random component in the energy landscape of the magnetization state, characterized by the correlation length (Voronoi grain size) and magnitude (disorder strength) of the disorder.In a material, such disorder can be caused by inhomogeneous distribution of elements in the sample, impurity material, defects, and other lattice imperfections.Introducing the disorder by parameter variation in the zero-anisotropy case can be done through variation of either M s or A or both.Here we investigate the effect of the M s variation only.According to our test runs, simultaneous variation of exchange stiffness A does not alter the magnetization dynamics qualitatively but the dynamical process is dominated by the spatial M s variation.
In each run, the magnitude B ext of the external field B ext = B ext x is swept from 100 mT (at which point the system is close to saturation with m x ≈ 1) to −100 mT.The run is considered finished once m x falls below −0.98 or the external field reaches −100 mT.This procedure thus produces half of the hysteresis loop where m x reverses from a value close to +1 to a value close to −1.
The simulation times of dynamic simulations are in practice limited to the microsecond range, hence seriously limiting the range of accessible field frequencies.Therefore, we mainly focus on quasistatic simulations, where the total energy of the system is consecutively minimized while altering B ext between each minimization in small steps ∆B ext .To validate the use of quasistatic simulations, we present in Fig. S1 (Supplemental Material [33]) results for both dynamical and quasistatic simulations, showing how the dynamic hysteresis loops approach the hysteresis loop obtained from the quasistatic simulations in the low frequency limit of the sinusoidal driving field.In what follows, we thus consider quasistatic simulations with |∆B ext | = 4 µT; see Fig. S2 [33] for the justification of the selected value.For an example of the half-loop pro-duced by this protocol, see Fig. 1(b).For each disorder strength, we collect statistics by running 40 simulation runs with different disorder realizations.

A. Hysteresis curves
Figure 2 shows all of the simulated hysteresis curves for each disorder strength.The curves for 5 and 10 % disorders are quite smooth, i.e., Barkhausen jumps are largely absent.The magnetization is being rotated in a smooth and continuous fashion as the applied field is driven from 100 mT to −100 mT, i.e., no infinite spanning avalanche takes place unlike in RFIM for weak disorder.Note also that the magnetization does not saturate to −1 before B ext = −30 mT.In contrast, the systems with the highest disorder strengths, 30 % and 35 %, exhibit bursty behaviour with more frequent and larger abrupt avalanches.The saturated state with m x = −1 is reached before B ext = −30 mT for most runs.The curves with 20 % disorder show features of both kinds: long smooth sections interrupted by occasional abrupt jumps.
Following this observation, let us make a distinction between the parts of the hysteresis curve where (1) the magnetization is rotated smoothly and continuously and (2) a large part of the magnetization is changed abruptly (Barkhausen jumps).We note that the tangential slopes of the smooth parts in between avalanches seem to have a characteristic magnitude.Also, the curves with 10 and 20 % disorder strengths seem to be divided into two separate subgroups, as well as the 30 % curves into four subgroups (Fig. 2).In the following, we will connect all of these features of the hysteresis curves to the disorder strength, the respective Barkhausen jump size distributions, and the magnetization reversal mechanisms, including the appearance of the 360 • DWs.

B. Barkhausen jump size distributions
Figure 3 shows the Barkhausen jump size distributions, that is the distributions of the absolute changes s ≡ |∆m x | of m x during each field step in the quasistatic simulations.Based on the renormalization group theory and previous studies of Barkhausen noise statistics [9], we expect that the tails of size distributions P (s) follow truncated power laws with a critical exponent τ , where g(x) is a cutoff scaling function, and s 0 a cutoff scale parameter.The tails of the distributions in Fig. 3 do behave according to Eq. ( 2) with a disorderdependent s 0 (see below for more details), but a bump exists in each distribution below the power-law part, introducing a characteristic length scale to the statistics.Increasing the disorder strength leads to translation of the bump to smaller jump sizes and possibly to broadening of the bump.We attribute the positions of the bumps to the largest magnitudes for the tangential slopes in the hysteresis curves where the magnetization reversal is smooth.Avalanches exceeding the characteristic size of the bumps, obeying Eq. ( 2), are taken to be the actual Barkhausen jumps.The disorder strength of 35 % seems to be close to the critical disorder strength since the tail of P (s) exhibits a clear power law behaviour.Even though the expected shape of Eq. ( 2) includes the cutoff function g(x), we employed a plain power law fit with g = 1 in the absence of properly resolved cutoffs due to limited statistics.Fitting a power law to the tail of P (s) for 35 % disorder results in τ = 1.75 ± 0.02.As illustrated in the inset of Fig. 3, the distributions for the other disorder strengths fall off from the power law trend in the order given by the disorder strength, i.e., s 0 in Eq. ( 2) is disorder-dependent.This trend is expected when approaching the critical disorder strength, and so we conclude that the critical disorder strength here is likely to be close to 35 %.In the literature, the reported values for the critical exponent τ for Py thin films are diverse: 1.65 [6], 1.33 [8], 1.45 [16], and 1.6 [36].Our estimate is close to those reported in Refs.[6] and [36].
Another interesting feature in the distributions is the almost perfect power law behaviour for jump sizes smaller than that of the bumps of the 5-20 % disorder distributions; this seems to be absent in the stronger disorder statistics.These very small magnetization changes correspond to tangential slopes of m x (B ext ) smaller in magnitude than the ones of the bump, and originate from the parts of the hysteresis curves where the magnetization is close to saturation.It is worth noting that the bump positions depend also somewhat on the field step used in the simulations; see Supplemental Material [33] and Fig. S2 therein.

C. Magnetization reversal processes
To understand the origin of the observed subgroups in the hysteresis curves we proceed to study the disorderdependent magnetization reversal mechanisms.For 5 % disorder (Fig. 4(a)), we observe that, interestingly, a structure of two 360 • DWs is formed at the late stages of the reversal process.This formation process proceeds via a gradual rotation of the magnetization in opposite directions within two "domains", such that the eventual 360 • DWs are formed via 90 • DWs at B ext = −1.59mT and 180 • DWs at B ext = -3.21mT.Due to topological protection [24], the 360 • DWs remain in the system until the end of the simulation, resulting in the absence of full negative saturation of m x for the B ext -values considered.Upon making B ext more negative, m x is slowly approaching −1 since the widths w of the DWs decrease with increasing the B ext (see Fig. 5 and the related dis- cussion below).Each of the runs with 5 % disorder exhibits the formation of two 360 • DWs, and the hysteresis curves of the different realizations overlap almost perfectly (Fig. 2).
For 10 and 20 % disorder, the two subgroups of the hysteresis curves in Fig. 2 correspond to two and four 360 • DWs per simulation box, respectively, and some of the 360 • DW structures disappear before the applied field reaches −25 mT; this is seen as m x reaching −1 in Fig. 2. The reversal mechanism in an example run where four 360 • DWs are formed is presented in Fig. 4(b).The mechanism for the formation of the DWs is similar to that of the 5 % disorder in Fig. 4(a), but there is more spatial variation in the local magnetization and the DWs appear more rough.As expected, even more roughness is present in the example cases with 30 and 35 % disor-ders shown in Figs.4(c) and (d), where the final magnetization structures are more complicated, and features resembling 360 • DWs can be recognized only in parts of the system.The magnetization configurations shown in both Figs.4(c) and (d) are fully reversed towards the field direction at the next field step.A closer look at the 30 % disorder strength hysteresis curves in Fig. 2 reveals that there are distinct subgroups of curves that correspond to varying number of DWs (up to 6 or 8) in the system.Four example runs with 30 % disorder are highlighted in Fig. S3  Thus, stronger disorder induces the formation of a higher density of 360 • DWs that also appear less stable -presumably disorder lowers the energy barrier to remove a 360 • DW, thus rendering the topological pro-tection less effective.For example, in the case of 30 % disorder, all the hysteresis curves collapse to −1 before B ext has reached -25 mT.Before collapsing, the DW width w behaves like w ∼ |B ext | −1/2 as demonstrated in Fig. 5.The inverse square root relation between w and the uniaxial anisotropy constant is a well-known result [37] for stationary 180 • DWs, and our numerical examination verifies that B ext affects the immobile 360 • DWs in a similar way as the uniaxial magnetocrystalline anisotropy affects 180 • DWs.The disorder strength does not seem to affect the w(|B ext |) relation.

IV. CONCLUSIONS
To conclude, our full micromagnetic simulations suggest that Barkhausen noise in Py thin films obeys neither predictions of nucleation nor front propagation models, but is instead a consequence of gradual formation of immobile 360 • DWs via a sequence of abrupt magnetization rotation events.The criticality exhibited by the system appears to be disorder-induced such that the cutoff avalanche size is disorder-dependent.However, the "infinite avalanche" typically observed in Isingtype models for weak disorder is absent here.Instead, we observe smooth hysteresis curves for weak disorder, with progressively larger Barkhausen jumps as disorder is made stronger so that the critical disorder strength is approached from below.Disorder also controls the morphology of the ensuing 360 • DWs such that for weak disorder almost straight DWs spanning the system are formed.Upon increasing the disorder strength, the DWs become increasingly rough, and finally form an irregular DW structure without clear, spanning 360 • DWs at the critical disorder strength.The property of the Py samples not having any magnetocrystalline anisotropy results in the magnetization reversal process happening in a sequence of gradual and localized rotations of the magnetization, highlighting the importance of using full micromagnetic simulations to properly capture the details of the magnetization reversal process [38], and calling for experimental verification of our results, e.g., using magneto-optical imaging [39].
where the first integral is over the volume and the second integral over the surface of the system, n being the unit surface normal vector.These formulations are implemented in mumax3 as given in Ref. [1].
Justification of quasistatic mode.We compare the quasistatic and dynamic runs by simulating the magnetic reversal process in similar systems with both modes.To reduce the computational cost of the dynamical simulations, we begin the runs by driving the system quasistatically with the field step |∆B ext | = 20 µT until m x is below 0.99, and the run is continued using the dynamic mode (time-integration of the LLG equation).Supplemental Fig. 1a shows halves of the hysteresis loops for simulations run in the quasistatic and dynamic modes.In the dynamic run, we use the RK45 algorithm [1,2] with an adaptive time step and apply the sinusoidal external field with frequency f = 100 kHz.The Gilbert damping constant has the value of α = 0.02.In the quasistatic run, the field step |∆B ext | = 4 µT is used.We see that the curves are qualitatively similar: the jumps are roughly the same size and they occur at roughly the same strength of the external field.This is a strong indication that the microscopic magnetization behaves locally the same way, that is the same magnetic structures and patterns are formed in the quasistatic and dynamical runs.We have confirmed that this is the case for the particular run by comparing the resulting magnetization configurations.
The details in the hysteresis curves are slightly different, however, as demonstrated in Supplemental Fig. 1b.The curves are shown for multiple frequencies in the dynamical runs.Focusing first on the 100 kHz curve, the most clear difference is that the Barkhausen jumps are abrupt in the quasistatic simulations, and they appear delayed in the dynamical curves due to the slowness of the response.After each jump, the dynamical curve experiences some vibration, which is then damped.The vibration is attributed to the damping precession of the local magnetizations, and its frequency is related to the material parameters in the LLG equation.
In Supplemental Fig. 1b, we further see that decreasing the frequency down from 10 MHz to 100 kHz leads to the dynamical curves approaching the quasistatic one.Based on this solid trend, we conclude that the intuitive hypothesis of the quasistatic mode corresponding to the limit f → 0 is valid.
To further support the use of the quasistatic mode, let us discuss the computational requirements of the two modes.The system size of 1024 nm × 1024 nm × 16 nm and frequency of f = 100 kHz, as used in the comparison above, is at the computational limit of the dynamical mode with dynamics governed by the RK45 algorithm: simulations of the dynamics of larger systems or smaller frequencies would lead to an unreasonably large requirement of computational time.However, for reliable Barkhausen jump statistics, we want to simulate larger system sizes to account for different emerging magnetic structures and their behaviour under changing applied field.On the other hand, the relevant frequency range in applications and experimental research is well below 1 kHz.Therefore, at this level of theory, we conclude that there are insurmountable challenges in modeling the Barkhausen noise statistics with dynamical simulations.We will run our main simulations with quasistatic mode, that can be used to mimic the dynamic runs with small frequencies.Also, larger system sizes can be simulated  with the quasistatic mode, although a field step has to be selected that is small enough for reasonable accu-racy and large enough for reasonable computational time.
Effect of the field step.In Supplemental Fig. 2, we show jump size distributions for the systems of 30 % disorder, computed quasistatically with two different field step magnitudes.The distributions include 40 individual runs for both the field steps.The distributions are plotted both against |∆m x /∆B ext | and against |∆m x |, where the latter corresponds to the figures in the main text.We see that the position of the bump is changed in the distribution for |∆m x |, but not for the distribution for the quantity |∆m x /∆B ext |.This supports the conclusion presented in the main text that the bump, as well as the part on its left hand side, corresponds to the smooth rotation of the magnetization.The part on the right hand side of the bump is preserved when the field step is changed, supporting the presence of the abrupt Barkhausen jumps in that size region.There are no essential differences in the shape of the Barkhausen jump part of the distribution.Therefore, we conclude that the field step of 4 µT that is used for the main results is sufficiently low for accurate results.
Reversal processes with 30 % disorder In Supplemental Fig. 3, different types of magnetization reversal processes in 30 % disorder systems are illustrated along with the respective hysteresis curves.The selected runs represent the four subgroups that can be distinguished in the group of the hysteresis curves.We attribute the different subgroups to different densities of the forming 360 • DWs.The large vertical jumps in the hysteresis curves correspond to vanishing of some of the DWs.

FIG. 1 .
FIG. 1.(a) Example of the spatial variation of Ms, sampled from the normal distribution with mean 800 kA/m and standard deviation of 30 % of the mean value.(b) Example evolution of mx with quasistatically decreasing Bext.

FIG. 2 .
FIG. 2. Hysteresis curves for mx(Bext) of all simulations.The external field Bext is decreased from 100 mT to −100 mT.The disorder strength is indicated in the upper left corner of each panel.

FIG. 4 .
FIG. 4. Examples of the magnetization reversal mechanism with (a) 5 %, (b) 20 %, (c) 30 %, and (d) 35 % disorder strengths.The arrows denote the local magnetization direction in every 64th grid point.The background color indicates the direction of the local magnetization m as given by the color bar; θ denotes the angle between the x-axis and m in radians.The Bext values corresponding to each snapshot are indicated in the top left corners.
FIG.5.360 • DW width w as a function of |Bext| in selected runs with varying disorder strength where two 360 • DWs are formed.w evolves approximately as |Bext| −1/2 (dashed line).The inset shows an example configuration where w is evaluated at Bext = −5.6 mT for 5 % disorder.w is taken to be the distance between the two points at which the magnetization profile mx(x), averaged over y, crosses 0.
arXiv:2209.10847v2 [cond-mat.mtrl-sci]6 Mar 2023 SUPPLEMENTAL FIG. 1.Comparison of the hysteresis curves between dynamic and quasistatic runs.The system size is 1024 nm × 1024 nm × 16 nm, and the disorder strength is 40 %.(a) Halves of the hysteresis loops with quasistatic and dynamic modes.(b) A more detailed look at a single magnetization jump along the curve, showing the quasistatic result and the corresponding results from dynamic simulations with a wide range of field frequencies f .Notice how the dynamic curves gradually approach the quasistatic one as f gets smaller.

2 .
Comparison of the Barkhausen signal distributions with different field steps for disorder strength of 30 %.(a) Plot against s = |∆mx/∆Bext| in units of mT −1 .(b) Plot against s = |∆mx|.The magnitudes of the field steps are given in the legend.The dashed line indicates the fit power law with the exponent τ = 1.75 in accordance to the main text.

SUPPLEMENTAL FIG. 3 .
Magnetic configurations with 360 • domain walls from selected runs with 30% disorder strength.In the runs, (a) 2, (b) 4, (c) 6 distinct domain walls are formed.In (d), more walls are beginning to form, but some of them collapse during the run; at −13.71 mT, 6 domain walls can be distinguished as shown in the inset.The inset coloring corresponds to that of the main text.The respective hysteresis curves are drawn as thick, black lines.SUPPLEMENTAL FIG. 4. Magnetic configurations in a validation run with open boundary conditions.The disorder strength is 30 %.The system size is 2 µm × 2µm with thickness of 16 nm.The field step is 0.1 µT.Except for the formation of closure domain patterns near sample edges, similar dynamics of the magnetization to those in Fig. S3 are observed, validating that the formation of 360 • DWs is not due to introduction of the periodic boundary conditions in the system.