gFORC: A graphics processing unit accelerated first-order reversal-curve calculator

First-order reversal-curves have proven to be an indispensable characterization tool for physics as well as for geology. However, the conventional evaluation algorithm requires a lot of computational e ﬀ ort for a comparable easy task to overcome measurement noise. In this work, we present a new evaluation approach, which exploits the diversity of Fourier space to not only speed up the calculation by a factor of 1000 but also move away from the conventional smoothing factor toward real ﬁ eld resolution. By comparing the baseline resolution of the new and the old algorithm, we are able to deduce an analytical equation that converts the smoothing factor into ﬁ eld resolution, making the old and new algorithm comparable. We ﬁ nd excellent agreement not only for various systems of increasing complexity but also over a large range of smoothing factors. The achieved speedup enables us to calculate a large number of ﬁ rst-order reversal-curve diagrams with increasing smoothing factor allowing for an autocorrelation approach to ﬁ nd a hard criterion for the optimum smoothing factor. This previously computational prohibitive evaluation of ﬁ rst-order reversal-curves solves the problem of over-and undersmoothing by increasing general readability and preventing information destruction.


I. INTRODUCTION
Hysteresis measurements are an essential tool to investigate the magnetic behavior of a large variety of systems and characterize them by basic parameters such as coercivity field, remanent, and saturation magnetization.Although a single major loop measurement already contains valuable details about the sample, it is possible to maximize the information obtained from magnetization measurements by systematically collecting data within the major loop in different magnetic paths or minor loops.The information obtained from comparing successive minor loops is not accessible with a mere major loop measurement and, following Preisach's model of hysteresis, 1 it reveals the system's coercivities and interactions at once.Two virtually congruent major loops can have a fundamentally different composition of minor loops. 2,3First-order reversal-curves (FORCs) 4 succeed at disentangling different magnetic contributions and allow us to distinguish between various magnetic phases within one single sample.This capability makes FORC a versatile technique applicable to a broad spectrum of samples.
In principle, collecting FORC data is possible with any device capable of measuring hysteresis loops. 28A measurement starts in positive saturation. 23The field is then decreased to the first reversal field H r , usually negative saturation.The measured data consist of the M-H minor loop starting at H r and ending at positive saturation.This procedure is then repeated for increasing reversal fields H r until the minor loop starts in positive saturation.Thus, the measured data can be displayed as a two-dimensional magnetization landscape, which not only depends on the applied field H but also on the reversal field H r .
The so-called FORC density is calculated in two steps.First, the FORC diagram ρ is calculated from the magnetization landscape M(H, H r ) by taking the mixed second derivative along H and H r ρ(H, H r ) ; À 1 2 The second part is a coordinate transformation of H and H r into a more intuitive coercivity and an interaction field, H c and H u , respectively, This transformation equals a rotation and expansion of the coordinate system. 23Frequently done but not strictly needed, the FORC density can be illustrated in H À H r -coordinates as well.
Numerically evaluating Eq. ( 1) can be quite challenging since derivatives are already susceptible to small amounts of noise, even worse for a mixed second derivative.Nevertheless, a finite differences approach has recently been published by Abugri et al., 29 where smoothing was not a valid option.A similar strategy is part of the FORCit software by Acton et al., 30 which calculates the derivative after a real space Gaussian filter.
The most prominent method for the calculation of FORC densities was proposed by Pike et al. 23 Instead of directly calculating the mixed second derivative as a difference quotient, they proposed to interpolate the magnetization data on a rectangular grid in order to locally fit a polynomial p of the form to the magnetization.The amount of data points included in the fit depends on the so-called "smoothing factor" SF.When substituting p into Eq.( 1) and calculating the mixed second derivative, it can be directly seen that only Àa 6 =2 remains and thus equals the FORC density at that given position.][33][34] However, this algorithm has some inherent disadvantages.
First of all, it calculates a two-dimensional fit at every point of the interpolation grid, which is in the order of 250 000 fits.That typically means a calculation time of up to 10 min on modern computers for a reasonable resolution.Moreover, since the smoothing factor defines the number of surrounding data points included in the fit on a rectangular grid (N Data ¼ (2SF þ 1) 2 ), it is not possible to choose noninteger smoothing factors.Thus, the resolution of a FORC diagram not only depends on the smoothing factor but also on the grid that is chosen for the data interpolation.The smoothing factor is virtually meaningless without the grid spacing.
There are multiple extensions and variations of this conventional fitting algorithm.The first one, "extended FORC," proposed by Pike, is able to capture fully reversible information in a "reversible ridge" at H c ¼ 0 by extending each minor loop. 31FORCinel" by Harrison and Feinberg 32 calculates FORC diagrams based on a locally weighted regression smoothing, often referred to as LOESS 35 algorithm.The distance weighted and, therefore, not squarelike smoothing kernel allows for a noninteger smoothing factor.
"VARIFORC" by Egli 33 is able to dynamically adapt the smoothing in order to meet the requirement of high resolution or noise reduction for different regions of the FORC diagram.
Cimpoesu et al. presented a versatile fitting tool for first-order reversal-curves ("doFORC" 34 ), which is able to use a great variety of different smoothing kernels for the local regression.However, the calculation of a large amount of fits persists throughout all of these different adaptations.
Another challenge is the choice of the optimum smoothing factor.FORC diagrams can be hard to read and are not always straightforward to interpret. 2,27For a correct interpretation, the choice of the smoothing factor is crucial.It is relatively easy to determine whether a FORC diagram is undersmoothed since it mainly consists of noise; however, it is much harder to identify a slightly oversmoothed FORC diagram.There have been approaches to numerically determine the optimum smoothing factor. 32,36owever, with several minutes computing time per smoothing factor, this can be a long and exhausting task to perform.New approaches that allow for a shorter computing time are, therefore, not only important for convenience but also necessary for a more reliable interpretation of FORC diagrams.
In this work, we present an alternative approach for the calculation of FORC densities based on Fourier filtering and derivatives.To verify the approach, we compare FORC diagrams of four different samples at identical field resolutions and over a large range of smoothing factors.Achieving a severe speed up of the FORC density calculation allows us to reiterate on previously computational prohibitive approaches to find the optimum smoothing factor of a FORC diagram and apply it to the presented measurements.

