Unmixing coercivity spectra from backfield curves

Background on coercivity spectra¶

Natural samples contain ferromagnetic materials dispersed in paramagnetic and/or diamagnetic matrix material. Ferromagnetic mineral assemblages can be composed of minerals with different compositions, shapes, and sizes, such that they have different coercivities. Coercivity refers to the resistance of a ferromagnetic material to becoming demagnetized in the presence of an applied field. It is common in the fields of rock and environmental magnetism to use backfield magnetic remanence curves to gain insight into coercivity and thereby the ferromagnetic mineral assemblage within samples. Alternatively, a similar analysis can be conducted on isothermal remanent magnetization (IRM) acquisition or alternating field (AF) demagnetization data.

A backfield experiment is typically conducted on a vibrating sample magnetometer (VSM). In such an experiment, a strong external field is applied to the sample in one direction to induce a saturation isothermal remanent magnetization (sIRM), followed by the application of a field in the opposite direction. The process is repeated with the sIRM being applied, followed by reverse fields that are progressively larger in each step of the experiment. The derivative of the resultant backfield curve is referred to as the coercivity spectrum. A common data analysis workflow is to fit components of this spectrum to ‘unmix’ different populations of magnetic minerals.

Programs developed for this analysis include MAX UnMix (Maxbauer et al., 2016), IRM UnMix (Heslop et al., 2002), GECA (Egli, 2003), and IRM-CLG (Kruiver et al., 2001). This notebook documents the unmixing framework implemented in rockmagpy, which includes an implementation of the MAX UnMix workflow alongside a complementary direct curve-fitting approach.

The component model¶

Robertson and France (1994) made the observation that the coercivity distribution of a single population of magnetic grains can often be closely approximated by a Gaussian distribution in $\log_{10}(B)$ space (a log-Gaussian). Egli (2003) and subsequent work showed that natural components can deviate from perfect log-normality, motivating skewed generalizations.

In rockmagpy, every unmixing component is a skew-normal distribution in $x = \log_{10}(B/\mathrm{mT})$ described by four parameters:

• contribution — the integrated area of the component in the units of the magnetization data (its share of the total remanence);

• location — the location parameter $\xi$ (equal to the mean only for symmetric components);

• dp — the scale parameter $\omega$ in $\log_{10}$ units (equal to the standard deviation only for symmetric components);

• skew — the Azzalini shape parameter $\alpha$ (0 = symmetric log-Gaussian; negative values skew the distribution toward low coercivities).

Because the skew-normal has an analytic cumulative distribution function (via Owen’s T function), the same component model can be fit in two different data spaces, and all methods are accessed through a common entry point, rmag.unmix_coercivity(x, M, method=...):

1. spectrum space (method='spectrum') — fitting the derivative $\mathrm{d}M/\mathrm{d}\log_{10}(B)$, the approach popularized by Kruiver et al. (2001), Egli (2003), and MAX UnMix (Maxbauer et al., 2016). This operates in the space where components are visually interpreted, but requires numerically differentiating (and typically smoothing) the measurements.

2. measurement space (method='curve') — fitting the measured remanence curve $M(B)$ directly with cumulative components, which avoids differentiation and smoothing altogether and works with the original noise structure of the data.

3. MAX UnMix emulation (method='maxunmix') — the spectrum-space fit combined with the MAX UnMix resampling uncertainty scheme of Maxbauer et al. (2016). Two implementation details differ from the R original: skewness follows the Azzalini parameterization rather than the Fernandez–Steel one, and component amplitudes are parameterized by integrated area rather than peak height, which maps directly onto the fraction of remanence carried by each component.

4. Bayesian posterior sampling (method='bayes') — nested sampling of the full posterior in measurement space (requires the optional dynesty package), which represents parameter trade-offs between overlapping components and returns the Bayesian evidence for choosing the number of components. This method, together with multi-start mapping of alternative solutions and the mineral-prior library, is documented extensively in the companion notebook coercivity_unmixing_bayesian.ipynb.

Because the parameterization is shared, results from the different methods are directly comparable, and agreement between the two fitting spaces is a good indication of a robust unmixing model. Derived quantities reported for every component include the mean coercivity B_mean_mT, the median and peak coercivities, the distribution standard deviation sd_log (the dispersion parameter DP of the literature), and the proportion of the total remanence. The method registry is user-extensible (rmag.register_unmixing_method; see the companion notebook).

Import scientific python libraries¶

Run the cell below to import the functions needed for the notebook.

MaxUnmix example¶

We start with the example data file that was published with MaxUnmix. rmag.process_backfield_data shifts the magnetization positive and converts the field values to $\log_{10}$(mT). A light smoothing (smooth_frac=0.075) is applied here for display of the spectrum; the degree of smoothing is a choice the user should test, since over-smoothing biases component widths (Maxbauer et al., 2016).

