Hysteresis Processing walkthrough notebook

This notebook provides a step-by-step walk through of the processing of hysteresis loops. Each processing step is illustrated and described. This walk through can enable more granular data processing and illustrates the decision tree that is automatically implemented in the hysteresis_processing.ipynb.

Install and import packages¶

Import local data in MagIC format¶

In this demonstration we will be using local data.

The data is from the following publication:

• Swanson-Hysell, N. L., Avery, M. S., Zhang, Y., Hodgin, E. B., Sherwood, R. J., Apen, F. E., et al. (2021). The paleogeography of Laurentia in its early years: New constraints from the Paleoproterozoic East-Central Minnesota Batholith. Tectonics, 40, e2021TC006751. https://doi.org/10.1029/2021TC006751

Get hysteresis data from the measurements dataframe¶

The method codes relevant to hysteresis loops are:

• LP-HYS for regular hysteresis loops

• LP-HYS-O for hysteresis loops as a function of orientation

• LP-HYS-T for hysteresis loops as a function of temperature

Inspect the original measurement level data¶

• Below we investigate all the available specimens that have hysteresis data in the database for this contribution.

• We extract the experimental data for an individual specimen.

Let’s take a look at the columns relevant to hysteresis data processing

• magn_mass is mass normalized magnetization (in Am^2/kg)

• meas_field_dc is the applied DC field (in T)

Visualize the raw data with plot_hyst_loop¶

• The function returns a Bokeh plot object which can be further customized before displaying

Loop processing¶

rockmagpy calculates the hysteresis loop parameters as available in the IRMDB software. Summary parameters can then be exported into the MagIC specimen data table.

The processing functions implement the protocol of Jackson and Solheid (2010) as realized in the hysteresis processing software of the Institute for Rock Magnetism (IRM). Paterson et al. (2018) present a related treatment of hysteresis data processing (implemented in HystLab) that is a complementary resource. This flowchart shows the decision tree for data processing.

Interpolate the raw data and re-grid them into an upper and a lower branch with exactly the same and symmetric field steps¶

• The applied DC field values through the upper and lower branches during experiments are usually not exactly symmetric about the origin. We need to correct for the slight mismatching between the field steps by splitting the loop into an upper branch and a lower branch, interpolating the branches and unifying their field axis to be the same and be symmetric about 0.

• The function rmag.grid_hyst_loop accomplishes this task.

Perform linearity test on the whole loop¶

• If the whole loop is linear, the sample is dominated by paramagnetic or diamagnetic materials and no ferromagnetic parameters can be estimated (the high-field slope is the result to report).

• If the whole loop is not linear, we should move on to further processing.

• The FNL statistic is a lack-of-fit F ratio from an ANOVA partition in which the upper branch and the inverted lower branch are treated as replicate measurements (Jackson and Solheid, 2010). FNL values above ~1.25 reject whole-loop linearity at 95% confidence for typical loop sizes.

• The resulting stats dictionary from the hyst_linearity_test function has a key loop_is_linear with a boolean value that signals whether the loop is linear or not.

• process_hyst_loop implements this as a decision-tree exit: a statistically linear loop terminates with only χ_HF reported (from the whole-loop regression) and the ferromagnetic parameters as NaN. Passing fit_linear_loop=True overrides the exit and processes the loop in full, which can be useful when a weak ferromagnetic signal near the noise level is of interest.

Loop centering¶

• In the example case, the loop is not linear, so we need to isolate the ferromagnetic component from the paramagnetic component.

• We need to center the loop to remove any offset in the field and moment axes that could be present in the data.

Quality factor Q¶

Once the loop has been centered, the quality factor Q can be calculated. Following Jackson and Solheid (2010), Q is the decimal log of the signal/noise ratio, where the noise is quantified from the mean square mismatch between symmetrically equivalent moments of the offset-corrected loop. The convention used is $Q = \log_{10}(1/\sqrt{1-R^2})$, matching the IRM’s processing software and HystLab (Paterson et al., 2018, eq. 4); note that the corresponding equation as printed in Jackson and Solheid (2010) omits the square root, although the values reported in that paper include it (see Paterson et al., 2018). A higher value indicates better quality data. Note this value can deviate from that currently calculated in the IRM VSM Ferret as that value is estimated from the FNL statistic.