II. EXPERIMENTAL SECTION A. Sample preparation
In this study, two types of samples are investigated, thin film stripe structures and a macroscopic permanent magnets.As a guide to the eye, the stripe measurements are labeled with green, light blue, and orange, and the permanent magnet measurements in violet.
The 1 Â 1 Â 2 mm 3 bulk permanent magnet is cut from a commercially available 1 Â 1 Â 1 cm 3 neodymium magnet (NdFeB) via wire electrical discharge machining.Its easy axis is aligned along the long edge of the cuboid.Thin film stripe structures were fabricated via photolithography and direct laser writing using the "LOR 3A" and "AZ ECI 3027" resist by "MicroChem" and "MicroChemicals," respectively.For exposure, the UV laser system "KLOÉ Dilase 250" was used.After development, 50 nm permalloy (Fe 20 Ni 80 , Py) and 2 nm aluminum (Al) were deposited via ion beam sputtering at a pressure lower than 10 À7 mbar.Further details can be found elsewhere. 27

B. Measurement methods
Permalloy sample measurements were conducted using "Durham Magneto Optics' NanoMOKE3," which makes use of the magneto-optical Kerr effect (MOKE).For more information on the MOKE-FORC setup, the reader is referred elsewhere. 27,37ulk measurements were carried out using "Quantum Design's MPMS 3" vibrating sample magnetometer SQUID (VSM-SQUID).A fast ramping 7 T magnet is able to saturate the NdFeB magnet while still reducing the measurement time by approximately a factor of five compared to a conventional SQUID.All FORC measurements presented consist of approximately 100 to 150 minor loops.The reversal field spacing was distributed nonuniformly in order to measure more minor loops in interesting field regions while keeping the measurement time at a minimum.An in-house written acquisition and processing software is able to collect the MOKE as well as the SQUID data and subsequently process it.

