Chasing Tails: Insights From Micromagnetic Modeling for Thermomagnetic Recording in Non‐Uniform Magnetic Structures

Paleointensities are key to understanding the formation and evolution of Earth and are determined from rocks which record magnetic fields upon cooling; however, experimental protocols for estimating paleointensities frequently fail. The primary reason is that laboratory protocols assume that rocks are dominated by uniformly magnetized, single‐domain grains, instead of much more common non‐uniformly magnetized grains. Our model for larger grains shows a multiplicity of stable domain states; with preferred states changing as a function of temperature. We show that domain state distribution depends on the thermal history of the sample—in nature and the laboratory. From numerical thermomagnetic modeling, we show that particles with non‐uniform domain states will theoretically fail standard experimental paleointensity protocols, preventing us from determining reliable ancient geomagnetic field intensities. We propose that recognizing this type of behavior, and the resulting bias, will yield more reliable paleointensity records, and a better understanding of the Earth.

• Micromagnetic models of elongate 100 nm magnetite grains suggest several stable domain states that change as a function of temperature • The preferred state is also a function of cooling rate and so blocking and unblocking of remanence depends on the thermal history • Such particles will likely have unequal blocking and unblocking temperatures causing Arai plot curvature and biased paleointensity estimates

Supporting Information:
Supporting Information may be found in the online version of this article.
Failure of the reciprocity of blocking and unblocking is frequently observed as a magnetic tail in paleomagnetic observations, with little understanding of the underlying causes and consequences. Experimental evidence (Dunlop & Ozdemir, 2000, 2001 has shown that for non-uniformly magnetized magnetite particles (that is grains above the critical SD grain size), as much as 90% of magnetic remanence can unblock either below or above the blocking temperature. The presence of low temperature tails, where T u < T b , implies that more remanence can be demagnetized at any temperature T than can be replaced by remagnetizing to the same temperature and external field strength. These then cause curvature of the Arai plots (Nagata et al., 1963) from which paleointensity estimates are determined. Attempts to explain such non-ideal (non-SD) behavior have relied on phenomenological models, for example, the work of Yu and Tauxe (2005), with no theoretical underpinning, or recourse to behavior of very large magnetic particles called multidomain (MD) grains (Fabian, 2000(Fabian, , 2001(Fabian, , 2003. However, MD grains are generally too magnetically weak and unstable to significantly contribute to many paleomagnetic recordings. Of far more importance is a class of non-uniform magnetic state called pseudo-single-domain (PSD) (Dunlop et al., 1974). As their name suggests, they share many of the characteristics of SD grains, but contain non-uniform magnetizations; they are often observed in the form of a single vortex (SV) (Roberts et al., 2017;Williams & Dunlop, 1989) ( Figure 1). SV magnetic states occupy a range of grain sizes at least two orders of magnitude larger than SD grains-from ∼80 to ∼2,000 nm Ferreira 2019) and so are likely to make a significant, if not dominant, contribution to the magnetic recording in many rocks.
Here we present a complete thermal simulation of the magnetic recording process in a non-uniform domain grain of magnetite. We determine how a thermal magnetic remanence is acquired and how the rate at which a grain cools in a weak magnetic field, such as Earth's, affects this recording. We simulate the standard laboratory demagnetizing protocol in order to recover the paleointensity and show that the predicted thermal behavior explains many commonly observed features found in experimental data with profound implications as to their reliability.
The continuum theory of micromagnetism was described by W. F. J. Brown (1963) and these ideas have been exploited in three-dimensional numerical micromagnetic models (Schabes & Bertram, 1988;Williams & Dunlop, 1989) which can determine the magnetic structure, remanence, and stability of SV domain states as a function of mineralogy, grain size, shape, and temperature (Ó Conbhuí et al., 2018). A surprising result from micromagnetic simulations is that SV domain states have thermal and temporal stabilities that often exceed those of SD grains (Muxworthy et al., 2003;Nagy et al., 2017;Winklhofer et al., 1997), yielding room-temperature domain states that can retain their structure over many billions of years, far exceeding the age of the Solar System. Yet, as shown by Nagy et al. (2017), some SV particles can support domain states that are aligned along "hard" magnetic-anisotropy axes-these states are unstable over geological time. More worrying is that these particles may have more than one stable axis of magnetization leading to a multiplicity of choice for magnetic directions in violation of the premise of Néel theory and a potential cause for the failure of reciprocity and additivity required for reliable paleointensity estimates (Thellier & Thellier, 1959).

Methods
In order to characterize the ability of SV grains to reliably record the geomagnetic field, we consider a grain of the most abundant naturally occurring magnetic mineral, magnetite. We examine a grain size of 100 nm ESVD with a 30% elongation along [100]. These characteristics were chosen to be just above the critical size for transformation from SD to SV, and with a moderate elongation that avoids the atypical domain stabilities found in grains with equidimensional morphologies. Such grains have six possible orientations of the SV core: two long axis aligned 〈100〉 and four short axis aligned 〈011〉 directions. Figure 1 shows three of the possible six SV orientations, the other three being anti-parallel. The essential thermomagnetic characteristics of this grain are contained within the calculated domain state transition phase diagram shown in Figure 2. This diagram illustrates the relaxation times (which are a function of the energy barriers between states and can be translated into probabilities) between every possible transition as a function of temperature. Combining this information with knowledge of the net magnetization of each domain state as a function of temperature, we determined the magnetization of a randomly oriented assemblage of such particles for changes in thermal history and varying external magnetic field.
In this study, we examine a grain size of 100 nm ESVD with a 30% elongation along ⟨100⟩ -which we refer to as the "long-axis" or axis of elongation. In order to build a thermal model, we first identified the local energy minimum (LEM) states that the grain can support from 20.0°C to 579.0°C, in steps of 0.5°C. At each temperature, a total of 100 random initial guesses were minimized using the open-source finite element micromagnetic modeling package MERRILL (version 1.3.5) (Ó Conbhuí et al., 2018) to map out each of the possible LEM states available. The compiled code for MERRILL is available at https://www.rockmag.org and temperature dependent material parameter values for magnetite's saturation magnetization, M S (T); magneto-crystalline anisotropy, K 1 (T); and exchange, A(T), are detailed in Ó Conbhuí et al. (2018). It was found that at each temperature, the most energetically favorable state was a SV structure with a large scatter in the number of possible directions across all temperatures (see Figure S1 in Supporting Information S1). A more detailed analysis, however, reveals that the most energetically favorable domain states are SV states aligned in the ⟨100⟩ direction (the axis of elongation, Figure 1a) and vortex states aligned in the ⟨011⟩ direction ( Figure 1b). When these states are used as starting states for LEM minimization (across all temperatures), directional scatter in SV states is significantly reduced ( Figure  S2 in Supporting Information S1). The minimized end members are then used in energy barrier calculations.
Energy barriers were calculated between each long-and short-axis SV combination, resulting in long-long, long-short and short-short energy barrier paths (Berkov, 1998;Fabian & Shcherbakov, 2018) using MERRILL and relaxation times/frequencies were derived from energy barriers using the Arrhenius equation, where τ 0 is the fundamental switching time (taken to be 10 −9 s), ΔE ij is the energy barrier between domain structure i and domain structure j (in Joules), k B is Boltzmann's constant and T is the temperature in degrees Kelvin. These relaxation times are presented in Figure 2.

The Thermal Model
In addition to energy barriers, net remanence vectors for short-axis and long-axis vortex domain states were calculated as a function of temperature. Graphs of the components of this vector are shown in supplementary materials (Figures S3 and S4 in Supporting Information S1). In our thermal model, we exploit the symmetry of the possible short and long net remanence vectors in the following way. First the short-axis vortex net remanence vector is rotated by 90° about the 〈100〉 axis three times; this results in four directions denoted ρ 0 , ρ 1 , ρ 2 , and ρ 3 . Next we reflect the long-axis vortex net remanence vector in the 〈010〉∧〈001〉 plane-this results in two additional directions denoted ρ 4 and ρ 5 and the total is six orientations for the net magnetization.
In order to calculate the thermal behavior we build a Markov model based on our observations of domain state and energy barrier calculations. Possible LEM states are summarized in Figure 3. Each of the possible domain states is represented by a blue arrow (as seen in the inset of Figure 1). Transitions between short-and long-vortex states are color coded using the convention in Figure 2, with double headed arrows depicting symmetric energy barriers (this only occurs for the short-short energy transition). No direct transitions were observed to occur between anti-parallel domain structures (such transitions always go via an intermediate state) and hence there is never an arrow between short-short and long-long directions that are separated by a 180° rotation. This behavior is illustrated in Figure 3 where the transition between long-axes antiparallel states ρ 4 and ρ 5 must always go via a short axis LEM state ρ 1 or ρ 2 . Similarly antiparallel short axes transitions such as ρ 0 and ρ 2 can occur either through another short axis state such as ρ 1 or ρ 3 , or via a long axis aligned state such as ρ 4 or ρ 5 .
The Markov model in Figure 3 is then used to construct the infinitesimal generator matrix Q (Fabian & Shcherbakov, 2018) with Relaxation times for switching between domain states as a function of temperature. The short-short transition is in green, the long-short transition is in blue and the short-long transition is in red. Direct transitions between anti-parallel domain structures never occur; such transitions always go via another structure (e.g., long to short to anti-parallel long), hence long-long transitions are not depicted. The blue shaded region is the "activation zone" where the magnetization is blocked.
and ΔE ij is the energy barrier between state ρ i and ρ j , M is the net remanence vector as a function of temperature (see Figures S3 and S4 in Supporting Information S1), ΔM ⋅H is the energy contribution due to a (weak) external  Figure 2), with double headed arrows indicating symmetric energy barriers. The green red and heavy blue arrows represent the energy barriers for short-to short-vortex, short-to long-vortex, and long-to short-vortex structures calculated using MERRILL. These energy barriers change as a function of temperature. The short-vortex structures ρ 0 -ρ 3 were observed to participate in transitions between each state perpendicular to it, whereas the long-vortex structures denoted by ρ 4 and ρ 5 are only seen to participate in transitions to short-vortex structures (we did not observe a direct transition between a domain structure and its anti-parallel equivalent in our simulations).
field H (Fabian & Shcherbakov, 2018), T is the temperature in degrees Kelvin and k B is Boltzmann's constant. The diagonal elements, σ ii have the form or more simply the diagonal elements of Q are the negative sum of the off diagonal elements in any given row. The quantities ΔE ij represent the short-to long-axis vortex energy barrier whenever i = 0, 1, 2, 3 and j = 4, 5 (red arrows in Figure 3), the long-to short-axis energy barrier whenever i = 4, 5 and j = 0, 1, 2, 3 (heavy blue arrows in Figure 3), and the short-to short-axis energy barrier whenever i = 0, 1, 2, 3 and j = 0, 1, 2, 3 (green arrows in Figure 3).
The occupancy (column) vector = [ 0 1 2 3 4 5] T contains the fraction of grains, or probability of occurrence, in each of the four short and two long domain states. This vector evolves as a function of temperature (and therefore time) and may be used to calculate the net remanence for an ensemble of grains with a fixed orientation-or equivalently with an applied field in a fixed direction. Mono-dispersions of grains with different orientations can then be simulated by appropriately rotating the applied field and taking the total sum of occupancy vectors.
Finally, our model evolves as a function of temperature, which is itself a function of time. We use Newtonian cooling to model decreases in temperature, where T 0 is the initial temperature at time t 0 , T 1 is a reference temperature at t 1 and T amb is the ambient temperature. On the other hand, we model temperature increases using linear heating whereby temperature steps Δt are discretized so that the temperature never changes by more than 0.1%. By piece-wise joining linear and Newtonian temperature sections, we can construct a temperature function T fun (t) which is versatile enough to model any laboratory protocol (see Figure S5 in Supporting Information S1). We split our model up into two major phases: a TRM acquisition phase where an assemblage of grains cools from the Curie temperature to T stop (the stopping temperature) and a demagnetization phase where sample is subjected to a standard laboratory protocol (e.g., Yu et al., 2004) with pTRM checks (repeated lower temperature steps in a field) and pTRM-tail checks (repeated lower temperature steps in zero field) as shown in Figure S5 in Supporting Information S1. The main model steps for both TRM acquisition and laboratory demagnetization are summarized in Figure S6 in Supporting Information S1.

Results and Discussion
Our simulations illustrate two fundamental differences between the SD theory used to interpret paleomagnetic observations and behavior characterized by the SV domain state(s). The first is that the grain magnetization strength changes with the orientation of the SV core (light blue arrows in the inset to Figure 2). In our example, the intensity of the long-vortex (Figure 1a) is approximately 1.4 times larger than the hard-aligned short-vortex direction, (Figures 1b and 1c). This is in contrast to uniaxial SD grains, for which intensity of magnetization is invariant with choice of magnetization direction. The second difference is that this SV grain has a multiplicity of domain states whose relative stabilities change with temperature. In SD grains, although it is possible to find locally stable states aligned along different grain axes, there will be a single preferred axis of the magnetization, independent of temperature.
The changing alignment of domain orientation preference with temperature (ADOPT) leads to very different magnetic recording behavior in assemblages of grains dominated by SV as compared to SD domain states. In particular, the proportion of grains in each domain state, and thus the magnitude and stability of the room temperature magnetic remanence, will depend upon the sample's thermomagnetic history and the rate at which the sample is cooled, especially through the "activation zone," shaded blue in Figure 2. For example, when a rock cools quickly (blue curve in Figure 2), the long axis aligned grains are the first to become stable (at ∼490°C at the 1-min relaxation time line in Figure 2). By the time the sample has cooled below ∼180°C, when short-aligned SV states are the most energetically preferable, the relaxation times to realign to short axes are in excess of 10 6 s (11.6 days), rarely attained in laboratory observations for most grains. For fast cooling rates we therefore expect room-temperature magnetizations to be dominated by domain states aligned with the long axes of the grain and to have relatively high magnetizations. Conversely, following the red curve in Figure 2 for slower cooling times, the domain states that become stable at the million year relaxation time just below ∼180°C are those aligned with the short axes, which then dominate the magnetization at room temperature. Such a magnetization is weaker than its fast cooled equivalent. These differences in how the final, room-temperature, TRM is distributed among short and long domain states is at the heart of much of the non-ideal magnetic behavior that we observe in the laboratory.

SV Behavior in Simulated Experimental Paleointensity Protocols
Paleomagnetic intensities are frequently determined using variations of the experimental protocols set out by Königsberger and used by Thellier-Thellier (KTT) (Koenigsberger, 1938;Thellier & Thellier, 1959). These involve repeated thermal cycling at increasing peak temperatures, and for each temperature cycle the original (natural) thermal remanent magnetization (NRM) is replaced by one produced in the laboratory. We can simulate the acquisition of an NRM and its subsequent laboratory treatment by incorporating the relaxation phase space of Figure 2, modified by a bias field that represents the geomagnetic or laboratory field, to create a thermal activation model similar to those used in SD models Shcherbakov et al., 2021).
The fundamental issue for paleomagnetic observations is that the NRM will be carried by a variety of domain states, not all of which are in their lowest possible energy configuration. As we have seen in the previous section, high-energy states at high temperatures can be trapped by the original slow cooling of the sample from above the T C , and subsequently impossible to remagnetize into the same occupancy of domain states by (faster) laboratory heating to temperatures below T C . Since the percentage of these high-temperature, high-energy domain states is inversely proportional to cooling rate, we expect to see larger remagnetization differences as the ratio of NRM to laboratory cooling rate increases. Unless quickly cooled materials are targeted (e.g., the studies of Pick and Tauxe (1993) and Cromwell et al. (2015)), the original NRM cooling rate for the geological sample is different, most often considerably slower than, that used in laboratory remagnetizations. Therefore, we expect a variety of commonly observed non-ideal magnetic behavior to occur; these often lead to biased paleointensity estimates as seen by, for example, Tauxe et al. (2021).
The behavior predicted from our simulations is illustrated in the Arai plots (Nagata et al., 1963) shown in Figure 4 for increasing NRM cooling rates. These were produced for a mono-dispersion of grains governed by the thermal behavior outlined in Figure 2 for different cooling times.
In these examples, we use the version of KTT protocol described by Yu et al. (2004) which was designed to reveal non-reciprocal pTRM blocking and unblocking in the experiment. Arai plots show progressive demagnetizations against partial TRM (pTRM) remagnetizations for increasing temperatures T i , so that the gradient (red lines) is proportional to the ancient geomagnetic field intensity. Arai plots are a fundamental interpretation tool in KTT-type experimental paleointensity studies, where a number of phenomenological selection criteria have been devised based on the deviation of Arai plot data from a single value gradient throughout the entire range of remagnetization temperatures (Paterson et al., 2014).
The sequence of simulated Arai plots show a clear trend as a function of cooling rate which reflects two related processes. The first is that already well described in SD theory Shcherbakov et al., 2021) where more slowly cooled samples remain in thermal equilibrium to lower temperatures resulting in more strongly magnetized samples. Such cooling rate effects in SD grains typically amount to ∼5% increase in remanence per order of magnitude decrease in the cooling rate (E. M. Brown, 1984;Halgedahl et al., 1980;Leonhardt et al., 2006;Nagy et al., 2021). In assemblages of SV grains, our model demonstrates a more pronounced effect caused by the ADOPT phenomenon, whereby different cooling rates will create populations of grains dominated by different domain states. Thus NRM cooling times exceeding ∼10 6 s or longer will result in increasing numbers of grains occupying short axis aligned SV states, whilst laboratory remagnetization typically occurs from some 40 min to 10 hr, and will result in a higher proportion of long axis aligned SV states. This realignment of domain states occurs continuously at all remagnetization temperatures, but predominantly where domain state relaxation times approach the 40 min laboratory cooling times, at temperatures between ∼160° and 190°C. This cooling-rate dependency may explain why Muxworthy et al. (2014) did not observe a discrepancy for mono-dispersions of near-identical SV particles in their experimental KTT study; their synthetic NRM used the same cooling rate as the laboratory TRM.
A fundamental prerequisite for multi-step paleomagnetic intensity observations is that rock-magnetic recorders must obey the law of reciprocity, whereby a sample magnetized (blocked) by cooling from a temperature T b is completely demagnetized (unblocked) at a temperature T u such that T b = T u . Our simulations demonstrate that this will never be precisely obeyed for grains that exhibit ADOPT, since NRM demagnetization (unblocking) by zero field heating from room temperature to a temperature T i will reactivate different domain state populations compared to remagnetization (blocking) by in field fast cooling from the same temperature.
The failure of reciprocity can be recognized in experimental protocols by checking residual magnetization after demagnetization of a pTRM emplaced at temperature T i , and is often referred to as a pTRM tail (Dunlop & Ozdemir, 2000;Shcherbakov et al., 1993Shcherbakov et al., , 2000Y. J. Yu & Dunlop, 2003). Most modern experiments incorporate such a test. Our simulations provide a theoretical model of pTRM tails in assemblages of SV grains, and demonstrates that the existence of a pTRM tail is an important indicator of both reciprocity failure and the presence of a multiplicity of domain states with different blocking temperatures. The presence of tails in our model is explained by reactivation of low remanence, short-aligned SV domain states (blocked during slow cooling) that switch to long aligned SV states at temperatures around 175°C; these cannot be reached by subsequent thermal cycling to the same temperature because of their high relaxation times. This remagnetization drift to lower remanence states for low-temperature pTRMs, and higher remanence states at higher pTRM temperatures will also produce a distinct curvature of the Arai plot. Such Arai plot curvature has been suggested as an indicator of poor magnetic recorders (Paterson, 2011) and has become standard in some experimental protocols (e.g., Cromwell et al., 2015;Leonhardt et al., 2004).
pTRM tails are not produced for the grain size investigated here at temperatures above about 270°C as it becomes possible to reactive all domain states, and thermal cycling at the same temperature should yield near-identical magnetizations. Instead, we point out one of the most noticeable features of the simulated Arai plots: the presence of high-temperature "hooks" produced when the NRM has been completely demagnetized by zero-field heating to a temperature T p but laboratory remagnetization at temperatures T > T p produces a pTRM of decreasing magnitude. According to SD theory (Néel, 1988), pTRMs emplaced above T p should always produce the same magnitude of pTRM. In our model of SV themoremanence, the presence of hooks reflects the progressive reactivation of high energy domain states that then switch to lower energy states which may have a different net magnetization. For the example we illustrate here, Figure 2 shows that pTRMs blocked at increasing temperatures first reactivate the short-short transition producing no domain state change, but at slightly higher temperatures short-long activation will remagnetize grains along the long axis which have larger magnetizations. At still higher temperatures, the high-remanence, long-axis aligned states are activated and realigned to the lower remanence short axis states. It is this last step that results in decreasing pTRM magnitudes with increasing pTRM emplacement temperatures beyond ∼300°C.
The results described here are for a mono-dispersion of one particular grain size and shape, chosen so that it is possible to nucleate a variety of different vortex states with different magnetizations. For idealized grain geometries, the grain-size range where SV states share multiple choices of axes for domain stability (MCADS) is expected to be narrow, and primarily near the transition between SD and SV states (Nagy, Williams, Tauxe, Muxworthy, & Ferreira, 2019;Nagy et al., 2017;Valdez-Grijalva et al., 2018). In equidimensional cuboctahedral magnetite, this size range is only ∼20 nm centered on a grain size diameter of 90 nm. The mean size and range of the MCADS zone will increase with grain elongation, and will vary with mineralogy, and likely be more prevalent in natural grains with irregular morphologies and multiple grain axes. The location of the "activation zone" will also be sensitive to grain size and move to higher or lower temperatures dependent on grain shape and geometry. Experimental evidence for high energy and high remanence MCADS grains is provided by the high ratio of thermal to anhysteretic remanence which exhibits a peak at exactly the predicted MCADS grain sizes (Paterson et al., 2016). Not all paleomagnetic samples will be affected by the blocking temperature reciprocity failure described here, but fine-grained rock samples, say in lava flows cooled over days to weeks and traditionally thought to contain ideal magnetic recorders, might ironically contain high percentages of MCADS grains, and so contribute to the high failure rate reported in experimental paleointensity studies. Our models indicate that pTRM tails and Arai curvature are good indicators of reciprocity failure and the presence of MCADS grains, and so provide a theoretical underpinning of the inclusion of curvature in the sample selection criteria.
The SV model of PSD thermomagnetic recording we have put forward provides some insight into why experimental determination of the paleomagnetic field strength has proved so difficult, but reassuringly indicates that careful experimental protocols can correctly identify poor recorders. The predictions account for commonly encountered non-ideal behavior in paleomagnetic intensity experiments such as pTRM tails and Arai plot curvatures. The magnetic domain state changes predicted here and characterized in our relaxation phase maps, implicate MCADS grains as the source of "fragile curvature" Yamamoto et al., 2022) where a long-period laboratory storage of a sample re-configures the relative room-temperature population of different domain states, and thus changes their behavior in subsequent thermal-cycling experiments. Our results also bring additional opportunities to developments in single-grain (Tarduno et al., 2001(Tarduno et al., , 2015 and micro-paleomagnetic studies (Fu et al., 2020;Out et al., 2022). These techniques attempt to extract paleomagnetic signals from samples of a few microns based on a detailed knowledge of the magnetic mineralogy. It is now possible to identify and eliminate MCADS grains from observations and if such a procedure can be developed into an actual paleointensity experiment, it will provide more reliable estimates of the ancient geomagnetic field.

Data Availability Statement
All results reported here were generated using the open source micromagnetic modeling code of (Ó Conbhuí et al., 2018). A complete guide to installation and use of MERRILL is described here: https://blogs.ed.ac.uk/ rockmag. The data required to reproduce our results (MERRILL script and geometries) is available at https://doi. org/10.5281/zenodo.7307935.