Unmixing of Magnetic Hysteresis Loops Through a Modified Gamma‐Cauchy Exponential Model

Quantifying the contributions of distinct mineral populations in bulk magnetic experiments greatly enhances the analysis of environmental and rock magnetism studies. Here, we develop a new method of parametric unmixing of susceptibility components in hysteresis loops. Our approach is based on a modified Gamma‐Cauchy exponential model that accounts for variable skewness and kurtosis. The robustness of the model is tested with synthetic curves that examine the effects of noise, sampling, and proximity (similar coercivities) of susceptibility components. We provide a Python‐based script, the Hist‐unmix, which allows the user to adjust a forward model of up to three ferromagnetic components as well as a dia/paramagnetic contribution. Optimization of all the parameters is achieved through least squares fitting (Levenberg‐Marquardt method), with uncertainties of each inverted parameter calculated through a Monte Carlo error propagation approach. For each ferromagnetic component, it is possible to estimate the saturation magnetization (Ms), saturation remanent magnetization (Mrs) and the mean coercivity ( Bc ${B}_{c}$ ). Finally, Hist‐unmix was applied to a set of weakly magnetic carbonate rocks from Brazil, which typically show distorted hysteresis loops (wasp‐waisted and potbellied). For these samples, we resolved two components with distinct coercivities. These results are corroborated by previous experimental data, showing that the lower branch of magnetic hysteresis can be modeled by the presented approach and might offer important mineralogical information for rock magnetic and paleomagnetic studies.


10.1029/2023GC011048
2 of 21 Analysis of IRM curves were very often performed through the symmetrical cumulative log-Gaussian approach of Robertson and France (1994), but they were later better explained by generalized Gaussian functions (Egli, 2003).For magnetic hysteresis, the shape of some curves typically suggests the presence of more than one magnetic component.These include wasp-waisted (constricted middles, near the origin of the coercivity axis), potbellied (spreading middles near the origin and slouching shoulders) and goose-necked (constricted middles and spreading shoulders) (Tauxe et al., 1996).In some cases, these hysteresis shapes have been considered as a fingerprint of some geological processes, such as remagnetization of carbonate rocks (Jackson & Swanson-Hysell, 2012).This evaluation, however, is usually done qualitatively, without quantitative identification and separation of magnetic components.
There are free-access interfaces to deal with magnetic hysteresis data, such as HystLab of Paterson et al. (2018), but the unmixing of distorted curves is not a focus on their work.There are also several ways to unmix magnetic mineral populations from magnetic hysteresis.Some authors model the magnetic properties of natural materials by assuming end members in a mixture, which could be either pure magnetic phases with different grain sizes or typical mineral sources in the study area or region or yet end members identified from the data itself (Heslop & Roberts, 2012;Jackson & Solheid, 2010;Thompson, 1986).Another approach requires the fitting of basis functions to the hysteresis loops.In this case, the linear combination of different basis functions representing the different magnetic populations should represent the bulk behavior of the magnetic assemblage (Heslop, 2015).The advantage of this approach is that it requires little to no a priori information, relying on the ability of a mathematical model to represent a physical phenomenon (Vasquez & Fazzito, 2020;von Dobeneck, 1996).
Many families of basis functions might have potential applications for modeling magnetic hysteresis data (especially those based on Preisach models).Some families include Gaussian and Cauchy distributions (Mörée & Leijon, 2023;Pruksanubal et al., 2002Pruksanubal et al., , 2006;;Stancu et al., 2005); hyperbolic functions (von Dobeneck, 1996); and sigmoid logistic functions (Jackson & Solheid, 2010).Cauchy distributions have many applications in mechanical and electrical theory, often referred to as Lorentzian distributions in the physics literature.
A simple solution for the unmixing of magnetic components by fitting Lorentzian curves to the lower branch of magnetic hysteresis loops has been proposed by Vasquez and Fazzito (2020).It considers the magnetization (M) acquired through the induction of an applied field (   ) as expressed by: The first term of the Equation 1 describes a linear magnetization acquired through an inducing field   , which is the dia/paramagnetic contribution to   (𝐵𝐵) .Consequently, the second (and non-linear) term represents the ferromagnetic contribution, while A is the total area under the  () curve.If   is equal to   , the ferromagnetic contribution will be zero, which is the very definition of coercive force.If   approaches the infinity, Equation 1will tend to  ∕2 , which is the magnetization saturation  () of   (𝐵𝐵) .Now, if Equation 1 is evaluated at zero field (  = 0 ), then saturation remanence  (rs) is also easily calculated.The magnetic susceptibility (κ) is sequentially computed as follows: In order to model the susceptibility components, one of the branches of a magnetic hysteresis is used to calculate a numerical derivative.Vasquez and Fazzito (2020) fitted the parameters of Equation 2 using a generic inversion routine through commercial and/or free-software and report coherent results in the unmixing of components from previously published data (Roberts et al., 1995) and from their own synthetic samples, but acknowledge that the simplicity of the model might fail to cover more complex scenarios.Such a case could arise from the contribution of fine SD-like particles (e.g., a Stoner-Wohlfarth assemblage- Stoner & Wohlfarth, 1991).A distribution of such grains might cause the reversible and irreversible segments of a lower branched magnetic hysteresis to be very different, which will result in an asymmetry.Furthermore, for viscous SD-like particles, the irreversible segment may abruptly start at  = 0 , leading to a discontinuous derivative (Egli, 2021).Neither of these cases can be explained by a symmetrical Lorentzian curve of the form of Equation 2 and would require a skewness control parameter, similar to the coercivity analysis of Egli (2003).It is also important to consider that Equation 2Validation: U. D. does not account for the approach-to-saturation behavior expected in high-fields (Fabian, 2006); therefore, an additional parameter is required to account for a variable kurtosis and susceptibility components with different tails.In situations where distorted hysteresis loops are observed, it is advantageous to utilize a model capable of accommodating variable curves with flexible morphology.This flexibility in the model allows for an accurate representation of the susceptibility component(s) involved.The presence of distinct populations of grains with different coercivities can give rise to various symmetries, resulting in characteristic variations like wasp-waists and potbellies within the hysteresis loops.Therefore, incorporating a model that considers variable kurtosis/ skewness may prove to be beneficial in capturing these distinctive features.
To achieve a more robust phenomenological model to unmix susceptibility components from magnetic hysteresis data, we introduce the use of generalized gamma-Cauchy exponential distributions (Alzaatreh et al., 2016).We present a Python-based (ipynb-file) open-source script (Hist-unmix) that can be used to perform unmixing of hysteresis curves (Bellon et al., 2023).A forward model of up to three susceptibility components is demonstrated, as well as the mathematical formulation to optimize initial parameters in our inverse model, with uncertainty estimates of the parameters determined through a Monte-Carlo error propagation.We also perform numerical tests on synthetic data to assess the sensitivity of a modified Gamma-Cauchy Exponential fit (mGC), evaluating the effect of (a) sampling, (b) signal/noise ratio, (c) similarity of components and the (d) ambiguity of the model.Finally, we test the Hist-unmix script on distorted hysteresis loops of Neoproterozoic remagnetized rocks from São Francisco craton (Brazil), comparing the information recovered from Hist-unmix with previous rock-magnetism/paleomagnetic data.