A. A new approach
The first step of the calculation is the interpolation of the measurement data onto a rectangular grid.This step automatically ensures that the algorithm is robust enough to also handle nonuniformly spaced reversal field distributions.
Since the Fourier transformation expects a periodic function, simply Fourier transforming each minor loop will cause artifacts after the inverse transformation.Therefore, each minor loop is extended by a flipped copy of itself, which results in a smooth periodic signal.Subsequently, the fast Fourier transformation (FFT) M H (H, H r ) of M(H, H r ) along H is calculated.H denotes the H-spectrum.
In Fourier space, any arbitrary filter can be applied; usually, a Gaussian works well.However, any frequency filter can be applied, e.g., known instrument noise frequencies can be eliminated selectively.The derivative is calculated by multiplying the filtered M * H (H, H r ) with iH in Fourier space where i denotes the imaginary unit.The inverse transformation (iFFT) of iHM * H (H, H r ) yields a smoothed @ H M(H, H r ).ρ(H, H r ) is obtained by repeating the same process along the H r -axis.The total calculation process is as follows: The advantage of this method is that it does not calculate a six parameter fit at every data point anymore.The algorithm only consists of Fourier transformations, a multiplication in Fourier space, and an inverse Fourier transformation.These operations are inherently fast and parallelized on graphics processing units (GPUs).The performance of the algorithm was tested on an "Intel Xeon E3-1505M v6" with an "Nvidia Quadro M1200" and an "Intel Core i7-8850H" with an "Nvidia Quadro P2000".The calculation time of a full FORC density, which used to be in the order of minutes, is reduced to approximately 0.2 s.

B. Comparing smoothing factors
To compare the results obtained from the proposed and the conventional algorithm, a common resolution is necessary.Since the proposed FFT algorithm expresses resolutions in absolute field units whereas the conventional fitting algorithm expresses resolutions in a combination of smoothing factors and field step size, a transformation is necessary.As it is nontrivial to compare the resolution of the squarelike smoothing kernel of the conventional fitting algorithm to the resolution given by the Gaussian in Fourier space, we instead compare the baseline resolution of both algorithms.For the conventional algorithm, the baseline resolution is directly given by the smoothing factor SF and the field step ΔH of the interpolation grid The baseline resolution of two neighboring Gaussians is given by their width. 38If the inverse width in Fourier space is 1=σ, it can be shown that the baseline resolution in real space is given by From Eqs. ( 4) and ( 5) and the condition for those two resolutions to be equal R proposed ¼ !R conventional , a smoothing factor equivalent SF eq of the new algorithm can be calculated from the Gaussian width σ and the field step given by the interpolation grid ΔH

