Estimating the Effect of Cooling Rate on the Acquisition of Magnetic Remanence

The effect of cooling rate on the magnetization of rocks must be accounted for when estimating ancient magnetic field strengths. Calculating this effect is not trivial, even for uniformly magnetized grains. Here, we present an open‐source package to compute the behavior of uniaxial single‐domain grains for different temperature and magnetic fields. We revisit the problem of thermal remanence acquisition as a function of cooling rate and find that our predictions are broadly in agreement with those of Halgedahl et al. (1980, https://doi.org/10.1029/jb085ib07p03690) but differ significantly from those of Dodson and McClelland‐Brown (1980, https://doi.org/10.1029/JB085iB05p02625). We also find that remanence acquisition curves correspond well with recent experimental observations. Cooling rate corrections made using our model are at the upper limit suggested by Halgedahl et al. (1980, https://doi.org/10.1029/jb085ib07p03690) but can reduce slightly for larger (single‐domain) grains, very slow cooling rates of the original thermal remanence and large field strengths.

2 of 10 of SD grains. Using Néel's theory (Néel, 1949), Dodson and McClelland-Brown (1980) examined the effect of slow cooling on the blocking temperature of ensembles of SD particles. A concurrent effort was undertaken by Halgedahl et al. (1980), who modeled the effect of cooling rate on the acquisition of paleointensity in several different cooling scenarios (regimes). Unfortunately, there is a mismatch in predicted thermal remanent magnetization (TRM) between these twin efforts (Santos & Tauxe, 2019).
In this study, we revisit the single-domain model of remanence acquisition from Néel's theory of elongate single-domain grains, referred to as Stoner-Wohlfarth grains after Stoner and Wohlfarth (1948). We take advantage of advances in numerical computation capability since the early 1980s and provide a fast and publicly available code, written in C++, that calculates TRM gained as a function of cooling in an external field. This code uses the Boost multiprecision library (Boost, 2021) to avoid possible numerical issues that arise when calculating the fractional alignments of noninteracting grains that make up our model. We then examine a number of cooling and field regimes and produce a new set of cooling rate correction curves and find that our results agree well with the cooling rate curves provided by Halgedahl et al. (1980) for the majority of the size elongation and field scenarios in this study.

The Stoner-Wohlfarth Model
The Stoner-Wohlfarth model (Stoner & Wohlfarth, 1948) describes the energy barriers that a simple uniformly magnetized uniaxial ferromagnetic grain, with ellipsoidal geometry, must overcome to change its magnetic state in the presence of an externally applied field ⃗ . The external field makes an angle ϕ with the grain's axis of elongation ̂ as shown in Figure 1. The magnetic state is the angle θ between the unit magnetization vector ̂ and the grain axis ̂ . An expression for the magnetic energy of the system envisioned originally by Néel (1949) is given by Dunlop and Özdemir (2001, p. 207) and can be written as with the temperature dependent constants C 1 (T), C 2 (T), and C 3 (T) given by Figure 1. The Stoner-Wohlfarth model of magnetization assumes an ellipsoidal grain oriented along the vector ̂ . The applied field ⃗ (of strength H) and the magnetization ̂ makes angles ϕ and θ (respectively) with ̂ . The model assumes that the magnetization will rotate within the ̂ − ⃗ plane, and so for an arbitrary angle θ we can always recover a threedimensional magnetization vector.
NAGY ET AL.
10.1029/2021GL095284 3 of 10 M s (T) is the saturation magnetization at temperature T, μ 0 is the permeability of free space and the particle volume is v. The strength of an externally applied field is H and its direction is given by ϕ, as described previously, with θ the direction of magnetization.
The demagnetizing factors of a prolate ellipsoid, with aspect ratio m, are defined in Cullity and Graham (2011, p. 54), with N a and N b corresponding to the demagnetizing factors along the long and short axes, respectively, as shown in Equations 5 and 6 It should be noted that in this study we quote elongation as a percentage as opposed to aspect ratio where elongation is defined by m = (elongation + 100)/100. This means that an aspect ratio of 1.3 corresponds to an elongation of 30%.
In order to find the critical points for the energy of a Stoner-Wohlfarth particle, we look for the values of θ where the first derivative of Equation 1 with respect to θ is zero, doing this gives ( , ) = 2 1( )cos sin + 2( )sin − 3( )cos = 0.
There is no general analytical solution for Equation 7 except for the special cases when ϕ = 0 and ϕ = π. However, we can numerically find the zeros by making the substitution θ = i log(x) where = √ −1 . This then transforms Equation 7 from a trigonometric one into the polynomial from which we can form the upper Hessenberg matrix (Press et al., 2007, p. 469) The eigenvalues of  are the zeros of the polynomial version (Equation 8), and are found using the Eigen linear algebra library (Guennebaud & Jacob, 2010). Then we calculate the critical θ-values, denoted θ k , by using our original transform = log ( ) . This results in θ k ∈ [ −π, π], with each θ k solving Equation 7. Any θ k values that have nonzero imaginary parts are discarded as these do not represent real magnetization directions.

Thermal Theory of Remanence
We briefly review the thermal theory of single-domain remanence with particular reference to the implementation details in our C++ code. We are interested in both a "cooled remanence" which solves the thermal equations with the assumption that grain assemblages spend only a finite amount of time at a given 4 of 10 temperature and the "equilibrium remanence" which is the theoretical limit for which a grain has spent an infinite amount of time at each temperature step.

Cooled Remanence
Once the critical values of Equation 1 are found, we can evaluate whether they correspond to energy minimum/maximum states by taking the second derivative of Equation 7, which results in When Equation 10 is positive for any critical value θ k , then we have found a local energy minimum (LEM) state and the critical value is denoted θ k,min . Likewise θ k values that make Equation 10 negative correspond to local energy maxima and are denoted θ k,max . The energy barrier for a two-state system is then given by We take the energy barrier as the transition energy between any two LEM states, θ k,min and θ j,min , as this represents the physical path that the magnetization would take when transitioning between any two LEM states. The isothermal transition rate matrix (Fabian & Shcherbakov, 2018), denoted P, may now be formed from the above energy barrier calculations by assuming that a grain population (given by the vector ⃖ ⃗ described below), has experienced the same field and temperature conditions for a given time Δt where T is the temperature of the grain in Kelvin, 1/τ 0 = 10 10 Hz is the attempt frequency (Dunlop & Özdemir, 2001), k B is Boltzmann's constant and "exp" is the matrix exponential function (see Appendix A1). Equation 12 may then be used to calculate an updated grain population ⃖ ⃗ +Δ according to For a monodispersion of grains, which is a population of grains with a single size and shape shape, we define a "population vector." Each element of the population vector represents the fraction of grains that occupy a particular magnetization state. This means that ρ k,t+Δt must sum to unity and that each element indexed by a specific LEM state k must be consistent with its predecessor ρ k,t−Δt . The normalized magnetization is then given by where ⃖⃖ ⃗( min) represents the conversion of the magnetization LEM angle, that solves the Stoner-Wohlfarth equations described above, back to a three-dimensional vector (see Appendix A2).

Equilibrium Remanence
To estimate the effect of cooling rate, we need to also estimate the equilibrium TRM, which is defined as the magnetization reached when an ensemble (population) of particles have experienced a given field and temperature for an infinite amount of time. The equilibrium population vector ⃖ ⃗eq components are given by Dunlop and Özdemir (2001, pp. 213) as 5 of 10

Cooling Models
The effect of cooling was calculated for a number of different cooling regimes with temperatures given by classical Newtonian cooling Here, ( 0, 0) is an initial time-temperature pair which we take to be t 0 = 0 s and T 0 is the Curie temperature of magnetite in degrees centigrade. The other known time-temperature point along the cooling curve, ( 1, 1) , is taken to be T 1 = 15.15 °C (since in the Newtonian cooling model the ambient temperature is an asymptote) and for t 1 -the time taken to reach T 1 -we use t 1 = 6 × 10 e s (with e ∈ 1, 2, …, 15) to give a range of "rapid" to "slow" cooling rates. Finally, the ambient temperature, T amb , is 15 °C.

Results and Discussion
The results in Figures 2 and 3 show that the magnitude of the TRM as a function of applied field in "rapid" and "slow" cooling regimes. The main difference between Figures 2 and 3 is that the applied field for Figure 2 is directed along the grain axis ̂ = ⟨1, 0, 0⟩ whereas the field in Figure 3 is directed along ⟨1, 1, 1⟩ , forming an angle of 54.7° with respect to the grain axis.
In both figures, it can be observed that TRM increases as a function of grain size, expressed as equivalent spherical volume diameter (ESVD), and remains approximately linear as a function of applied field. In all cases, the TRM acquired increases from rapid to slow cooling times as is evident from the way the solid lines fan out from left to right as the cooling times become longer. We expect this is because for slow cooling, the grain has more opportunity to equilibrate with the external field, allowing a stronger magnetization to be acquired. It may also be observed that TRM drops (the solid lines fan out less) as the particles become more elongate. In order to explain this effect, it should be noted that, upon cooling, the earlier a TRM acquisition curve departs from its equilibrium behavior, that is, its blocking temperature, the smaller its room temperature remanent magnetization will be. For highly elongate grains, the rapid increase of energy barriers upon cooling results in a higher blocking temperature and so lower TRM as can be observed in Figures 2 and 3.
Energy barriers to domain switching in Stoner-Wohlfarth particles for fields parallel to the grain axis are in general higher than at other angles. For small fields similar to the Earth's field, the grain's magnetization will always lie along its elongation axis and so the difference between the two possible states is higher in the field-parallel case as opposed to some other angle. This means that TRM acquisition is more efficient when the field is applied parallel to the elongation axis, as in Figure 2 when compared to the case when the field is applied at an angle ( Figure 3). In addition to TRM efficiency being a function of cooling rate, the curvature of equilibrium (dashed) lines is also greater for grains with the field directed parallel to the grain axis ( Figure 2) as opposed to grains with the external field directed at an angle to the grain axis ( Figure 3). For example, the 80 nm dashed line in Figure 2 reaches its saturation value at ∼100 μT, whereas the same line in Figure 3 reaches the saturation value at the much higher field of ∼175 μT. Figure 4 shows the results of our modeling along with the predictions of Halgedahl et al. (1980) (dashed black line) and Dodson and McClelland-Brown (1980) (dotted black line) along with experimental data from Santos and Tauxe (2019) and other authors (detailed in the caption of Figure 4). Our numerical calculations are for a collection of grains with no fabric, which is a monodispersion of grains over a random distribution of directions (with respect to applied field). The majority of the grain size and elongations correspond well with the predictions of Halgedahl et al. (1980). The most noticeable exception being the 30 nm 30% elongate grain (light blue line) which is border-line superparamagnetic since its volume and elongation are relatively small. Figure 5 shows TRM acquisition curves for the complete time range for a study that goes well beyond that seen in Figure 4 with an assumed laboratory cooling time of 1,000 s to a maximum cooling time of 190 Ma. A population of grains with a strong fabric (a monodispersion of grains that are all aligned with the applied field) are shown along with a set of predictions for high field strength of 210 μT. We see in Figures 5a and 5b 6 of 10 that the spread of TRM acquisition for slow cooling is relatively small in weak fields. This is not the case for stronger fields shown in Figures 5c and 5d where there is a much greater spread. This illustrates that, in weak fields at least, elongation and grain size have little effect on TRM acquisition. It should also be noted that the highly elongate grains (red) show only relatively minor variations in all field regimes. TRM acquisition is affected by grain size, though much less so in elongate grains. This is most clearly seen in the strong field regimes in Figures 5c and 5d with both the parallel field and intermediate field showing that for slow cooling, the TRM recorded decreases as a function of grain size (we see the darker lines taking on shallower gradients). It may also be observed that in the larger grains under strong field conditions, there is a slight curvature. This is again most apparent in the 30% elongate grains, indeed the 30 nm 30% elongate grain (lightest blue) plateaus for slow cooling. As observed previously, this grain size is just on the cusp of being superparamagnetic and at a particular cooling rate the "cooled" TRM acquisition curve achieves equilibrium. The threshold for superparamagnetic behavior is when the magnetization reaches equilibrium with the external field over the time span of observation. In the case of the 30 nm 30% elongate grain, the relaxation time is short enough when cooled slowly, for its thermal-magnetic behavior to achieve equilibrium, meeting the definition of superparamagnetism. In principle all cooling rate curves should eventually plateau, if the cooling rate is slow enough (see Figures 1 and 2 in Dodson & McClelland-Brown, 1980). A final observation is that grains with strong fabric and no fabric show small differences. These differences are a drop in the ratio of TRM gained since the gradient of each line becomes very slightly shallower from strong fabric to no Figure 2. Thermal remanent magnetization (TRM) acquisition versus applied field for cooling from the Curie temperature (580 °C) to 15.15 °C as a function of grain size (equivalent spherical volume diameter (ESVD)) and elongation for "rapid" (t 1 = 6 × 10 3 s) and "slow" (t 1 = 6 × 10 15 s) cooling regimes. TRM has a value of 1.0 when all particles are aligned with the field direction. Field strengths range from 30 to 210 μT and are aligned parallel to the grain elongation axis . Solid lines show TRM acquired through cooling, whereas the dashed lines show equilibrium TRM (infinite cooling time). 7 of 10 fabric; and an increase in the TRM ratio gained when the grain hits its equilibrium behavior (lightest blue line). The shallower gradient is due to the fact that in simulated monodispersions with no fabric, there are many grains that have smaller energy barriers since the field is at an angle to the axis of elongation. The light blue line plateaus later (with higher ratio of TRM) for the same reason and so the cooling effect is reduced. This effect can also be seen by comparing the plateau between the weak field and the strong field in samples between the same fabric (i.e., between Figures 5a and 5b and Figures 5c and 5d) since in stronger fields grains hit their equilibrium behavior sooner.

Conclusions
In this study, we have presented a model for calculating the TRM acquisition as a function of field and cooling rate and have found good correspondence with experimental data for single-domain grains. Of the previous published models we find that our results are very close to the predictions of Halgedahl et al. (1980). Our model also demonstrates subtle variations in recording as a function of grain size and shape; however, we also show that there is relatively little variation in remanence acquisition as a function of field strength and direction (at least for weak fields like the Earth's).
The source code for the model that we have presented is freely available at https://github.com/Lesleis-Nagy/ sd-cooling (version 1.0.1 was used in this study). Currently, it based on simple Stoner-Wohlfarth modeling described; however, the thermal theory of remanence described in this study is also applicable to grains with much more complicated magnetizations and switching mechanisms. For more realistic grains, we with  (̂min) the the 3 × 3 rotation matrix given by  Dodson and McClelland-Brown (1980) (dotted black line) for the complete time range in this study with (a) an assemblage in weak field with strong fabric, (b) an assemblage in a strong field with strong fabric, (c) an assemblage in weak field with no fabric, and (d) an assemblage in strong field with no fabric. The color scheme is the same as in Figure 4. Cooling rates are calculated with respect to a laboratory reference cooling time of 1,000 s. In order to apply a cooling rate correction, simply divide the sample age (in seconds) by the laboratory reference time and take the base 10 logarithm, after that the ratio of remaining TRM can be read off from the graph (depending on the grain characterization of the sample).
 ̂ , , = (1 − cos ) + sin (A10) for the off-diagonal components. For the special case where ⃖⃖⃖ ⃗ and are parallel, we assume that the magnetization is also parallel with ̂ .

Data Availability Statement
Raw data used in this study are archived with https://doi.org/10.5281/zenodo.5086636. The specific version of the code used to produce the data in this study is archived with https://doi.org/10.5281/zenodo.5085691 and updated at https://www.github.com/Lesleis-Nagy/sd-cooling.