Interactive choice of initial parameters¶

An important aspect of the MaxUnmix process is for the user to make an initial guess for the parameters associated with the coercivity spectra that can then be optimized to fit the spectra. Initial choices matter, because coercivity unmixing is a non-convex problem in which differently seeded fits can converge to different solutions (demonstrated explicitly near the end of this notebook).

rmag.coercivity_unmixing_interactive provides sliders for each component’s peak field (mT, on a log scale), proportion of the total remanence, dispersion (DP, in $\log_{10}$ units), and (optionally) skew. The sliders are seeded by automatic peak detection, and the live curves show the model implied by the current slider values. Pressing Fit from these values runs the optimizer starting from the slider settings and overlays the result.

For live figure updates run %matplotlib widget first (requires the ipympl package); with the inline backend the widget still works but the figure only refreshes when a cell is re-run.

The returned handle is a dictionary whose 'initial_parameters' entry always reflects the current slider values and whose 'result' entry is filled once Fit has been pressed.

MaxUnMix fit with resampling uncertainties¶

method='maxunmix' fits the spectrum from the chosen initial parameters and then applies the MAX UnMix resampling scheme, faithfully reproducing the Maxbauer et al. (2016) workflow: each replicate recomputes the spectrum from a random 95% subset of the magnetization points (drawn without replacement) and refits it, restarting from the best-fit parameters perturbed by 2% Gaussian noise; parameter means, standard deviations, and 95% percentile intervals over the converged replicates are reported, and the fit plot shows the corresponding uncertainty bands.

Case study: red bed backfield data in MagIC format¶

The rest of this notebook works with backfield curves from hematite-bearing siltstone intraclasts of the ~1 Ga Freda Formation (Bad River, Wisconsin) that were studied in:

Swanson-Hysell, N. L., Fairchild, L. M., & Slotznick, S. P. (2019). Primary and secondary red bed magnetization constrained by fluvial intraclasts. Journal of Geophysical Research: Solid Earth, 124. Swanson‐Hysell et al. (2019)

These specimens contain a narrow high-coercivity detrital hematite component and a broader lower-coercivity pigmentary hematite component, making them a nice unmixing test case. A subset of the study’s MagIC-format rock magnetic data is included in the example_data folder of this repository (the full contribution can be downloaded from MagIC with ipmag.download_magic_from_id). We load the contribution and filter the measurements for experiments with the LP-BCR-BF method code (lab protocol - backfield coercivity remanence - backfield); columns magn_mass and treat_dc_field contain the mass-normalized magnetization and the applied back field respectively.

Processing an example experiment¶

We process the backfield curve for specimen BRIC.22c. Here we leave the data unsmoothed (smooth_frac=0) — the fitting approaches below are designed to work with raw measurements, and smoothing choices are one of the subjective steps that unmixing results should not hinge on.

Unmixing the coercivity spectrum (method='spectrum')¶

With a number of components specified (rather than an explicit initial-parameter table), initial guesses are estimated automatically by peak detection with rmag.estimate_coercivity_components. The fit is a bounded trust-region least squares (via scipy.optimize.least_squares). The returned dictionary contains the fitted parameter table ('params', with linearized standard errors se_* and derived coercivities), the fitted curve, residuals, and fit statistics (RSS, $R^2$, AIC, BIC).

Here we fit two symmetric log-Gaussian components (vary_skew=False). While the library default is to allow skew (vary_skew=True), fixing the skew at zero is often a sensible starting point: skew is best varied only when there is a clear data-driven or mineralogical justification, since skewed components can trade off against one another.

Fitting the measured curve directly (method='curve')¶

Numerical differentiation amplifies measurement noise, and the degree of smoothing applied before differentiating can bias the widths and amplitudes of fitted components (Heslop et al., 2002; Maxbauer et al., 2016). The 'curve' method sidesteps this by fitting the backfield curve itself with cumulative skew-normal components plus a constant offset — no derivative is ever taken. There is empirical and theoretical support for log-Gaussian coercivity components, and since the cumulative model is the exact integral of the spectral model, nothing is lost by fitting in measurement space; the parameters remain those of the log-Gaussian (or skew-normal) coercivity distributions.

Note that for a shifted backfield curve, which spans twice the saturation remanence, each fitted contribution is twice the remanence carried by that component; proportion values are unaffected.

