Micromagnetic simulations of first-order reversal curves in nanowire arrays using MuMax3

We perform simulations of magnetic reversal in a 3 × 3 array of nanowires using MuMax3 micromagnetic simulation program. We record a series of first-order reversal curves (FORCs) that form distinct branches of ascending minor curves depending on the initial magnetization state. We calculate the FORC distribution, which shows 9 positive primary peaks, representing single reversals of the 9 simulated nanowires. The primary peaks form an interaction field distribution (IFD), a common feature in experimental FORC distributions due to demagnetizing interactions. The FORC distribution also contains positive and negative secondary peaks due to differing magnetization during reversal. We demonstrate the use of MuMax3 simulations to relate FORC distribution features to visualized magnetic configurations.


I. INTRODUCTION
Magnetic nanowire arrays are promising nanostructures for applications in high-density magnetic recording, 1,2 microwave devices, 3,4 spintronic devices, 5,6 and sensors. 7,8 The magnetic reversal fields of these arrays can be mapped by the first-order reversal curve (FORC) distribution. 9 However, large interaction fields can reshape the FORC distribution of arrays and influence the overall switching behavior. 10,11 Better understanding of nanowire array switching is needed to mitigate or utilize interaction effects in future devices.
Previous simulations of nanowire arrays have modeled grids of point-dipoles, 12 grids of point-elements, [13][14][15] or collections of hysteretic elements, 11 which exhibit FORC distribution features similar to those observed experimentally. These models allow simulation of arrays composed of many nanowires but at the cost of approximating or neglecting intra-wire interactions and inhomogeneous states. Our goal has been to simulate the reversal processes in nanowire arrays and to map changes in magnetization to features in the FORC diagram. We have used micromagnetic simulations of 3-D nanowires to investigate how the internal magnetization reversals, including mixed or inhomogeneous states, of the nanowires in the system appear in FORC distributions. Since full micromagnetic simulations are more computationally intensive than simulations of simplified models, we are limited to small arrays. In this study, we have simulated 3 × 3 and 4 × 4 nanowire arrays. Here, we report results for the 3 × 3 array, noting that the 4 × 4 array exhibits similar results.

II. MICROMAGNETIC SIMULATION
Simulation has been performed with the open-source program MuMax3. MuMax3 simulates magnetization dynamics in nanoand microscale ferromagnetic structures. 16 The magnetization at all points in the model is solved simultaneously and over time from the Landau-Lifshitz-Gilbert equation. 17 The program computes finitedifference solutions in CUDA-enabled GPUs. 16 Prior to simulation, we have tested and verified MuMax3 by matching the published results for the Standard Problems. 18 The 3 × 3 array of nanowires is composed of identical cylinders with diameter (D = 80 nm) and length (L = 1 μm). The chosen length is in the regime for which modeling nanowires as pointdipoles may be limited. 19 The wires are arranged in a 2-D square ARTICLE scitation.org/journal/adv grid with equal spacing between nanowire centers (r = 200 nm). Every nanowire has the same exchange constant (Aex = 1.3×10 −11 J m −1 ), saturation magnetization (Ms = 800 kA m −1 ), and damping (α = 0.008). Therefore, we have not introduced differences in reversal behavior of the nanowires themselves, which would otherwise alter the FORC distribution. 15,20 The system is discretized into 4 nm × 4 nm × 4 nm cells, width smaller than the exchange length (lex = √ 2Aex/μ 0 M 2 s ≤ 5.69 nm). Thermal effects are not accounted for.
After the model is constructed, and before solving the Landau-Lifshitz-Gilbert equation, the demagnetization kernel is calculated for each finite-difference cell in the simulation. MuMax3 does this by assuming a uniform magnetization for each cell, directly calculating the demagnetization field, solving for the demagnetizing field kernel, then convolving with the magnetization of a particular state. 16 When determining the kernel, MuMax3 allows for periodic boundary conditions that extend the system. This significantly reduces run-time due to the lack of zero-padding required for the Fourier transforms. 16 This also prevents interactions being weaker for nanowires along the edges or corners of the simulation window for better representation of an infinite array. 15 Our simulations repeat the system four times on each side of the simulation window.