High Q values (e.g., 2 or greater) indicate that deviations from inversion symmetry due to all possible sources (noise, drift, and inherent asymmetry) are small. Conversely, low values indicate that at least one of these is significant. When Q falls below ~0.3 (signal/noise ~2) it becomes difficult to compute meaningful parameters for the loop, and when Q < 2 the automated processing skips the horizontal (field) offset correction because it cannot be robustly determined.

Loop drift correction¶

• Sometimes spurious changes occur in measured signal strength on time scales comparable to that of the loop measurement. Drift is typically manifested by failure of loops to close, by lack of even symmetry in M_rh, and/or by failure of M_rh to decrease to zero in large (positive or negative) fields.

• In many cases a correction is not strictly necessary as the drift signal is small. The default approach is to examine the smoothed error curve M_e (the mismatch between the upper branch and the inverted lower branch) and determine whether a high-field drift correction or an upper branch drift correction is appropriate.

• If you zoom in on the ~1 T region of the loop above, you will find that there is a slight difference between the beginning and end points.

• The Me_drift_correction function corrects the drift by subtracting a smoothed M_e curve either in the high-field region (when the largest errors occur above 75% of the maximum field) or over the entire upper branch of the loop.

Calculate the M_r (remanent magnetization), M_rh (remanence component), M_ih (induced component), M_e (error), B_rh (field corresponding to 1/2 M_r)¶

• The function calc_Mr_Mrh_Mih_Brh calculates these values and curves.

• M_rh = (M⁺ - M^-)/2 and M_ih = (M⁺ + M^-)/2 where M⁺ and M^- are the upper and lower branches; for a well-behaved loop, M_rh is an even function of field and M_ih is odd (Jackson and Solheid, 2010).

• Below we print out the M_r value and B_rh value.

Test loop saturation (high-field linearity test)¶

• Next we test the loop saturation. This test allows us to determine the appropriate approach for fitting the high-field portion of the loop and determining M_s.

• The lack of saturation can cause a loop to deviate from linearity at high fields. Although such a deviation from linearity can be subtle to the eye, it can have a significant effect on estimates of M_s. As developed in Jackson and Solheid (2010), a linearity test is conducted over high-field windows starting at 60%, 70%, and 80% of the maximum field. Each window’s lack-of-fit F ratio (FNL60, FNL70, FNL80) is compared with the critical value of the F distribution at the window’s degrees of freedom (F_critical60 etc., at alpha=0.05 by default, as in Jackson and Solheid, 2010, and HystLab); the p-values are returned as well. If every window rejects linearity, a linear fit to the high-field data would be inappropriate for correcting the paramagnetic or diamagnetic slope and an approach-to-saturation fit is required; otherwise the linear fit is taken over the widest window that is statistically linear (the saturation_cutoff, None when no window is linear). Because the drift correction subtracts a smoothed version of the branch mismatch that the test estimates its pure error from, process_hyst_loop passes the fraction of the pure-error degrees of freedom that the smoothing leaves (pure_error_df_fraction, reported by Me_drift_correction); windows in which the loop was measured at fewer than 12 points are reported as untestable rather than tested on interpolated points.

• In the example case, the loop is not saturated at 60%, 70%, and 80% of the maximum field range. Therefore, an approach to saturation fit needs to be used for isolating the ferromagnetic loop and then calculating the M_s and the high field susceptibility values.

Fit high-field slope and get the intercept as M_s¶

• In this case, the loop is not saturated at high fields and we should use an approach to saturation fit to calculate the M_s value.

• But to show the difference between a linear fit and a saturation fit, we will also perform a linear fit on the high-field data and use it to correct the loop. Since no window is statistically linear here (saturation_cutoff is None), the linear fit is taken over the window from 70% of the peak field, the routine default of most instrument software.

• We use the linear_HF_fit function to estimate the high-field susceptibility (and the M_s intercept in the case where the loop is saturated). Note that the returned chi_HF is the fitted slope converted to SI units (multiplied by $\mu_0 = 4\pi \times 10^{-7}$).

• We then use the resultant chi_HF to correct the loop using the hyst_slope_correction function (which applies the inverse unit conversion internally).

Perform approach-to-saturation fit and estimate a more appropriate M_s¶

• In this example case, a linear fit would underestimate the M_s value.