The two data spaces agree closely on the narrow, high-coercivity component (mean coercivity ~720 vs ~705 mT, essentially identical dispersion) but not on the broad, low-coercivity component: its mean coercivity differs substantially between the two fits (~242 mT in spectrum space vs ~425 mT in measurement space), along with its dispersion (DP ~0.57 vs ~0.74) and its share of the remanence (~25% vs ~33%). This is the first concrete instance of a theme that runs through the rest of this notebook and the companion Bayesian notebook: a broad, low-amplitude component is intrinsically poorly constrained, and two nearly equivalent fits can partition it quite differently. Agreement between the two data spaces is a useful robustness check precisely because disagreement flags which components are, and are not, well determined — here the narrow high-coercivity component is trustworthy while the broad low-coercivity one should be treated with caution.

How many components?¶

The number of components is the most consequential modeling choice. rmag.compare_unmixing_models provides the standard statistical tools: information criteria (AIC/BIC, where lower is better) and sequential F-tests against the next-simpler model.

The second component produces an enormous improvement (large F, vanishing p-value, a large drop in AIC/BIC), unambiguously supporting two components over one. But notice that the criteria still improve for a third component — in this fit it is even the AIC/BIC minimum, and the F-test finds it significant (p ≈ 0.0005) — despite the two well-separated mineral populations already being captured. This is a general and important trap: with high-resolution, low-noise curves, AIC, BIC, F-tests, and the Bayesian evidence tend to keep rewarding extra components well beyond the mineralogically meaningful count, because a real coercivity distribution is never exactly log-Gaussian and each added component mops up a little more systematic misfit. (Here the fourth component is finally rejected — it worsens both AIC and BIC and the F-test is insignificant, p ≈ 0.86 — but on even lower-noise data the criteria would keep favoring more.) Left unchecked, these criteria would happily fit spurious components that subdivide a single mineral population.

The fix is an explicit parsimony rule: prefer the simplest model that is already adequate, rather than the one that minimizes a criterion (Egli, 2003; Heslop, 2015). rmag.select_n_components implements this. Its default ‘parsimony’ criterion retains an added component only if it reduces the residual by at least a set fraction (default 2%) of the baseline single-component misfit — so the second component (which removes ~99% of the misfit) is kept while a third (which removes only ~1%, below the 2% threshold) is rejected, regardless of how low the noise is.

The rss_improvement_fraction column makes the choice transparent: the second component explains the overwhelming majority of the baseline misfit, while the third and fourth explain almost none. The reduced_chi2 column shows the complementary view — adding the second component collapses the reduced chi-square by roughly two orders of magnitude; the residuals that remain sit about an order of magnitude above the estimated measurement noise (a real spectrum is never exactly a sum of log-Gaussians), and further components chip away at that small systematic misfit rather than capture additional mineral populations. (For data with a trusted, independently measured noise level, criterion='chi2' selects the simplest model whose reduced chi-square reaches ~1 instead.) The parsimony choice of two components matches the detrital plus pigmentary hematite populations documented for these red beds — a case where the statistics and the mineralogy agree, but where the parsimony rule, not raw criterion-minimization, is what prevents over-fitting.

Uncertainty: linearized standard errors versus the bootstrap¶

Two uncertainty estimates are available:

1. Linearized standard errors (se_* columns in the parameter table) come from the covariance of the least-squares fit. They are fast, but they assume the model is locally linear in its parameters and they degrade badly when components overlap or parameters sit near bounds — conditions that are common in coercivity unmixing.

2. Bootstrap resampling (rmag.unmixing_bootstrap) refits the model to many resampled versions of the data and reports percentile intervals of every parameter and derived quantity. This captures the full nonlinearity of the problem, including asymmetric uncertainties and parameter trade-offs, and is the recommended approach for reporting. (method='maxunmix' used above is this same machinery with the MAX UnMix resampling settings.)

Validation against saved MAX UnMix results¶

Backfield curves for four of the red bed specimens, together with the two-component results saved from the MAX UnMix Shiny app (Bh and DP in $\log_{10}$ units and the observed contribution OC of each component), are included in the example_data folder. We fit each curve with symmetric components seeded automatically and compare against the MAX UnMix output. Exact agreement is not expected — MAX UnMix fits skewed components with a different skew parameterization, its Bh is a distribution location rather than a mean, and its input files were trimmed slightly differently — but the recovered component structure should be consistent.

The symmetric two-component fits recover the same structure as MAX UnMix: a broad lower-coercivity component carrying roughly a quarter to a third of the remanence and a narrow component near 700–830 mT carrying the rest, with dispersions closely matching the saved DP values for the narrow component. The largest differences are in the broad component of specimens where the saved MAX UnMix models used strong skew (e.g., BRIC_20c, S = 1.35): a strongly skewed broad component and a symmetric broad component at lower field describe these spectra almost equally well. This kind of trade-off is inherent to coercivity unmixing and is why consistent initialization across a specimen suite — and reporting uncertainty — matter more than any single best-fit number.