Forward Model
To achieve a forward model for the first derivative of a lower branched magnetic hysteresis, we propose the use of the probability density function of a gamma-Cauchy exponential distribution (GC(α, β, θ)).If a random variable follows a gamma distribution with parameters α and β, a GC(α, β, θ)'s probability density function is defined as (Alzaatreh et al., 2016): In Equation 3,   has the role of a dispersion parameter (such as in the symmetrical Lorentzian functions) and  Γ() is the gamma function of   .The advantage of using GC (α, β, θ) functions is that their morphology can be symmetrical, right or left skewed, and cover a wide range of kurtosis (Alzaatreh et al., 2016).Since Equation 3will peak in the mode of   (which will coincide with its mean), we added a term to represent the coercivity  () in a gamma-Cauchy distribution.To improve convergence, a scale factor (I) is further included, which represents the contribution ratio of each ferromagnetic component.Our modified gamma-Cauchy exponential function,  GC(, , , , ) for magnetic susceptibility becomes: to work directly with the ferromagnetic contribution, a forward model is then simply given as:

Inverse Model
Whilst we have arbitrarily chosen to model the lower branch, it is of course assumed that the lower and upper branches are symmetrical and centered.If not, some preprocessing must be performed to achieve more coherent results.) are determined using the Levenberg-Marquardt method (Aster et al., 2013;Gavin, 2022), as: where is the Jacobian matrix of   (  (0) + ∆(0) ; I is an identity matrix with the same dimensions as ⋅ ;  (0) is a damping factor and  ∆(0) is calculated as: Where the first iteration begins by adjusting the parameters so that (1) =  (0) + ∆(0) .Obtaining analytically might result in singular matrices, which is a problem that can be avoided when these derivatives (  ∕ ) are here computed by making small adjustments (   ) to each parameter, and evaluating their effect through a numerical central difference finite approach.We define a correction criterion  () in order to evaluate if the adjusted parameters   (+1) better explain the observed model   than   () : 1. is updated using the corrected parameters   (+2) =  (+1) In the criteria above,   is the damping factor that will be updated by each step-scaling factor   .Both of these start with the same initial value of 0.1, as in the fixed approach of Hagan and Menhaj (1994).Iterations (i) will proceed until a convergence criterion is reached: If the user has previous knowledge of the coercivity component values in the sample (i.e., from other magnetic experiments), it might be useful to constrain these   values.When dealing with more than one component, the user might constrain one of the two coercivities and let the other be optimized (or even constrain them all, if required).Care in this approach is required since the model may produce biased results due to the constraints.Inverting a component with   = 0 (i.e., a superparamagnetic population) might also cause numerical issues when calculating the Jacobian matrix (such as singular matrices), so it is useful to constrain the coercivity to zero when minimizing the parameters of a superparamagnetic component.
The inverted parameters acquired from the optimization routine can be used to calculate a new forward model.Besides calculating the coefficient of determination  (  2 ) , we can test if this final model (which is a 1D-array) is statistically similar to the observed data.For that we can use a Two-Tailed F-test, considering a null hypothesis that the variance of the data and the variance of the calculated model ) can be distinguished at a 95% confidence interval.