• A more appropriate M_s value can be estimated using a nonlinear approach-to-saturation fit with the hyst_HF_nonlinear_optimization function.

• By default (fit_type='IRM'), the model used at the IRM is fit:

  $$M = \chi_{HF} H + M_s + a_1 H^{-1} + a_2 H^{-2}$$

  (1)

  Alternative model forms of $M = \chi_{HF} H + M_s + \alpha H^{\beta}$ (Fabian, 2006) are available with fit_type='Fabian' (free $\beta$) and fit_type='Fabian_fixed_beta' ($\beta = -2$).

• These models are linear in their parameters apart from $\beta$, so they are solved by least squares (Jackson and Solheid, 2010, section 7) with the sign constraints $a_1, a_2, \alpha \le 0$ that describe an approach to saturation imposed as bounds; no starting values are involved. For the free-$\beta$ form, $\beta$ is searched on a grid over $[-2, -1]$ with a linear solve at each value, as Jackson and Solheid (2010) and HystLab do; $\beta$ is confined to that interval because the $H^{\beta}$ and constant basis functions become indistinguishable as $\beta \to 0$.

• The high-field fitting problem is ill-conditioned (the basis functions are strongly correlated over a typical window), so the uncertainty of M_s and χ_HF is estimated by the bootstrap of Jackson and Solheid (2010): with n_bootstrap set, the high-field points are resampled with replacement and refit that many times, and the standard deviation (Ms_se, chi_HF_se) and 2.5–97.5% range (Ms_ci95, chi_HF_ci95) of the refits are returned.

• The returned Fnl_lin statistic (Jackson and Solheid, 2010, eq. 21) tests whether the nonlinear fit significantly improves on a linear fit; it is compared with the F distribution at the fit’s degrees of freedom and its p-value is returned as Fnl_lin_p (Jackson and Solheid quote ~3–3.5 as the critical value for typical numbers of points).

Test loop closure¶

• When there are significant high coercivity components in the sample such as hematite or goethite, the hysteresis branches will not be overlapping at the highest applied fields. The unsaturated fraction then biases the high-field fit, and so the M_s and χ_HF estimated from it, which is why the test is applied after the fit: the opening is expressed relative to the fitted M_s.

• The loop_closure_test function quantifies the opening with the statistic $f_{open}$ (HF_Mrh_fraction): the mean of the field-reflection average of M_rh, $[M_{rh}(+B) + M_{rh}(-B)]/2$, over the high-field window (80--99% of the peak field by default), divided by the fitted M_s. M_rh is an even function of field for a well-behaved loop and falls to zero once the branches have merged, so $f_{open}$ is the unswitched irreversible moment within the window as a fraction of the ferromagnetic moment, which is also the approximate relative error the opening induces in M_s: zero for a closed loop, positive for an open one, and negative when the branches cross (usually residual drift). The same quantity relative to M_r (HF_Mrh_fraction_Mr, the fraction of the saturation remanence still unswitched) is returned as well; it does not depend on the fit but overstates the opening for material with a small M_r/M_s. Called directly, loop_closure_test estimates the noise from the odd part of M_rh, $M_{rh}(+B) - M_{rh}(-B)$, which is zero for a symmetric loop and carries only the noise. Because gridding correlates neighboring points and the drift correction’s smoothing adds a slowly varying component to the even part that the odd part does not see, process_hyst_loop instead obtains the standard error by propagation (closure_mean_se_by_noise_propagation): the measurement noise is estimated on the measured branches, white noise of that size is generated at the measured fields, gridded and passed through the same drift correction the loop received, and the scatter of the resulting window means is the standard error (the odd-part value is returned alongside as HF_Mrh_mean_se_odd_part). When the approach-to-saturation fit leaves M_s on its bound — a loop so far from saturation that no model describes its high-field branches — $f_{open}$ is undefined and the opening is sized against the largest moment the loop reaches, M_max (HF_Mrh_fraction_Mmax); whether any high-field M_rh is detectable above the noise is decided first (if none is, the loop is closed), and since M_max includes the paramagnetic moment the sized fraction can only confirm that a loop is open, never declare it closed; the results record which normalization the verdict used (closure_normalization). Because the value depends on the window, the window is returned with it.