Sensitivity to initial parameter choices¶

Coercivity unmixing is a non-convex optimization problem: differently initialized fits can converge to different local minima that describe the data almost equally well, particularly when skew is allowed to vary. The cell below fits the BRIC.20c spectrum twice with skewed components — once seeded by automatic peak detection and once seeded near the solution published from the MAX UnMix program — and compares the outcomes.

Both fits describe the spectrum nearly equally well, yet they partition the remanence differently — a reminder that with skewed components the decomposition is not unique.

Assessing uncertainty and non-uniqueness¶

The bootstrap intervals above measure repeatability under measurement noise, but they are conditional on the single solution the optimizer found — they do not reveal whether a genuinely different decomposition would fit comparably well, and they cannot incorporate mineralogical prior knowledge. Because coercivity unmixing is an ill-posed problem, quantifying this uncertainty is essential, and rockmagpy provides a dedicated toolkit for it:

• rmag.unmixing_multistart launches the fit from many initializations and reports the distinct solutions the spectrum admits, with their statistical weights;

• method='bayes' samples the full posterior with nested sampling, giving the complete uncertainty representation and the Bayesian evidence for choosing the number of components;

• rmag.mineral_priors constrains the decomposition with a library of characteristic mineral coercivity ranges;

• rmag.plot_unmixing_posterior, plot_unmixing_ensemble, and plot_unmixing_tradeoff visualize the uncertainty.

These methods, and the practical guidance for using them, are documented in depth in the companion notebook coercivity_unmixing_bayesian.ipynb, which also covers extending the framework with custom methods.

Batch processing¶

If a suite of specimens has similar components (both in terms of number of components and their mean/dispersion) then the same initial parameters (or the same automatic seeding strategy) can be used for all of them. rmag.unmix_backfield_experiments processes every LP-BCR-BF experiment in a measurements table: data preparation, fitting with any registered method, optional resampling uncertainties, and collection of the results into a tidy components table (one row per experiment and component). A common initial-parameter table — or a per-experiment dictionary of tables, e.g. assembled with the interactive widget — can be supplied through the initial_parameters argument.

Recording results in the MagIC specimens table¶

rmag.add_unmixing_to_specimens_table records the unmixing results following the MagIC data model: with the default mode='rows', one specimens row is added per component with the controlled-vocabulary columns rem_cmf (component median field, tesla), rem_cd (component dispersion, $\log_{10}$ units), and rem_n_comp (number of components), plus the complete parameter set (proportions, confidence intervals) as JSON in the description column. An alternative mode='description' leaves the table structure unchanged and stores the full model as JSON on the existing specimen rows (readable back with rmag.parse_specimen_description). The resulting table can be written out with pmag.magic_write for upload to MagIC.

Next steps¶

For rigorous uncertainty quantification — mapping non-unique solutions, sampling the full Bayesian posterior, constraining the fit with mineralogical priors, and visualizing all of it — see the companion notebook coercivity_unmixing_bayesian.ipynb. It also shows how to register custom unmixing methods with rmag.register_unmixing_method.

References¶

Egli, R. (2003). Analysis of the field dependence of remanent magnetization curves. Journal of Geophysical Research, 108, 2081.

Heslop, D., Dekkers, M. J., Kruiver, P. P., & van Oorschot, I. H. M. (2002). Analysis of isothermal remanent magnetization acquisition curves using the expectation-maximization algorithm. Geophysical Journal International, 148, 58–64.

Heslop, D. (2015). Numerical strategies for magnetic mineral unmixing. Earth-Science Reviews, 150, 256–284. https://doi.org/10.1016/j.earscirev.2015.07.007

Kruiver, P. P., Dekkers, M. J., & Heslop, D. (2001). Quantification of magnetic coercivity components by the analysis of acquisition curves of isothermal remanent magnetisation. Earth and Planetary Science Letters, 189, 269–276.

Maxbauer, D. P., Feinberg, J. M., & Fox, D. L. (2016). MAX UnMix: A web application for unmixing magnetic coercivity distributions. Computers & Geosciences, 95, 140–145. Maxbauer et al. (2016)

Robertson, D. J., & France, D. E. (1994). Discrimination of remanence-carrying minerals in mixtures, using isothermal remanent magnetisation acquisition curves. Physics of the Earth and Planetary Interiors, 82, 223–234.

Swanson-Hysell, N. L., Fairchild, L. M., & Slotznick, S. P. (2019). Primary and secondary red bed magnetization constrained by fluvial intraclasts. Journal of Geophysical Research: Solid Earth, 124. Swanson‐Hysell et al. (2019)