C. Verification
In the following, we verify the proposed algorithm against standard data sets.Using Eq. ( 6), FORC diagrams of the conventional and proposed algorithm become comparable via the smoothing factor and the smoothing factor equivalent. Figure 1 displays three different verification measurements.A more extended verification against simulated data sets can be found in the supplementary material.
The first column displays the three different sample geometries.The first two measurements [(green, light blue, Figs.1(a)-1(f)] are conducted on the same thin film permalloy stripe sample (150 μm Â 10 μm Â 50 nm) for two different field geometries.The third measurement [ purple, Figs.1(g)-1(i)] is conducted on bulk NdFeB (1 Â 1 Â 2 mm 3 ).The second column displays the major and minor loops of the corresponding sample.
The two FORC densities of the first sample [green, Figs.1(b) and 1(c)] are virtually identical.Both have a peak at the coercivity expected from the major loop and at H u ¼ 0. However, the two FORC diagrams do not perfectly match.First, the noise level of the FFT algorithm is slightly superior but secondly and more important is that the peak shape is different.While the proposed FFT algorithms result in a round peak shape, the FORC peak of the fitting algorithm adopts the shape of the smoothing kernel.Obviously, a square shaped smoothing kernel will cause sharp peaks in the FORC density to have a squarelike shape.There are several approaches to avoid such artifacts.This can, for example, be done by taking a different smoothing kernel 34 or by applying a distance weighted fitting algorithm. 32The FFT algorithm automatically avoids such artifacts by using a Gaussian filter in Fourier space.
To verify the correct processing of fully reversible information, the second measurement (light blue) was conducted on the same sample as the first on (green) but with a different field geometry.Whereas the field was applied along the easy axis in the previous case, the field is applied along the hard axis of the stripe for the second measurement as seen in the sample column in Fig. 1.
Such a hard axis system will only experience continuous rotation of the magnetic moments without any irreversible switching as seen in the major and minor loops [Fig.1(d)].The extended FORC technique 31 ensures that also fully reversible information at This extended FORC technique is also implemented in the FFT algorithm.It ensures that information of fully reversible processes is displayed as part of a so-called reversible ridge 31 at H c ¼ 0. It also ensures that the integral over the entire FORC density equals the saturation magnetization, a necessary condition that follows from the way the FORC density is calculated [cf.To prove that the FFT algorithm can reliably recover the correct FORC diagram even for complex sample compositions, a third measurement ( purple) was conducted on a macroscopic NdFeB sample.This sample consists of a hard magnetic component at H c % 1:5 T, a soft magnetic component at H c % 0 T, and a positive-negative peak pair at H c % 0:7 T and H u % À0:6 T, which originates as a consequence of the interaction between the hard and the soft magnetic component in the sample. 27Further interpretation of the FORC diagram is beyond the scope of this work.The important point is that the FFT algorithm is able to recover all peaks in the FORC diagram reliably.
To further validate the FFT technique not only for different samples but also for various smoothing factors, an additional thin film permalloy stripe sample has been measured.The sample consists of 150 μm long permalloy stripes of alternating width (10 μm and 30 μm) with a thickness of 50 nm and a spacing of 10 μm.
To measure this sample with MOKE, the laser spot is slightly defocused in order to capture not only the magnetization of one single stripe but to include the full structure into the magnetization measurement.
A sketch of the sample and the major and minor loop measurement can be seen in Figs.2(a  The FORC diagrams are characteristic for a two component system with two different coercivities.Two peaks on the H c -axis correspond to the respective coercivities of the system and a positive-negative peak pair at negative H u arises as a consequence of the interaction.A thorough analysis of such FORC diagrams is once again beyond the scope of this work and can be found elsewhere. 27he FORC diagrams of the FFT [Fig.2(c)] and the fitting algorithm [Fig.2(d)] are almost identical over a large range of smoothing factors.For the two smaller smoothing factors (SF ¼ 1, 5) the FORC densities are basically identical.For larger smoothing factors (SF ¼ 10, 15), they are still similar; however, the peaks of the fitting algorithm start to adopt the square like shape if the fitting kernel is much larger than the peak itself.
Each FORC diagram has an optimum smoothing factor at which noise is sufficiently suppressed while no information is destroyed, in this case SF % 5. Since the FFT algorithms allow for an arbitrary width of the Gaussian Fourier filter, the smoothing factor can have any real value SF eq !À0:5, where À0:5 equals no smoothing [σ ¼ 0 Oe in Eq. ( 6)].
Judging whether a FORC density is well-smoothed and not over-or undersmoothed can be challenging.There have been approaches to find hard criteria for the optimum smoothing parameters. 32,36However, they rely on the calculation of a large number of FORC densities over an extended range of smoothing

Journal of Applied Physics
factors, which is very time consuming if done with the conventional fitting algorithm.

D. Exploring new possibilities
Reducing the time needed to calculate a FORC diagram by a factor of approximately 1000 allows for the exploration of science, which was computational prohibitive before.We use the autocorrelation approach described in Secs.4-6 in Ref. 36, which is based on the autocorrelation of the residual.The residual for a given smoothing factor is calculated as the difference of the unsmoothed magnetization landscape and the magnetization landscape of the corresponding smoothing factor that is obtained by integrating the FORC density along H and H r , The residual is the difference between the actually measured magnetization landscape and the landscape that corresponds to the smoothed FORC density.
If the residual is perfectly flat the corresponding FORC equals the unsmoothed FORC and contains a lot of noise.When increasing the smoothing factor, the residual will become noisier but there will not be any significant features in it.If the smoothing is further increased, the residual will start to form features.This is the point at which the smoothing starts to destroy features of the FORC diagram and separated peaks start to merge.
By calculating the sum over the autocorrelation matrix of the residual for a large number of different smoothing factors, these three regions can be distinguished.The blue graph in Fig. 3 displays the normalized sum of the autocorrelation matrix according to Eq. ( 12) of Ref. 36 (same sample as in Fig. 2).Each of the 500 data points corresponds to the calculation of a FORC diagram and the calculation of the residual autocorrelation over the ten nearest neighbors.In the optimum case, this plot consists of three different regions.The slope of the correlation in the first region (left gray box) corresponds to the reduction of noise in the FORC diagram.The flat part and second region (green box) corresponds to a wellsmoothed FORC diagram and is the main region of interest.In this region, the noise reduction by smoothing is sufficient and the signal is not reduced yet.In the third region (right gray box), the correlation starts to increase dramatically, which indicates the merging of peaks and the loss of information.
The optimum smoothing factor can be identified by looking at the red graph, which displays the normalized derivative of the residual correlation, and taking the minimum of it.This ensures that the identified smoothing factor is neither too close to the point of noise nor to the point of signal reduction.
Although it takes a considerable amount of computational time, this is a hard criterion for the optimum smoothing factor of a FORC diagram and also outputs a range of acceptable smoothing factors as well.This has been facilitated by the speedup of the FORC calculation without which it would have been impossible to calculate FORCs for so many smoothing factors in the first place.

IV. SUMMARY
A new approach for the calculation of FORC densities based on Fourier filtering and multiplications in Fourier space has been presented.This approach reduces the calculation time per FORC density by a factor of 1000 compared to the conventional six parameter fitting algorithm. 23By comparing the baseline resolutions of the fitting algorithm and the new Fourier approach, the well-established concept of smoothing factors is comparable to real field smoothing resolutions.This is a thorough step toward a universal comparability of FORC diagram resolutions.
The new approach has been verified by comparing three different sample geometries for identical smoothing resolutions.We find excellent agreement between the fitting algorithm and the newly developed Fourier algorithm for all presented samples.

Journal of Applied Physics
To further verify the approach, we compared FORC densities of the same samples over a large range of smoothing factors.All FORC densities, whether under-or oversmoothed, show an exceptionally good agreement with the conventional algorithm.
The calculation speedup allows us to evaluate FORC densities over a large range of noninteger smoothing factors with small spacing.This allows for the visualization of the evolution of a FORC densities as a function of smoothing factor that highlights over-and undersmoothing.Furthermore, the reduction of the calculation time allows us to reevaluate a previously computationally prohibitive technique to calculate the optimum smoothing factor.The speedup allows us to calculate the autocorrelation function of the FORC residual as a function of smoothing factor, which gives a hard criterion for the smoothing factor needed for the optimum calculation of the corresponding FORC diagram.

SUPPLEMENTARY MATERIAL
See the supplementary material for a movie of the FORC density evolution as a function of smoothing factor.The supplementary material also includes two more data sets of the autocorrelation of the data sets shown in Figs.1(a)-1(c) and 1(g)-1(i) and an example of the FFT algorithm evaluating data generated by the Preisach model with and without noise.A MATLAB implementation of the algorithm for the evaluation of FORC densities can be found at https://gitlab.gwdg.de/MoKeteam/gFORC/.

FIG. 1 .
FIG. 1.Comparison between the conventional fitting and the proposed FFT algorithm: three different sample geometries are shown (green, light blue, purple).The first column displays sample sketches and the applied field axis.The second column displays the major and a few selected minor loops.The color of the FORC density has been mapped onto the minor loops, which enables direct assignment of FORC density peaks to their origin in the minor loops.The third and fourth columns display the FORCs calculated with the fitting and the FFT algorithm, respectively.The field resolutions used for the FFT smoothing are 1.06 Oe (green), 10.8 Oe (light blue), and 0.18 T (purple).Positive signal in the FORC density is represented by red color, and negative contributions are illustrated in blue.The signal intensity is encoded in color saturation.

Figure 2 (
c) shows the FORC diagrams of the FFT algorithm and Fig. 2(d) displays the results of the fitting algorithm.

FIG. 2 .
FIG. 2. Comparison of FORC densities with different smoothing factors: (a) Sample sketch.(b) Major loop measurement with a few selected minor loops.(c) Series of different smoothing factors (FFT).Real field resolutions are SF eq ¼ 1 ≙ 0:29 Oe, SF eq ¼ 5 ≙ 1:07 Oe, SF eq ¼ 10 ≙ 2:05 Oe, and SF eq ¼ 15 ≙ 3:03 Oe.Positive signal in the FORC density is represented by red color, and negative contributions are illustrated in blue.The signal intensity is encoded in color saturation.(d) Series of different smoothing factors (fit).See supplementary material for a video of additional smoothing factors.

FIG. 3 .
FIG. 3. Identifying the optimum smoothing factor: the blue graph displays the residual correlation as a function of smoothing factor.The red graph displays its derivative.