Monte Carlo Error Propagation
With the considerable number of model parameters related to each ferromagnetic component, it is useful to simulate a collection of perturbed solutions to evaluate the statistical confidence of the model solutions.In our approach, we use a Monte Carlo error propagation method (Aster et al., 2013).We assume that our final inverted model produces parameters and then compare its difference with   inv by calculating an empirical covariance estimate: Where q is the number of parameters.Finally, the 95% confidence interval of   inv is computed as (Aster et al., 2013):

Workflow
Figure 1 shows the general workflow for the Hist-unmix package.The first step comprises the filtering of the lower branch of the hysteresis loop.We note that numerical derivatives through finite-differences method are strongly affected by noise, in a way that even small disturbances can cause large spikes.To reduce these effects, we apply a simple moving average  ( filter to the lower branch hysteresis curve: where (L) is the interval used to calculate the mean.This value will depend, logically, on choices made by the user and the number of data points.When applying polynomial or Gaussian smoothing techniques to the entire data set, individual data points can have a substantial impact on the smoothing result, and the presence of outliers, caused by random errors, can easily affect them.Given that this paper introduces an unmixing technique based on a modified Gamma-Cauchy model, we choose not to generate a smoothing curve using models that assume a different data distribution.Instead, we use a simple low-pass filter, as described in Equation 13, as it avoids introducing potential biases by solely relying on the actual distribution of experimental points.The original data and the smoothed curves are sequentially shown, so the user can visually inspect if the filtering step was effective.
For more complex and careful smoothing/preprocessing analysis, we recommend HystLab (Paterson et al., 2018).
The para/diamagnetic component  0 is sequentially estimated from a linear regression of the high-field irreversible section of the smoothed lower branch.The gradient of the smoothed curve is normalized by its maximum value  ( ) and subtracted from  0∕ to facilitate the adjustment of the curves.Sequentially, the user should choose how many ferromagnetic components (C) will be fit to the data.
Path 1 in Figure 1 requires estimation of a forward model, by providing the mean coercivity  () , the deviation (   ), the parameters α and β, and the scale factor (I).The coercivity  () must be specified within the values of the applied field (   ), while   of most of the curves will vary from zero to one (  GC functions, however, allow larger values to be tested).
GC functions can yield a large range of α and β values, but we set their initial input equal to 1 (a symmetrical approach).  parameter will normalize the contribution of the different components and its first estimation is performed automatically when the user selects the number of components.Path 2 determines a straightforward inverse model where the user provides initial guesses without first adjusting the forward model.A Monte Carlo error propagation (Aster et al., 2013) is performed to obtain the covariance of the inverted parameters and their 95% confidence interval, as well as the coefficient of determination (   2 ) and F-test.This consists in taking the final model produced by the inversion routine and carrying on the calculation of a new set of    .To calculate    , we disturb the final inverted parameters by adding random noise drawn from a normal distribution again.However, these random normal distributions will now have a mean centered on each parameter and their standard deviation will now follow a reduced chi-squared statistic of the final inverted model Where q is the number of parameters.

Magnetization Saturation (M s ) and Saturation Remanent Magnetization (M rs )
To calculate M s and M rs we rely on the definite integral of the susceptibility we can approximate M s and M rs of a given ferromagnetic component   through a numerical integration using Simpson's rule (Otto & Denier, 2005) as: Where   + is the maximum positive applied field.Because the quality of numerical integration strongly depends on the horizontal spacing (dB), a one-dimensional cubic interpolation is applied to the gradient data prior to the application of Equations 15 and 16.
The maximum field applied during a hysteresis procedure may be insufficient to induce saturation magnetization.The magnetization in high-fields  (ℎ ) can be expressed as (Fabian, 2006): where   and  Φ are negative constants (called alpha and beta in Fabian's work), for which: (a)  Φ = −2 in homogeneously magnetized defect free-materials; (b)  Φ = −1 for superparamagnetic particles (SP); and (c)  −1 < Φ ≤ 0 for assemblages of particles with closely spaced defects (Fabian, 2006 and references therein).Susceptibility components of Equation 2 are classified as  Φ = −1 curves, which is not appropriate for most of the natural samples.If the maximum applied field is enough to achieve an approach to saturation regime,  Φ must be smaller than zero (Fabian, 2006).As Equation 4 results in ferromagnetic susceptibility components, we remove the induced magnetization of the dia/paramagnetic contribution of Equation 17 (0 ⋅ ) to obtain the high-field ferromagnetic susceptibility  (ℎ ) as the derivative: To obtain   and  Φ , we can follow the same inversion routine described in Section 2.2 by simply changing the susceptibility terms of Equation 6.For example, we calculate a synthetic model based on Equation 18, while considering an applied field from 0.6 to 7T and  Φ = −2 and  = −2.6 (N = 100).These parameters are similar to those modeled in an example given by Fabian (2006), where it is experimentally observed that magnetization reaches saturation near 5T.By using Equation 18, we observe the same as  ℎ tends to zero in the same field values (Figure 2a).In our inversion procedure,  Φ and   converge to the same values either for a model with the whole curve (100 points), or, limiting the field values between 0.6 and 1T (N = 7, Figure 2b), showing that the lower field values within the approach to saturation regime may control these parameters.
Nevertheless, if one decides to use this approach in the observed data, noise might decrease the effectiveness of the optimization of   and  Φ .However, as we apply this high-susceptibility validation test to the unmixed components obtained from Equation 6, that is not an overall issue.For a given ferromagnetic component, if  Φ < 0 , we consider the   obtained from Equation 15 a valid saturation magnetization.If not, we can correct it using the inverted   and  Φ parameters.

Model Sensitivity
We tested the sensitivity of our model using a series of synthetic curves.Five base curves were generated (C 1 to C 5 in Figure 3a) with distinct parameters (Table 1), as well as a number of bimodal combinations, each with 1,000 field values (   ) between -1T and 1T.Coercivity values were simulated within known ranges of typical magnetic minerals (O'Reilly, 1984).We have varied α, β,   , and I to produce curves with distinct tails and symmetry.Since these parameters represent only ferromagnetic components, we neglect the dia/paramagnetic slope (  0 ).
Normally distributed random numbers (   = 0.0 Am 2 ,  = ±5 ⋅ 10 −6 Am 2 ) were added to the synthetic curves, to simulate measurement noise.Measurement errors might vary according to the measurement routine, the sensitivity of the equipment as well as the intensity of the magnetization.First, we optimized parameters of the synthetic models with one ferromagnetic component following Path 1 (Figure 1), and sequentially did the same for the bimodal curves as well.For the latter, we have added a small dia/paramagnetic component  (0).For both cases, the inversion approach produced optimized parameters whose forward model produced coefficients of determination R 2 greater than 0.9 (Table 2, and Figure 4) and indistinguishable variances at 95% confidence (Two-tailed F-test).Inversion of  0 for the unimodal curves return non-zero values, but their magnitude compared to the ferromagnetic susceptibility is negligible.
For the bimodal models (the curves with two ferromagnetic components), inverted curves successfully represent the synthetic data as well.The dia/paramagnetic contribution for the high-field irreversible segment explain well the displacement of the base level either for a strong paramagnetic (e.g., coming from a fabric enriched in biotite) or diamagnetic influences (e.g., coming from a calcium carbonate matrix).
To further test our model sensitivity, we examined the influence of the (a) signal-noise ratio, (b) sampling of the hysteresis curves, (c) the level of contribution to the total magnetic susceptibility and the proximity and dispersion of the components.
Since the data used to fit the  GC functions are the gradient of the magnetization, small perturbations might strongly affect the dispersion.In order to test the sensitivity of the models to the proximity of different magnetic  2), where the two components are so close that susceptibility appears as a single peak.
In this case, even curves with a high signal/noise ratio (  ≈ 0.95) can lead to a high dispersion (compare   -values in A and B scenarios, Figure 5a).However, a moving average filter seems to be very effective to remove random noise, in a way that simply choosing the L-value of five (L = 5, Equation 13) resulted in a good fit, with   2 > 0.9, Figure 3. Synthetic models produced using Equation 4. In the case of a single ferromagnetic component (a), the dia/ paramagnetic slope was set to zero (check Table 1).Further examples are linear combinations of these to produce curves with two (b) and three (c) components.Random noise was added to all the curves to represent measurement errors. 10.1029/2023GC011048 9 of 21 although the error of the less noisy data is smaller.We used the same  1 + 3 case to investigate if the two components could still be detected when reducing the sample size (n) from 1,000 points to 500 points and then to 200 points (Figure 5b).The errors increase as the number of points decrease, even though the inversion procedure satisfactorily recovered the parameters in all cases (   2 >0.9, Figures 5a and 5b).
For the  1 + 3 case, the parameters are very distinct.However, in mixing cases like  1 + 2 (Figures 5c and 5d) where there is an overlapping of distributions with similar parameters, the ambiguity of the model would allow other solutions with similar residuals.This is a recurrent problem that arises with basis function solutions to the unmixing problem, and that also affects generalized Gaussian approaches for IRM unmixing (Egli, 2003;Maxbauer et al., 2016).In our case, constraining the coercivity of the  2 component allowed us to obtain good estimates of the two distributions with small residuals in the sensitivity test for noise similar to that obtained for the  1 + 3 mixture.However, without a priori information (acquired from other kind of magnetic experiment) the constraining of the coercivity value could not be justifiable.Otherwise, we would recommend the simplest model to explain the observed data.Similar issues are seen as we increase the number of components in the sample, exemplified by the two cases shown in Figure 3c.In the case of the  1 + 3 + 5 mixture, the resulting morphology of the curve allows a clear distinction of at least three components and inversion of  , , and  curves result in a fitting with indistinguishable parameters of those that form the original data (Figure 6a).
For the  3 + 4 + 5 case, the mixing of the most coercive fractions produces a broad peak.Since the position of the component of smaller coercivity is more evident, one could adjust two other components to explain the rest of the spectrum (Figure 6b) with an almost negligible residual.However, it is also possible to explain the same curve with a composition of only two components (Figure 6c) with similar quality of fit.Still in this case, increasing the number of components to three (considering   component fixed) will limit the coercivity of the other two components to a single minimum region (Figure 6b)'.However, the objective function of the  3 + 4 + 5 case with only two components (fixing the other parameters) shows that local minima might be present (Figure 6c)'.Still, our procedure to calculate a    vector (revisit Section 2.4) allowed us to avoid the local minimum in Figure 6c'.Nevertheless, assuming that more than two components explain the susceptibility data should only be considered in cases where a priori information is available, or if the shape of the curve clearly indicates their respective contributions.
From modeling (with both two and three components) we observe that there are many factors that could influence the minimum contribution that a component should have so it can be recovered by the numerical model.First, we estimate that the intensity of the smaller component should not be less than two orders of magnitude lower than the largest component.However, it will also depend on (a) the standard deviation (θ) of such population, since lower standard deviations will generate sharper peaks, making it easier to identify a signature in the magnetic hysteresis and (b) the contribution of random noise, since it could significantly blur the signature of a weak component, potentially making its detection impossible.
Finally, we evaluate the presence of SP as one of the susceptibility components.As shown by Tauxe et al. (1996), potbellied and wasp-waisted magnetic hysteresis can be generated by mixing SP and stable SD particles.To examine this, we construct a ferromagnetic mixture as the sum of an assemblage of SP  ( = 0T) with a higher coercive fraction (i.e., SD-magnetite,   = 0.07T ), and another one with a ferromagnetic low coercive fraction (i.e., MD-magnetite,   = 0.002T ), all with the same dispersion.This is the most extreme scenario, since Note.Coercivities   (T) ranging within known values for terrestrial magnetic minerals.

Table 1 Synthetic Ferromagnetic Components (C)
reproducing the same parameters only by varying the coercivity will make the identification of a superparamagnetic fraction a hard task because the difference in coercivity is very small.
We can evaluate changes in the curves with two components by varying their contributions (by adjusting I) to the final synthetic curve.As the contribution of C SD increases, the SP particles become less significant (Figure 7a) but one can still identify that such curve is not perfectly matching the purely SP component.
The same is valid if C SP is mixed with the lower coercivity component in the same proportions (Figure 7c), but in this case, it becomes hard to distinguish the SP component even if its contribution is equal to that of C MD .
When we calculate the second derivative of the lower branch of these hysteresis curves, this observation becomes even clearer.For C SP + C SD mixing cases, the derivative curve will not cross at zero field (Figure 7b), indicating the presence of a magnetic population with larger coercivity.Meanwhile, because C SP and C MD components coercivities are very close, the second derivative of their mixture crosses zero much closer to the origin (Figure 7d).Nevertheless, if there is a priori information of the presence of SP particles then constraining one component to have zero coercivity enhances the identification of the remaining fractions.
Magnetic properties of Sete Lagoas and Salitre formations are very similar (D'Agrella-filho et al., 2000;Trindade et al., 2004): (a) wasp-waisted/potbellied magnetic hysteresis, (b) contradictory Lowrie-Fuller/ Cisowski tests (Cisowski, 1981;Jackson, 1990), (c) anomalously high hysteresis ratios, and (d) tri-axial thermal demagnetization (Lowrie tests) with similarly behaved components.Although these formations belong to different basins and their sampling sites are separated by almost 600 km, they bear very similar paleomagnetic directions.Thermal demagnetization of these samples commonly yields up to three components (A, B and C) with very similar unblocking intervals (Figures 8a and 8e).
Each magnetic component may be correlated to a particular mineral assemblage depicted in the Lowrie test.The Lowrie test consists of the stepwise thermal demagnetization of three IRM acquisitions along three orthogonal axes: hard (1.3 T), intermediate (0.3 T) and soft (0.1 T).Samples from both Sete Lagoas and Salitre formations show similar behavior in these diagrams (Figures 8d and 8h).
The soft component shows sluggish decay up to 400°C, a common behavior for  2. multidomain magnetite.However, there is a steep decay of the soft component at 500°C, probably associated to the C-component of the thermal demagnetization which can be attributed to stable PSD/SD magnetite.Contrastingly, medium and hard components of the Lowrie test are stable up to 250°C (Figure 8d) and rapidly decay at 320°C.This is close to the Curie temperature of monoclinic pyrrhotite.This mineral may correspond to the B-component identified in the Sete Lagoas and Salitre formations.
The magnetic signature of these carbonates is interpreted, as suggested from Pb isotopic data (D' Agrella-filho et al., 2000;Trindade et al., 2004), as a result of a large-scale remagnetization throughout the São Francisco Craton, caused by the percolation of orogenic fluids during the final stages of the Gondwana assembly.Therefore, the B and C-components of both basins would be contemporary and the result of craton-wide chemical remagnetization.The fact that these rocks present more than one stable direction, likely carried by different magnetic minerals with contrasting magnetic properties, makes them an interesting case study to apply the Hist-unmix script.In this section, we have selected samples from each of these formations (Sete Lagoas and Salitre) and measured magnetic hysteresis curves to test our numerical model., where constraining the coercivity of one of the components using a priori information will produce very similar models to the observed data.

Experimental Methodology
Eight samples of the Sete-Lagoas (BB) and Salitre (IR) formations (each) were selected.First, small fragments (  ≈ 1 cm³) were cut from the typical cylindric samples used in paleomagnetic investigations, using a non-magnetic saw.Then, each sample was bathed-in an acid solution (HCl, 10%) for about 5 s to get rid of any superficial contamination, put into an ultrasonic bath (20 min) with ultra-pure water to neutralize any remaining reaction and/or get rid of surface impurities.Samples were consecutively dried in a silica desiccator (at 25°C).A precision balance was used to measure the mass of the samples in order to normalize the subsequent magnetic measurements.
Magnetic hysteresis was performed with a vibrating sample magnetometer (MicroMag 3900 Series VSM), using the discrete measurement mode from -1T to 1T, so the lower branch of each sample totalized 500 data points.
Processing followed the steps provided in Section 2.4 (Path 1), not constraining the coercivity for any of the curves and allowing 300 models  (   ) to run for each of the hysteresis loops.

Modeling With Hist-Unmix
Data from both the Sete Lagoas and Salitre formations have typical hysteresis of multiple mixed components.Samples from Sete Lagoas have constricted middles (wasp-waisted, Figures 9a and 9b) while Salitre samples have spreading middles (potbellies, Figures 9c and 9d).It is worth noting that although these are carbonate rocks, the paramagnetic contribution completely overcomes the diamagnetic response of calcite and dolomite.This paramagnetic contribution (Figure 9e) is probably caused by the presence of terrigenous (essentially Fe-bearing clay-minerals) in these rocks.To avoid any bias, the lower branches of the hysteresis curves were smoothed using small L-values (Equation 13, L < 5).None of the samples could be simply fitted by a single susceptibility component without inducing large errors.The models were calculated assuming two magnetic components (e.g., Figures 11a and 11b) and resulted in   2 > 0.98 with indistinguishable variances from a two-tailed F-test.
Boxplot distributions compiling the results of the inversions are shown in Figure 10.Both Sete Lagoas and Salitre samples show magnetic components with distinct coercivities.For the Sete Lagoas formation, the component with the lowest coercivity (C a ) has a median  ≈ 1.7 mT, with minimum and maximum values of  ≈ 1.0 and 11.0 mT (Figure 10a), with an asymmetric distribution.For the component with the highest coercivity (C b ), the median is 50 mT, with maximum and minimum values of 260 and 15 mT, respectively (Figure 10a).Saturation magnetization (M s , Figure 10b) is similar for both components, which implies that they contribute almost equally to the whole susceptibility spectrum.The shape of the susceptibility curves, however, are quite distinct.C a components have a small dispersion (θ), being constricted to the region around the median, while C b components have greater dispersion, spreading throughout a wide range of coercivities.For Salitre formation samples, the C a components also have an asymmetric distribution, with median coercivity value of  ≈ 0.6 mT and minimum and maximum When employing a superparamagnetic (SP) fraction with the same properties as those of the SD and MD fractions (only varying   ), it becomes difficult to distinguish the SP contribution in both cases.Constraining the coercivity of one of the components to zero allows the user to test if (mathematically) a SP population could explain part of the observations.For the SP populations,   is fixed at 1. values  ≈ 0.098 and 11 mT, respectively (Figure 10d).Bulk coercivities of C b components are mostly higher than those of the Sete Lagoas samples.Minimum and maximum values are  ≈ 95 and 244 mT, respectively and the median is 200 mT (Figure 10d).
For both Sete Lagoas and Salitre formations, coercivity boxplots of C a are quite short and match the expected values for magnetite.We suspect that the smallest coercivity values may also arise from a population of superparamagnetic grains.Although the C b component could be related to more than one high coercivity mineral, such as hematite or pyrrhotite, the contribution to remanence is comparable to or higher than that of C a (Figures 10c  and 10f).Since the remanence of hematite is much smaller than that of magnetite, it must exceed 95 wt% of the magnetic population to influence the magnetic parameters of an assemblage formed by the hematite + magnetite mixing (Frank & Nowaczyk, 2008).Such a high proportion of hematite in these samples would contradict previously published thermal demagnetization data (Figures 8a and 8b) as well as the Lowrie tests shown in Figures 8d  and 8h.This implies that the higher coercivity phase is likely to be monoclinic pyrrhotite.
Most of the modeled curves did not yield a significant asymmetry, so that a simple Lorentzian model (such as those from Vasquez & Fazzito, 2020) could have successfully explained the observed data as well.Nevertheless, some curves (e.g., Figure 11a) might require a more complex model that accounts for distinct degrees of kurtosis and skewness, which is better accommodated by the modified gamma-Cauchy exponential function.
We can evaluate the C a and C b components in terms of their M rs /M s ratios, which are consistent with the expected range between the SD and MD structures (Day et al., 1977;D. J. Dunlop, 2002).Smaller grain sizes tend to have higher M rs /M s ratios.The C a component (whose M rs /M s ratios are below 0.2 and are greater than 0.02) would correspond to larger grain sizes within the PSD threshold (the yet poorly understood vortex state).However, it could also be represented by a mixture of MD + SP particles.The M rs /M s ratios of both components vary broadly because of the authigenic origin of these particles.The compositional heterogeneities in the sedimentary column affects how much iron is available within a region.This might imply different sizes of particles in different locations (depending on how fast the chemical reactions occur and the thermodynamic favorability of their growth).
If C a is a mixture of MD + SP particles, the presence of coarser grains (MD) is supported by the small   values modeled for this component, which could explain the viscous component observed in the thermal demagnetization procedures (Component A, Figures 8a and 8e).
Cb component (whose M rs /M s ratios are usually greater than 0.2) would correspond to either a mixture of SP + SD particles (following the SP + SD mixing trends) or could represent a population with a mixture between equidimensional SD particles + the smallest particles in the PSD range.Therefore, the assemblage of particles forming the C b component is probably the most stable carries of remanence in these carbonate rocks.Some of the M rs /M s ratios of the C b component exceed the 0.5 threshold.In non-equidimensional grains, where the magnetization is strongly controlled by uniaxial shape anisotropy, the M rs /M s ratio for an SD particle is 0.5.But in equidimensional particles, whose magnetization is controlled by magnetocrystalline anisotropy, the M rs /M s ratio can be significantly higher (e.g., 0.866 for magnetite -Dunlop, 2002).
The hysteresis properties of remagnetized carbonate rocks usually plot along a power law trend controlled by cubic magnetocrystalline anisotropy (Jackson & Swanson-Hysell, 2012).This behavior was originally attributed to an authigenic origin for magnetite, resulting in equidimensional grains lacking significant shape anisotropy (Jackson, 1990).Jackson and Swanson-Hysell (2012) have shown, however, that such an interpretation is not necessarily correct.They attribute M rs /M s ratios above the 0.5 threshold in the previous work of Jackson (1990)    as experimental bias caused by a maximum applied field not being enough to saturate the samples (which was around 0.3 T in most of the samples) and experimentally show that shape anisotropy was actually dominant in their remagnetized carbonate samples.Furthermore, these power law trends (when below the 0.5 threshold) might also match the SD + SP mixture trends (as compared with Dunlop, 2002).However, in our work, we apply a maximum field of 1T and provide a high-field saturation test following Fabian (2006) to attest that both C a and C b components are saturated in our maximum applied field.Euhedral and spheroidal iron oxides have been detected in our samples through previous SEM-EDS studies (D' Agrella-filho et al., 2000), so we suggest that a considerable amount of these could indeed contribute to the anomalous M rs /M s ratios calculated for the C b component.
The magnetic data suggest that the major cause of the distorted hysteresis loops in the Sete Lagoas and Salitre formations are populations of magnetic minerals with distinct coercivities.These different populations can be different magnetic minerals, for example, magnetite and pyrrhotite, or different grain sizes of the same magnetic mineral.Tauxe et al. (1996) have shown (numerically) that when related to SP + SD mixtures (of the same magnetic mineral), wasp-waisting requires quick saturation of the SP contribution at low fields.Meanwhile, potbelling requires a low initial slope and saturation of the SP fraction at higher fields.For instance, the high-frequency dependence of magnetic susceptibility reported by previous works suggests that SP particles likely contribute to the magnetic mineralogy of these carbonates.But as argued in Section 3, the hysteresis loops are affected only when the fraction of SP is significantly high, which might be the case for part of the C a components with the lowest coercivity values.Nevertheless, the distribution of coercivities of C b in the Salitre formation is more skewed to higher coercivities (median = 200 mT, Figure 10a), while in the Sete Lagoas it is skewed to the lower coercivities (median = 50 mT, Figure 10b).Although there is no conclusive evidence, this difference in the distribution of sizes within the C b component might explain what is causing the prevalence of wasp-waisted in one formation and potbellies in the other one (Sete Lagoas and Salitre formations, respectively).
An important clue to understanding the remagnetization in these carbonate rocks comes from further information obtained from modeling with Hist-unmix: the significant paramagnetic component present in samples from both the Sete-Lagoas and Salitre formations, which surpasses the ferromagnetic contribution.This paramagnetic contribution is likely due to the high content of clay minerals in these rocks (Callaway & McAtee, 1985;Potter et al., 2004).Clay transformations (smectite-to-illite) are known to release Fe-ions in the medium, which might allow the growth of authigenic ferromagnetic phases (Katz et al., 1998;Tohver et al., 2008) responsible for chemical remagnetization.Therefore, investigating the origin of this large paramagnetic response might help to better constrain the geological processes responsible for the large-scale remagnetization in these two basins of the São Francisco Craton.

Conclusions
We have presented a Python-based open-source code to perform a parametric unmixing of magnetization curves, in order to separate susceptibility components of hysteresis loops.Our phenomenological model is based on a modified gamma-Cauchy exponential function, whose advantage lies in their capacity to explain variable shapes, including symmetrical, right or left skewed curves, and covering a wide range of kurtosis.
Hist-unmix is an easy to use Python script that includes a pre-processing interface, where the lower branch hysteresis data is filtered using a simple moving average.Forward models allow the user to adjust up to three ferromagnetic components and estimate dia/paramagnetic contributions.The parameters controlling each component can be subsequently optimized through a Levenberg-Marquardt method.The mean coercivity of ferromagnetic components can be fixed using a priori information, in order to constrain the solutions.Uncertainty of each optimized parameter is estimated for the final inverse model using a Monte Carlo error propagation (following the reduced chi-squared statistic of the inversion procedure) and its variance is compared to the observed data in order to verify if they are distinguishable at 95% confidence level (Two-tailed F-test).We also implement a test to verify (and correct, if necessary) the magnetization saturation values of each component by modifying the high-field saturation approach of Fabian (2006).
Our numerical routine was applied to separate susceptibility components from wasp-waisted and potbellied curves from Neoproterozoic remagnetized carbonate rocks.The inversion results clearly distinguished two ferromagnetic components: a less coercive (C a ) and a more coercive (C b ) one.Together with the results of Lowrie The experimental procedures in this paper were performed at USPMag lab at Instituto de Astronomia, Geofísica e Ciências Atmosféricas (IAG) at Universidade de São Paulo (USP).This study was funded in part by the Coordenação de Aperfeiçoamento de Pessoal de Nível Superior-Brasil (CAPES)-Finance Code 001, and by Grants 21/00861-2 and 16/06114-6, São Paulo Research Foundation (FAPESP).W.W. would like to acknowledge support from the Natural Environmental Research Council through Grants NE/V001233/1 and NE/ S011978/1.We would like to acknowledge the valuable discussion with Dr. Ramon Egli and Dr. Karl Fabian, which greatly improved this work.Lastly, we thank Dr. David Heslop and Dr. Mark Dekkers for their critic review.

Figure 2 .
Figure 2. Synthetic high-field susceptibility curves.(a) The inversion procedure effectively recovers mostly identical parameters for the whole synthetic curve (going from 0.6 to 7T, N = 100).(b) Optimization of parameters using only a small portion of the synthetic curve (bluish area in a, N = 7) efficiently recovers the same parameters.

Figure 4 .
Figure 4. Unmixing of susceptibility curves with two ferromagnetic components.The inversion procedure was performed by first creating a forward model to be used as an initialization for the optimization step. and  are the models calculated from the inverted parameters.Model parameters are given in Table2.

Figure 5 .
Figure 5. Sensitivity tests in synthetic models.(a) Varying the contribution of the random noise and (b) the size of the sample for the  1 + 3 case.In scenarios A and B, the noise scale (   ) or the number of samples (n) is varied.  and   are the resulted models for each of these.For the  1 + 2 case, the same tests are performed (c and d), where constraining the coercivity of one of the components using a priori information will produce very similar models to the observed data.

Figure 6 .
Figure 6.Three-component case inversion.(a) The shape of the curve indicates the presence of at least three different components, which are inverted readily using the Hist-unmix script.However, for the  3 + 4 + 5 case, three (b) or two (c) components can explain the data.The log of the objective function for variable coercivities (  () and  () ) while fixing  () and the other parameters indicates a single minimum (b').However, when assuming a two-component case for the  3 + 4 + 5 curve and fixing all of the other parameters with exception of the coercivities (  () and  () ), local minimum arises (c').Nevertheless, our inversion procedure reaches an apparent global minimum in both cases (white square).

Figure 7 .
Figure 7.Testing the sensitivity of the model for mixtures of superparamagnetic fractions with higher coercivity populations.When employing a superparamagnetic (SP) fraction with the same properties as those of the SD and MD fractions (only varying   ), it becomes difficult to distinguish the SP contribution in both cases.Constraining the coercivity of one of the components to zero allows the user to test if (mathematically) a SP population could explain part of the observations.For the SP populations,   is fixed at 1.

Figure 8 .
Figure 8. Paleomagnetism and magnetic mineralogy of Sete Lagoas (BB) and Salitre (IR) formations.(a) Zijderveld diagram of a thermally demagnetized sample from the Sete Lagoas Formation, (b) the mean-site directions of C-component and (c) B-component.In (d) Lowrie-test results for a sample from the Bambuí formation.(e-h) are the equivalents for the Salitre Formation.Data are taken from D'Agrella et al. (2000) and Trindade et al. (2004).

Figure 9 .
Figure 9. Characteristic magnetic hysteresis of carbonate samples for Sete Lagoas (a and b, BB samples), and Salitre (c and d, IR samples) formations.Samples are not corrected for diamagnetic/paramagnetic contributions, since these are accounted for in our model.Boxplots (e and f) are based on the 8 samples from each formation and indicate the modeled contributions of paramagnetic  (0) and ferromagnetic  (  ) fractions, respectively for Sete Lagoas and Salitre formations.

Figure 10 .
Figure 10.Boxplot distributions of the low (C a ) and high (C b ) coercivity components of samples from the Sete Lagoas (a to c) and Salitre (d to f) formations, obtained after modeling with Hist-unmix.Diamonds are statistical outliers.

Figure 11 .
Figure 11.Examples of the inversion procedure for samples of the Sete Lagoas (a and a') and Salitre (b and b') formations, showing the lower and higher coercive components (C a and C b , respectively).The paramagnetic contribution is represented by the separation of the ferromagnetic components (blue and green lines) from the whole susceptibility spectrum.(c and d) are the M rs /M s ratios (calculated) for the C a and C b components.
that faithfully represent the ferromagnetic data and introduce random noise  () drawn from a normal distribution centered at   inv and a given standard deviation.The perturbed models are calculated through Equation 5 with a new set of perturbed parameters inv n-times.Sequentially running the inversion procedure (Section 2.2) enables optimization of    .If this procedure is repeated n-times, we can produce an average model of disturbed solutions  (   ) To avoid getting stuck in local minima, the user can create N new array of inputs