III. FORC SIMULATION
The simulated external magnetic field is varied similarly to the applied field in experimental FORC measurements. First, a field sufficient to saturate the system is applied (Hs = 119 kA m −1 ). The descending major curve is simulated while the 3-D magnetization and initial external field (HR) are recorded. Unlike experiment, where saturation and the descending sweep are performed between each FORC, we initialize each FORC using these saved configurations to reduce total run time by 50%.
A particular FORC, or ascending minor curve, is then measured by recording the normalized magnetic moment (m) as the external field (H) is increased to positive saturation at Hs. The net moment and the external magnetic field are recorded as two 3-D vectors, although only the components parallel to the nanowire axis are analyzed. This process is repeated for 60 values of H R between +Hs and −Hs. This provides a data set describing the magnetic moment as a function of the reversal field and the external field [m(HR, H)], known as the FORC family ( Fig. 1). Since each FORC is an independent simulation from known initial states, multiple FORCs can be simulated in parallel to further reduce total run time.
The low H portion of some curves forms a branch, a visible vertical gap indicating a greater m than that of the curve labeled "Low" at the same H. The branches demonstrate hysteresis, a dependence of m on the initial state at HR as well as on H commonly seen in FORC measurement. Above a particular value of H and below a value of HR, the FORCs merge with Low. We refer to this threshold (H, HR) as the "merge point." There are 9 distinct hysteretic branches, labeled "A" to "I," which we have verified on an (H, HR) plot after subtracting the Low curve for improved contrast. Each branch is composed of several overlapping FORCs (Fig. 1 inset). The curves have similar slope at H values between merge points. The sloped component of m is dependent on H but independent of HR (discussed below).
At fields above the merge point, some FORCs exhibit temporary separation from the Low curve. The states along different FORCs but at the same H have nearly equal m except near the merge points.

IV. FORC ANALYSIS
The analysis method of the simulated FORC data is identical to that of experimental FORC data. The Preisach model describes magnetic reversal as the collective behavior of elementary hysteretic elements called "hysterons." 21 The FORC density, representing change in m due to a single hysteron reversal, is calculated as the mixed second derivative of FORC data, 9 ρ(H, HR) = − 1 2 This density calculated over the (H, HR) half-plane (H > HR) is the FORC distribution. 9 The FORC distribution is often displayed as a contour or density plot in the basis of the interaction field [Hu = (H + HR)/2] and the coercive field [Hc = (H − HR)/2]. These new fields represent the center (interaction) and the hysteresis width (coercivity) of each hysteron. 20 The resulting FORC distribution [ρ(Hc, Hu) calculated using Eq. (1)] is shown in Fig. 2. The hysteretic branch features of Fig. 1 can be described by their merge points, labeled by branch in Fig. 2. Each merge point has a corresponding positive primary peak in the calculated FORC distribution. Positive and negative secondary peaks are seen at similar H but lower HR.
Dobrotȃ and Stancu have observed comparable FORC distributions of a simulated array of 1600 nanowires, simulated as point-elements. Individual nanowires reverse at different field values between FORCs, creating multiple positive and negative contributions across the (H, HR) plane with fine field steps. 15  Stancu verify the net contribution of a single nanowire is a unit of m, identified by a single positive primary contribution at highest HR. 15 We observe similar notable contributions to m (merge points) from 9 primary peaks, possibly one for each nanowire. The observed hysterons form a ridge parallel to the Hu axis. This feature is comparable to an intrinsic field distribution (IFD) in experimental nanowire array FORC distributions. The IFD has been attributed to geometry-dependent demagnetizing interactions between nanowires in the array. 19,22 In the moving Preisach model, hysterons are shifted from their intrinsic positions on the Hc axis by the mean interaction field (H int ). 10,11 This model has been used to interpret experimental results. 11,19 The intrinsic field that causes hysterons to reverse is Hc = H + mHint. 10 Linear fit of the merge point fields versus moments provide intrinsic hysteron coercivity (Hc = 1.18 ± 0.07 kAm −1 ) and maximum mean interaction field (H int = −111.3 ± 0.1 kAm −1 ), or effective array axial demagnetizing field (N axial = −H int /Ms = 0.1391 ± 0.0001) with small root mean square error (RMSE = 0.4 kAm −1 ). The small error in the linear trend suggests primary reversals occur with similar coercivity and an interaction field that is uniform between elements.
We do not observe peaks forming a clear coercive field distribution (CFD), which often appears in experimental FORC distributions. 22 From the simulation study by Dobtrotȃ and Stancu, this can be attributed to our number of nanowires being small compared to the number of field steps. 15 The individual contributions produce finer features in the FORC distribution instead of statistical distributions. 15 As mentioned, the constant-slope component of m in Fig. 1 is independent of HR. As expected from Eq. (1), the sloped component does not have a corresponding positive density in the FORC distribution, which only shows hysteretic or irreversible change in magnetic moment. 9 Instead, the sloped component represents a large reversible change in magnetic moment.