• The verdict (closure_state) is 'open' when $f_{open}$ exceeds the openness_tolerance (2% of M_s by default) by more than two standard errors, 'closed' when it is below the tolerance by more than two standard errors, and 'indeterminate' when the data do not resolve the question. The signal-to-noise ratio (SNR) and high-field area ratio (HAR) statistics of HystLab (Paterson et al., 2018, section 4.5) are also returned, and can be selected as the decision rule with criterion='SNR_HAR'; they are not the default because they measure whether an opening is detectable above the noise rather than how large it is.

• Loop closure does not stop processing. process_hyst_loop computes every parameter for every loop, returns closure_state, HF_Mrh_fraction and its standard error alongside them, and prints a -W- line for loops that are open or of unresolved closure. Whether a given opening is acceptable is left to the analyst; when writing results to a MagIC specimens table, add_hyst_stats_to_specimens_table(..., exclude_open=True) withholds the slope-dependent parameters (M_s, B_c, χ_HF) of open loops while keeping the slope-independent M_r and B_rh.

• In the example case, the loop is closed.

Now let’s use the approach-to-saturation fit result to correct the loop.

Plot the ferromagnetic component together with the previous intermediate products¶

• You may click on the legend items to hide the curves

• The preferred final ferromagnetic component is shown in pink

Inspect the differences between applying a linear high field slope fit vs. an approach-to-saturation fit¶

Calculate the quality factor of the ferromagnetic component loop¶

Qf is the quality factor Q recomputed on the slope-corrected (ferromagnetic component) loop (the $Q_f$ of Jackson and Solheid, 2010). The high-field slope correction affects the mean square signal strength but not the noise, so Qf quantifies the signal/noise of the isolated ferromagnetic loop.

Calculate B_c¶

• This calculation is based on the ferromagnetic component loop.

• There will be a difference if you use the loops before high-field slope correction.

• The B_c value is calculated as the average of the lower and upper branch.

• The value shown is in unit of Tesla.

The process_hyst_loop function¶

The process_hyst_loop function applies all of the above steps using the decision tree to process the data. The summary table includes the parameters developed above as well as the hysteresis shape parameter $\sigma_{hys} = \ln[E_{hys}/(4 B_c M_s)]$ of Fabian (2003), where $E_{hys}$ is the area of the ferromagnetic loop. Options include NL_fit=True to force an approach-to-saturation fit, centering_protocol='iterative' for weak ferromagnetic loops superimposed on a strong linear background, fit_linear_loop=True to fully process a loop that passes the whole-loop linearity test (the one decision-tree exit), openness_tolerance to set the level of $f_{open}$ (as a fraction of the fitted M_s) at which a loop is flagged as open, quality_threshold to set the quality factor below which a loop is flagged as low quality, and n_bootstrap (default 1000; 0 to skip) for the number of bootstrap refits behind the uncertainties Ms_se and chi_HF_se of the approach-to-saturation fit. The summary table shows M_s with its standard error, and the results record which high-field fit was used (hf_fit).

References¶

Fabian, K. (2003). Some additional parameters to estimate domain state from isothermal magnetization measurements. Earth and Planetary Science Letters, 213(3-4), 337–345. https://doi.org/10.1016/S0012-821X(03)00329-7

Fabian, K. (2006). Approach to saturation analysis of hysteresis measurements in rock magnetism and evidence for stress dominated magnetic anisotropy in young mid-ocean ridge basalt. Physics of the Earth and Planetary Interiors, 154, 299–307. https://doi.org/10.1016/j.pepi.2005.06.016

Jackson, M., & Solheid, P. (2010). On the quantitative analysis and evaluation of magnetic hysteresis data. Geochemistry, Geophysics, Geosystems, 11(4), Q04Z15. https://doi.org/10.1029/2009GC002932

Paterson, G. A., Zhao, X., Jackson, M., & Heslop, D. (2018). Measuring, processing, and analyzing hysteresis data. Geochemistry, Geophysics, Geosystems, 19(7), 1925–1945. https://doi.org/10.1029/2018GC007620

Swanson-Hysell, N. L., Avery, M. S., Zhang, Y., Hodgin, E. B., Sherwood, R. J., Apen, F. E., et al. (2021). The paleogeography of Laurentia in its early years: New constraints from the Paleoproterozoic East-Central Minnesota Batholith. Tectonics, 40, e2021TC006751. https://doi.org/10.1029/2021TC006751