V. VISUALIZATION
3-D visualization allows us to deduce the mechanisms responsible for the three observed FORC features: hysteretic branches (primary peaks at merge points), variations between FORCs (secondary peaks at equal H and lower HR), and non-hysteretic slope (all H but not visible in the FORC distribution). Fig. 3 shows the 3-D visualization comparing the magnetization of FORCs starting in branches D and Low near the merge point D. From Fig. 2, FORCs from initialization to merge point D should differ from the Low curve by the contribution of peak D.
In images L(i) to L(iv) of Fig. 3, nanowire β2 evolves from a mixed state to saturated up. Nanowire γ2 appears nearly saturated in L(i) but relaxes to a mixed state in L(iv), so we consider this  Table I. The images have varying number of nanowires saturated up (blue) and down (red), demonstrating hysteresis. The labeling convention of nanowires is overlaid in D(i); for example, the center nanowire is referred to as "β2." ARTICLE scitation.org/journal/adv nanowire to be mixed in L(i) to L(iv). Nanowires α3, β1, and γ1 are saturated in the initial state. The number of nanowires fully saturated up along branch Low changes from 3 to 4. From D(ii) to D(iii), nanowire β2 evolves from saturated down to a mixed state. The number of nanowires saturated up in branch D remains constant at 4. Therefore, merge point D and its primary peak represent a net increase in the number of saturated nanowires by one, showing agreement with the study by Dobrotȃ and Stancu despite simulating fewer nanowires. 15 Images of L(ii) to L(iv) and D(ii) to D(iv) show states at H near that of the merge point D. The magnetizations in wires β2 and γ2 are visibly different between images L(iv) and D(iv). Despite this, the net moment values of these states are nearly equal as seen in Table I and the Fig. 1 inset. Dobrotȃ and Stancu observe differences between FORCs at the same H due to the coercivity of the particular reversed nanowire, resulting in secondary features in the FORC distribution. 15 Although we simulated identical nanowires, the nanowire reversal field may be influenced by different inhomogeneous states before reversal.
Domain growth in mixed-state nanowires is gradual, seen between each image of both branches in Fig. 3. Domain growth even occurs at fields between the merge points or observed primary peaks. Thus, the resulting m accounts for the sloped component of the FORCs which, as discussed earlier, does not result in FORC peaks. The imaged states show primary peak D only describes the late stages of reversal (saturation) in a single nanowire, β2 in the case of the Low curve. This is a newly described detail enabled by micromagnetic simulation.
3-D magnetization of states along the Low curve is shown in Fig. 4. States have been chosen such that the number of nanowires saturated up increases by 1 in each consecutive image. The imaged measurements are shown as points in Fig. 2. Most of these images are states at H between those of the 9 primary peaks, supporting the hypothesis that secondary peaks and a primary peak  Table II and labeled in Fig. 1.  15 Discrepancies in H may be accounted for by contributions of secondary peaks, related to variations in switching field between FORCs.

VI. CONCLUSION
Micromagnetic simulation allows combining 3-D visualization with standard FORC analysis. The FORC distribution shows 9 primary peaks, which are also seen as the family of FORCs forming separate branches. From visualization, these peaks are found to represent saturation of a single additional nanowire. Gradual domain growth in mixed-state nanowires is not represented in the FORC distribution. The peaks occupy different Hu positions, forming an IFD, attributed to demagnetizing interactions between elements in the array. 12 Our simulation suggests that identical nanowires in a region of an array would reverse one-by-one at different external fields. This verifies common explanations for the IFD in FORC distributions of nanowire arrays in experiment 12,22 and simulation. 12,15