Size Ranges of Magnetic Domain States in Tetrataenite

Paleomagnetic studies of meteorites provide unique constraints on the evolution of magnetic fields in the early solar system. These studies rely on the identification of magnetic minerals that can retain stable magnetizations over ≳4.5 billion years (Ga). The ferromagnetic mineral tetrataenite (γ''‐Fe0.5Ni0.5) is found in iron, stony‐iron and chondrite meteorite groups. Nanoscale intergrowths of tetrataenite have been shown to carry records of paleomagnetic fields, although the effect of magnetostatic interactions on their magnetic remanence acquisition remains to be fully understood. Tetrataenite can also occur as isolated, non‐interacting, nanoscale grains in many meteorite groups, although the paleomagnetic potential of these grains is particularly poorly understood. Here, we aim to improve our understanding of tetrataenite magnetization to refine our knowledge of existing paleomagnetic analyses and broaden the spectrum of meteorite groups that can be used for future paleomagnetic studies. We present the results of analytical calculations and micromagnetic modeling of isolated tetrataenite grains with various geometries. We find that tetrataenite forms a stable single domain state at grain lengths between 6 and ∼160 nm dependent on its elongation. It also possesses a magnetization resistant to viscous remagnetization over the lifetime of the solar system at 293 K. At larger grain sizes, tetrataenite's lowest energy state is a lamellar two‐domain state, stable at Ga‐scale timescales. Unlike many other magnetic minerals, tetrataenite does not form a single‐vortex domain state due to its large uniaxial anisotropy. Our results show that single domain and two‐domain tetrataenite grains carry an extremely stable magnetization and therefore are promising for paleomagnetic studies.

The ability to conduct reliable magnetic studies of early solar system materials is contingent on the presence of magnetic carriers capable of retaining a magnetic remanence and preserving that remanence over >4.5 billion years (Ga). These recorders tend to adopt a single domain (SD) (Néel, 1955) or single vortex (SV) (Einsle et al., 2016;Shah et al., 2018) magnetic domain state. SD and SV states are typically the most resistant to remagnetization because of their higher coercivities, leading to magnetic relaxation times greater than the age of the solar system. In meteorites and lunar samples, the ferromagnetic minerals kamacite (α-Fe 1−x Ni x for x ≲ 0.05) (Garrick-Bethell & Weiss, 2010;Gattacceca et al., 2014;Weiss et al., 2010), pyrrhotite (Fe 1−x S for x < 0.17) (Weiss, Fong, et al., 2008), and magnetite (Fe 3 O 4 ) (Fu et al., 2014;Gattacceca et al., 2016) are found to hold paleomagnetic records. According to analytical calculations and micromagnetic modeling, kamacite and magnetite grains are in the SD state from ∼15 to ∼300 nm (Butler & Banerjee, 1975a;Muxworthy & Williams, 2015;Shah et al., 2018), and ∼30-∼1,000 nm (Butler & Banerjee, 1975b;Muxworthy & Williams, 2006;Muxworthy et al., 2003;Nagy et al., 2017), respectively, depending on the elongation of the grain. Above that size range, magnetite and kamacite occupy an SV or multidomain (MD) state. The SP to SD transition of pyrrhotite has been calculated to be 46 nm (Dunlop, 2021), but the SD to MD transition is not well characterized. Bitter pattern observations of pyrrhotite's domain states place its SD to MD transition at 1.6 μm (Soffel, 1976) and hysteresis properties suggest the mineral is in the single domain for grains <3 μm (Clark, 1984). Pyrrhotite in the SV state has not been reported in natural samples, but the formation of vortices has been invoked as a potential explanation for the decrease in saturation remanence exhibited by pyrrhotite upon heating (Dunlop, 2021). Analytical calculations have also been performed for cubic greigite (Fe 3 S 4 ), which found the SD state to occur for grains 60-250 nm (Ricci & Kirschvink, 1992), while micromagnetic modeling of non-interacting cubic greigite grains places the SD range between 46 and 107 nm (Muxworthy et al., 2013). More recent micromagnetic modeling of octahedral and cuboctahedral greigite indicates the SD to SV transition occurs at 50-56 nm (Valdez-Grijalva et al., 2018).
Little work has been undertaken to determine the domain states and grain size range for which tetrataenite magnetization is stable against external fields or viscous remagnetization at ambient temperatures. Einsle et al. (2018) studied the magnetic remanence acquisition of CZ islands by simulating the ordering process of tetrataenite using micromagnetic modeling. The authors found that SV taenite islands with lengths 82-92 nm would transition to SD tetrataenite via an intermediate two-domain state. The two-domain state becomes an SD state due to the annihilation of the domain wall through the presence of an external field and/or island-island magnetostatic interactions. Once such tetrataenite grains are placed in the SD state, a field >1 T is needed to re-nucleate the domain wall. Einsle et al. (2018) focused primarily on remanence acquisition mechanisms of the cloudy zone and the effect of magnetostatic interactions. However, the magnetic behavior of isolated, non-interacting tetrataenite has yet to be addressed and would provide an apt comparison to the behavior of closely packed, interacting tetrataenite grains. Additionally, non-interacting tetrataenite grains are found in other meteoritic microstructures such as plessite (Figure 1b), a microscale intergrowth of kamacite and taenite or tetrataenite, with microcoercivities

Analytical Calculation and Micromagnetic Simulation Methodology
Here we describe the methodology for our analytical calculations and micromagnetic modeling of the domain states in tetrataenite. We start by discussing the analytical formulas based on Néel theory to predict the superparamagnetic (SP) (relaxation time <4.5 Ga) to SD transition and SD to MD transition. Analytical relationships provide an intuitive manifestation of the principles of energy conservation. However, Néel theory is unable to predict the existence of an SV state or other metastable states. Therefore, we also conduct micromagnetic calculations to assess the possible existence of an SV state and the critical size range over which it might form. Together, micromagnetic modeling and analytical calculations provide a holistic view of the sizes and shapes of individual tetrataenite grains that are magnetically stable and therefore useful for paleomagnetism.

Analytical Calculations
To determine the SP to SD threshold size for tetrataenite at various axial ratios, we followed Equation 6 in Evans and Mcelhinny (1969). This is a modification of the Néel-Arrhenius equation, which defines the relaxation time as the ratio of magnetic anisotropy energy, to thermal energy, : where is the relaxation time, C is the frequency factor ∼10 9 Hz, v = 3 2 is the grain volume where is the length of the long axis, H c is the microscopic coercivity, 0 is the permeability of free space, M s is the saturation magnetization (1.39 × 10 6 A m −1 for tetrataenite; Néel et al., 1964), k B is Boltzmann's constant, and Tetrataenite islands decrease in size as distance from the clear tetrataenite rim increases due to increasingly lower Ni concentrations. Image taken on a Merlin Zeiss FEG-SEM at the MIT Material Research Laboratory. (b) X-ray photoelectron emission microscopy (XPEEM) image of duplex plessite, an intergrowth of taenite/tetrataenite and kamacite formed from decomposition of martensite with precipitates >200 nm (Goldstein & Michael, 2006;Zhang et al., 1993), in the Bacubirito ungrouped iron meteorite taken at the Advanced Light Source of the Lawrence Berkeley National Laboratory. Gray scale is indicative of nickel content. The individual tetrataenite grains stand out from the kamacite matrix due to their higher concentrations of Ni. (c) Backscattered electron microscopy image of a pyroxene grain with metal blebs in the Acapulco primitive achondrite. The tetrataenite appears brighter compared to the kamacite due to its higher Ni content. Image taken on a JEOL JXA-8200 Superprobe electron microprobe Department of Earth, Atmospheric, and Planetary Science at MIT.
T is temperature in K. In tetrataenite, magnetocrystalline anisotropy dominates over shape anisotropy. This is because the microcoercivity due to magnetocrystalline anisotropy [ = 2K∕ 0 where K = 1.37 × 10 6 J m −3 (Einsle et al., 2018) is the magnetocrystalline anisotropy constant] is larger than the microcoercivity due to shape anisotropy [ = ∆ where ∆ is the difference in demagnetization factors along the short and long grain axes (Equation S3 in Supporting Information S1)] at all axial ratios. In our calculations, we choose to consider the effects from both magnetocrystalline anisotropy and shape anisotropy, because = 1.57 × 10 6 A m −1 and = 1.39 × 10 6 ∆ A m −1 where 0 ≤ ∆ ≤ 0.5 remain within an order of magnitude. We define the total microscopic coercivity H T = + and substitute this for H c in Equation 1. This definition of H T assumes that the shape and magnetocrystalline easy axes are aligned. This may not be the case for all scenarios, but a sum of the two anisotropies places a lower bound on the SP-SD transition. Rearranging Equation 1, we find: We define the critical relaxation time for the SP-SD transition to be = 4.5 Ga since we are concerned with magnetic records from meteorites that are stable over the lifetime of the solar system. Conventionally, the SP-SD transition is defined as stable magnetization on the laboratory time scale (100 s; Tauxe, 2010), though since the grain length is dependent on ln( ) 1∕3 , the SP-SD transition at 4.5 Ga and 100 s is only a factor of ∼2 different.
We determined the SD-MD threshold size following Dunlop and Özdemir (1997), who define the threshold size to be the grain length at which the energies of an SD state and a two-domain state are equivalent. Above the critical SD size, it is energetically favorable for the grain to nucleate a domain wall rather than retaining a state of uniform magnetization with high magnetostatic energy (Dunlop & Özdemir, 1997). Following Equation 5.28 in Dunlop and Özdemir (1997) and substituting for the domain wall energy [4( ex) 1∕2 (Lilley, 1950), where A ex = 1.13 × 10 −11 J m −1 (Einsle et al., 2018) is the exchange energy], we find that the critical length of the long axis of a cuboid can be expressed as: where N is the demagnetization factor along the elongated axis. The demagnetization factor is dependent on the axial ratio, ranging from 1/3 when the grain is equant to zero when the grain is an infinitely elongated rectangular prism (Aharoni, 1998). The demagnetization factor is considered only along the elongated axis as it is the minimum energy direction of magnetization. Demagnetization factors for Equations 2 and 3 are calculated using Equation 5 from Aharoni (1998) for a cuboid with the two minor axes being the same length (see Supporting Information S1).
Assuming the lamellar shape of the domains persists to larger grains sizes, the analytical solution for the transition from a two-domain grain to a three-domain grain, or for any higher domain state, can be obtained from the extension of Equation 5.27 in Dunlop and Özdemir (1997) which defines the transition from an SD state to a two-domain state. By analogy, the transition from two domains to three domains occurs when the energy of the two-and three-domain states are the same, and more generally when the energy of an n-domain grain is the same as an n + 1-domain grain. Using the approximation in Dunlop and Özdemir (1997) Equation 5.3 that the demagnetizing energy of an n domain grain is 1/n times the demagnetization energy of an SD grain, the transition from an n domain grain to an n + 1-domain grain occurs when where w and u are the width and height of the grain, respectively, and are each equal to if the grain is a rectangular prism. Substituting v = 3 2 and solving for , we find: The SV state is not predicted by Néel theory, and while attempts have been made to analytically calculate the transition from single domain to a circular spin state (for example, Butler & Banerjee, 1975b), the circular spin state modeled has no vortex core that can account for the remanence of the grain. However, the SV threshold can be determined using micromagnetic modeling.

Micromagnetic Modeling
We conducted micromagnetic modeling of individual tetrataenite grains using the Micromagnetic Earth Related Robust Interpreted Language Laboratory (MERRILL version 1.6.4) software (Conbhuí et al., 2018). MERRILL is an open-source software package that uses finite-element modeling to simulate the three-dimensional magnetic structure of a tetrahedrally-meshed grain. MERRILL finds the magnetization state by solving for the orientation of an array of dipoles that has the minimum magnetic energy associated with exchange, anisotropy, Zeeman, and internal demagnetization fields. A full description of MERRILL can be found in Conbhuí et al. (2018).
Cuboid grain meshes were created using the meshing software MeshRRILL (Conbhuí et al., 2018) with a tetrahedral element size of 3 nm, which is the magnetic exchange length of tetrataenite as defined by Equation 12 in Conbhuí et al. (2018). The element size was limited to the exchange length to ensure that spatial variations of the magnetization within the particle were correctly accounted for in the calculation of the total free magnetic energy. Magnetic parameters for tetrataenite were taken from Einsle et al. (2018) and Néel et al. (1964), with K, A ex , and M s defined in Section 2.1. All simulations were calculated at room temperature T = 293 K. While meteorites likely experienced temperatures ∼150-280 K while in the asteroid belt and transfer to Earth , grains stable at 293 K will be almost always be stable at the lower temperatures barring phase transformations since there is less thermal energy. We are not aware of any reports of tetrataenite crystallographic phase transitions at such low temperatures.
To determine the SP-SD threshold size, each element of the mesh was initialized with a randomly oriented magnetic moment of equal magnitude. Simulations were conducted with the magnetic easy axis oriented parallel to the elongated [001] axis to be consistent with our analytical calculations. Once initialized, MERRILL minimized the total energy of the system using a modified conjugate gradient method (Fabian & Shcherbakov, 2018). For each grain size and axial ratio modeled, 100 least-energy magnetization (LEM) solutions were completed, each with a different randomized initial state, to calculate the least-energy state of that geometry (i.e., to identify the global minimum energy state). To calculate the energy barrier for viscous remagnetization, we considered two different end-member magnetization states with the lowest energy but opposite moment directions. These two configurations were passed into MERRILL, which used a nudged elastic band (NEB) method to determine the minimum-action path between the end-member states (Fabian & Shcherbakov, 2018). Our calculations considered 100 magnetization steps between the two states. The energy barrier to remagnetization was calculated by the energy difference between the saddle point on the minimum-energy path and the initial state (∆ ). The relaxation time was then determined using the Néel-Arrhenius equation (Equation 1) with T = 293 K. As discussed above, grains stable at this temperature will be stable at any colder temperatures the tetrataenite experienced. The critical relaxation time was set at 4.5 Ga to ensure that any magnetization would be retained over the lifetime of the solar system.
We determined the threshold size for a single-to-nonuniform domain (e.g., single vortex or MD) state by identifying the grain size at which the least-energy state from 100 random initial guesses transformed from a single-domain state to a nonuniform state. The resistance of these nonuniform states to viscous remagnetization was determined by calculating NEB paths and relaxation times in a similar manner to the SP-SD threshold.
Simulated hysteresis loops for grains at various axial ratios with sizes above the SP-SD threshold, including nonuniform states, were created using MERRILL to determine the behavior of the states in response to an external field. The least energy state identified for each grain geometry was used as the initial configuration. The energy of the system was then minimized in the presence of subsequent increasing external fields along the easy axis direction. The resulting least-energy state after each minimization was the initial configuration for the following external field step. External fields were increased in steps of 10 mT from 0 mT to the saturating field (depending on the geometry of the grain). The external field was then decreased in steps of 10 mT to the saturating field value in the opposite direction and increased again to reach the positive saturating field.
All simulations were conducted on the MIT Engaging Cluster using Centos 7 cores with 128 GHz RAM. The SP-SD threshold size was determined for axial ratios between 0.1 and 1 in steps of 0.1. The memory demands of MERRILL on our computation facilities limited our calculations of the threshold size for the single-to-nonuniform domain state to axial ratios between 0.55 and 1 in steps of 0.05. The largest grain modeled was a 120 nm cube with a volume 1,728,000 nm 3 comprising 269,855 elements.

Analytical and Micromagnetic Modeling of Tetrataenite Domain States
The domain state behavior of tetrataenite as a function of axial ratio and length is shown in Figure 2. Both the analytical calculations of the SP-SD transition and micromagnetic modeling show a monotonic increase of the threshold size with decreasing axial ratio. This trend follows a dependence of A −2/3 as shown in Equation 2. Because of the dominant magnetocrystalline anisotropy, the threshold volume of the grain is nearly constant with axial ratio-in contrast to the SP-SD threshold for kamacite and magnetite. Magnetite's high saturation magnetization compared to its magnetocrystalline anisotropy constant (4.80 × 10 5 A m −1 and −1.35 × 10 4 J m −3 , respectively; Dunlop & Özdemir, 1997) leads to this mineral being dominated by shape anisotropy at nearly all axial ratios. Thus, the SP-SD threshold size subsequently decreases with increasing elongation to A = ∼0.6 before growing again with further elongation (Butler & Banerjee, 1975a). A similar trend in the SP-SD threshold was found for kamacite (Butler & Banerjee, 1975b), which has a high saturation magnetization compared to its magnetocrystalline anisotropy constant (1.715 × 10 6 A m −1 and 4.8 × 10 4 J m −3 , respectively; Dunlop & Özdemir, 1997). The transition from SP to SD for tetrataenite results in a sharp increase in relaxation time. For a cube (SP-SD threshold of 6 nm), the relaxation time jumps from ∼2,560 years at 5 nm to 3.34 × 10 8 Ga at 6 nm due to the increase in the energy barrier when the SD state is reached (Movie S1).
Our micromagnetic simulations show that unlike magnetite and kamacite, which transition through a SV state prior to a MD state as their grain size increases (Muxworthy & Williams, 2006Shah et al., 2018), tetrataenite does not form a SV state but instead transitions directly from a SD state to a two-domain state ( Figure 3). This is due to the substantially higher magnetocrystalline anisotropy constant of tetrataenite compared to kamacite and magnetite (28 and 100 times larger, respectively). The large magnetocrystalline anisotropy prohibits the individual dipole moments associated with each element node from being canted relative to the magnetic easy axis. This prevents the formation of a vortex core, which would be favored if exchange interactions outweighed the anisotropy, forcing the least energy state instead to be two domains separated by a nucleated domain wall. The absence of a SV state in tetrataenite agrees with the results from Einsle et al. (2018), who found that tetrataenite grains formed a two-domain state prior to the destruction of the domain wall due to island-island interaction fields and/or external fields.
Unlike for the SP-SD threshold, we find a mismatch between analytical and micromagnetic modeling solutions for the SD to MD threshold ( Figure 2). The analytical solution for a cube (79 nm) initially overestimates the transition size compared to the modeling results (64 nm). As the elongation of the grain increases (A decreases), the analytical SD-MD critical length grows at a slower rate compared to the modeling results. Thus, for A < 0.79, the analytical solution underestimates the threshold determined via modeling. Mismatches between micromagnetic modeling and analytical results were reported for magnetite (Muxworthy & Williams, 2006) and iron (Muxworthy & Williams, 2015). One suggested explanation for the discrepancy is that analytical solutions assume a uniform SD state transitioning to a two-domain state, which does not accurately represent the "flowering" state seen at larger SD grain sizes or the nonuniform state (Muxworthy & Williams, 2006).
We tested the consistency of our micromagnetic results by running simulations for tetrataenite spheres. The SD-to-two-domain threshold for a sphere is 80 nm diameter, higher than the critical length for a cube (Figure 2b). However, the domain transition occurs at a comparable volume (sphere transition at 102% volume of the cube), highlighting that the transition occurs at a critical volume and the length is therefore shape dependent. We present our results for cuboids to be consistent with previous analytical and micromagnetic modeling work (for example, Butler & Banerjee, 1975a; Muxworthy & Williams, 2015). The analytical solutions follow Equations 2, 3, and 5, and the superparamagnetic-single domain (SD) threshold is presented for uniaxial and shape anisotropy acting in the same direction. Red shapes show tetrataenite grains in the Bacubirito meteorite, with triangles representing SD grains, circles representing two-domain grains, and squares representing grains with greater than two domains (Section 3.3).
(b) Close up of blue-dashed box in (a). Colored rectangles represent cloudy zone tetrataenite islands sizes for selected meteorite groups Bryson et al., 2015Bryson et al., , 2017Bryson et al., , 2019Einsle et al., 2018;Goldstein et al., 2009Goldstein et al., , 2014Goldstein et al., , 2017Maurel et al., 2020Maurel et al., , 2021Nichols et al., 2020;Uehara et al., 2011;Yang et al., 2010). Since images of cloudy zone islands show a range of axial ratios, we allow each meteorite group to occupy the full width of the axial ratio domain, denoted by the doubled-sided arrows. The widths of reported tetrataenite rims are represented by the black bar. Green hexagon shows the SD-to-two-domain transition for spherical tetrataenite based on micromagnetic calculations.

Magnetic Stability of Tetrataenite Domain States
The magnetic stability of a grain is a measure of its resistance to changes in its magnetization state, such as direction or structure, which might be initiated by changing external fields or from viscous relaxation. We tackle the question of tetrataenite magnetic resistance to field and viscous decay separately in the following subsections.

Magnetic Relaxation Times of Tetrataenite
To determine the magnetic relaxation time of SD and two-domain tetrataenite, we calculated NEB paths using MERRILL to quantify the energy barrier between two states. For all tetrataenite grains that fall into the SD stability field (Figure 2), the energy barrier between anti-aligned magnetizations along the magnetic easy axis is associated with a relaxation time >4.5 Ga. For a 30 and 50 nm cube, the NEB path (Figures 4a and 4b) quantifies the energy barrier to be 3,095 and 8,280 k B T, respectively, corresponding to enormous relaxation times of ∼1 × 10 1,344 and ∼9 × 10 3,595 yr, respectively. The magnetization of the grains are flipped via the nucleation of a domain wall in one corner of the grain, which subsequently propagates through the grain, creating a metastable two-domain state in each NEB path (Figures 4a and 4b; Movies S2 and S3). Due to the presence of a metastable state, another NEB path was calculated between the initial SD state and the metastable state to increase the resolution of the energy barrier. Because the energy barrier between the metastable two-domain state and the SD state is ∼1,000 k B T (∼6 × 10 417 yr), a grain geometry that normally falls in the SD stability field (Figure 2) but is placed in a two-domain state will be highly resistant to thermal fluctuations relaxing it to the SD state.
For an 80-nm cube that falls into the two-domain stability field (Figure 2), we modeled the energy barriers between combinations of two-domain and SD states. The energy barrier between an initial two-domain state and a higher-energy SD state is ∼25,000 k B T, though the path involves multiple metastable states (Figure 4f, Movie S4). Recalculating the energy barrier between the two-domain state and the first metastable state, which is a two-domain state with a curved domain wall, results in a ∼9,000 k B T barrier. Modeling the energy barrier between two perpendicularly aligned two-domain states also shows the presence of metastable states with energy barriers consistent with relaxation times >4.5 Ga (Figure 4e, Movie S5). Thus, two-domain tetrataenite is stable against viscous relaxation.
An 80 nm cube placed in an SD state, which has a higher energy than the two-domain state, can transition through multiple metastable states to either an anti-aligned SD state or a two-domain state separately (Figures 4c and 4d; Movies S6 and S7). The first metastable state for both NEB paths is a two-domain state with a curved domain wall that segments the grain into two domains that occupy ∼12.5% and ∼87.5% of the grain volume each (Figures 4c  and 4d inset). The energy barrier between the initial SD state and the metastable state is > 4,000 k B T. Therefore, a grain with a size and shape that would fall into the two-domain stability field but is placed in an SD state is resistant to viscous decay over time periods relevant for paleomagnetic studies.

Simulated Hysteresis Loops
The hysteresis behavior of SD tetrataenite follows Stoner-Wohlfarth theory (Stoner & Wohlfarth, 1948). The grain occupies a uniform state until a critical field is applied antiparallel to the magnetization, resulting in a reversal of the magnetization direction ( Figure 5a, Movie S8). Barring any heating, it is therefore unlikely that a tetrataenite grain in the SD stability field would be found in a two-domain state despite being that state being resistant to thermal relaxation to the SD state (Section 3.2.1). The microcoercivity of tetrataenite increases with decreasing grain size for a constant axial ratio (Figure 6), suggesting that higher fields are necessary to remagnetize smaller tetrataenite grains. The microcoercivity of tetrataenite is also influenced by the axial ratio, with elongated grains requiring higher fields to flip their magnetization direction. This is expected as elongation in the direction of the magnetic easy axis inhibits the ability of the individual element magnetization to coherently rotate in the direction of an external field due to shape anisotropy.
Traditionally, MD grains are considered poor paleomagnetic recorders in part due to the ease of domain states changing in response to external fields by domain wall motion (Nagy et al., 2019). However, two-domain tetrataenite shows resistance to domain wall displacement at external fields strengths relevant for planetary or planetesimal dynamos (<300 μT; Maurel et al., 2021). Application of a 1 mT external field creates a <2% change in magnetization from domain wall movement (Figure 5b). However, the presence of stronger external fields causes the domain with magnetization parallel to the field to grow, displacing the domain wall in a direction perpendicular to the applied field. The field-aligned domain continues to grow at the expense of the antiparallel domain as the field increases until a critical field, H W , at which the domain wall is destroyed, and the grain becomes uniformly magnetized (see Movies S9 and S10). This behavior was observed by Einsle et al. (2018) for tetrataenite in the cloudy zone, who found that the domain wall could be destroyed by island-island interaction fields or an external field, though our results show that the required external field intensity is higher than that expected for most planetary or planetesimal dynamos. The critical field H W is dependent on the axial ratio and grain size, with increasing grain sizes and larger axial ratios possessing higher H W values. For a two-domain tetrataenite grain that transitions to an SD state due to an external field, the reduction of the external field to zero does not place the grain back in a two-domain state. Instead, once the field is eliminated, the initially two-domain grain will continue to occupy a uniform state. This suggests that the uniform state is highly stable even though it is not the least-energy state for the grain geometry. This pattern was observed for all two-domain grains modeled and is a potentially important aspect of cloudy zone tetrataenite magnetization evolution (Section 4.2). Conversely, a SD state cannot be placed in a two-domain state through the application of an external field.

Observations of Two-Domain Tetrataenite
Tetrataenite grains are found in meteoritic plessite, which forms from the decomposition of martensite into an intergrowth of kamacite and taenite/ tetrataenite (Supporting Information S1, Goldstein & Michael, 2006). The coarseness of plessite structures is primarily dependent on the bulk Ni content of the meteorite and its cooling rate (Buchwald, 1975;Goldstein & Michael, 2006). Duplex plessite is a coarse-grained plessite, where taenite/ tetrataenite grains are separated by several μm of kamacite unlike in the nm-scale cloudy zone. No magnetostatic interactions are expected between tetrataenite grains in plessite, but we cannot rule out the existence of exchange interactions between tetrataenite grains and surrounding kamacite based on current knowledge. The tetrataenite grains can be hundreds of nm to a few μm in size.
We collected X-ray photoemission electron microscopy (XPEEM) images of duplex plessite in the ungrouped iron meteorite Bacubirito using Beamline 11.0.1 at the Advanced Light Source, Lawrence Berkeley National Laboratory (Supporting Information S1). Our X-ray circular magnetic dichroic (XMCD) images produced from the XPEEM data (Stöhr et al., 1998) show the magnetization of kamacite and individual taenite/tetrataenite grains along the X-ray beam direction (Figure 7). We applied a magnetic field of increasing intensity in steps of ∼50 mT up to 324 mT and collected XMCD images at each step. We find that the uniformly-magnetized individual taenite/ tetrataenite grains have coercivities between 250 and 324 mT ( Figure S4 in Supporting Information S1 and Figure 5). Moreover, the grains have an average 48 wt% Ni content (measured using electron dispersive spectroscopy). While Ni composition alone cannot distinguish between taenite and tetrataenite, reported coercivities for taenite formed from the annealing of ∼30 nm tetrataenite in the Santa Catharina ataxite are less than 10 mT (Dos Santos et al., 2015). Therefore, we interpret the measured coercivities of the Ni-rich grains to be indicative of tetrataenite and not taenite.
XMCD images show that the kamacite is MD (Figure 7). The tetrataenite grains are either SD or MD, though some grains have XMCD signals that cannot be distinguished from the kamacite. Four grains observed show uniform magnetization and have lengths and axial ratios that place them in the two-domain stability field in Figure 2a close to the SD-two domain transition according to analytical calculations. There are multiple reasons why these grains appear SD but fall into the two-domain stability field: (a) The analytical calculations do not accurately reflect the SD-two domain transition, which is suggested by the mismatch between the MERRILL results and the analytical calculations. Grains of this size and geometry were unable to be modeled using MERRILL due to computational limits; (b) The grains were initially two-domain but later placed into an SD state due to an external field; (c) The grains are actually in a two-domain state, but the cut surface of the sample only allows the observation of one domain. This could be due to the easy axes and domain wall planes of the grains lying parallel the sample surface yet too far under the surface to be observed with XPEEM. XPEEM only provides the magnetization from the top 5 nm (Bryson, Herrero-Albillos, et al., 2014) of the sample.
A two-domain grain with a length of 400 nm and axial ratio of 0.57 was observed (Figure 7), showing that the two-domain state is achievable in nature under non-interacting conditions and relevant grain geometries. The size and axial ratio of the grain places it in the two-domain stability field ( Figure 2). Moreover, finding a two-domain state implies that the meteorite was not substantially remagnetized at fields above the coercivity of the grain (Section 3.2.2). This rules out the possibility that the aforementioned SD grains were initially in a two-domain grain and had their domain wall destroyed by an external field.
Other grains with clear XMCD signals are in other MD states, consisting of more than two domains or having more complex domain structures than the described two-domain grain in the previous paragraph (Figure 7). One grain has a size and shape comparable to that of neighboring uniformly-magnetized grains, while two other MD grains are very elongated (A < 0.25) and are larger than the other grains imaged (Figure 2). The observation of an MD grain of a size similar to the size of uniformly magnetized grains is inconsistent with our understanding of individual tetrataenite domain states as a function of grain geometry. We see three possible explanations: (a) Exchange interactions with the kamacite matrix that are currently incompletely understood; (b) Our micromagnetic modeling does not account completely for the geometry of the grain. Note that the domain threshold sizes are dependent on the morphology of the grain (Section 3.1); (c) The grains are larger than they appear on the cut surface of the sample.  XMCD images of the tetrataenite grains after the application of a 324 mT field show that the magnetization of the majority of the grains are aligned in a similar direction relative to the X-ray beam (Figure 7c). The uniformly-magnetized grains still retain a uniform magnetization and the higher-domain MD grains, while still showing multiple domains, appear to be mostly magnetized in the same direction as the SD grains. However, the previously identified two-domain grain shows a complex domain structure consisting of three to four domains with curved domain walls. This is unexpected based on hysteresis loops produced by micromagnetic modeling, which suggested that the grain should remain in a two-domain state until the domain wall is destroyed, creating an SD structure. We do see metastable states with curved domains in our NEB paths in Figure 4, and therefore propose that the grain is in a metastable state or that our micromagnetic modeling cannot account for all interactions between the external field and the grain magnetization.

Ideal Tetrataenite Grain Geometry
Tetrataenite is stable against viscous relaxation and external field remagnetization in both its SD and two-domain states. The size range of the SD stability field varies with axial ratio; for equant grains it ranges from 6 to 64 nm. Below 6 nm, tetrataenite is in a SP state. Above 64 nm, tetrataenite occupies a two-domain state. However, in the presence of external field on the order of 100 mT, a two-domain grain can be placed into a uniform state through the destruction of its domain wall (Section 3.2, Figure 5). Additionally, the magnetization of a two-domain grain or a SD grain in the two-domain stability field is stable over the lifetime of the solar system. Therefore, we propose that isolated, non-interacting tetrataenite is stable even at sizes above the SD range in the absence of heating or strong external fields, though we are unable to conduct micromagnetic models to place an upper bound on the stability range once tetrataenite forms more than two domains.
Tetrataenite grains in cloudy zones, the adjacent rims, and in plessite range in size from tens of nm to ∼10 μm (Goldstein et al., 2009(Goldstein et al., , 2014Nichols et al., 2020Nichols et al., , 2021. We can compare the grain sizes reported in previous studies with the SP-SD-two domain thresholds (Figure 2) to determine their magnetic stability. Among the meteorite groups that have been the subject of paleomagnetic studies, the fast-cooled IVA irons are the only group with cloudy zone tetrataenite grains all falling in the SD region (Nichols et al., 2020). Two other meteorite groups, the IVB and IIIAB irons, have islands that fall exclusively in the SD stability field but have not been studied paleomagnetically. Otherwise, tetrataenite grain sizes in cloudy zones fall either partly inside or entirely outside the SD stability field. For example, the IAB iron meteorites possess island sizes that fall partly in the SD range for elongated grains while the mesosiderites contain islands that fall exclusively outside the SD stability field. Mesosiderite tetrataenite islands are reported to exhibit MD behavior (Nichols et al., 2020) and IAB islands close to the clear tetrataenite rim appear SD (Bryson, Herrero-Albillos, et al., 2014). However, the spatial resolution of XPEEM (30 nm; Stöhr et al., 1998) prevents the determination of individual island domain states.
While we are unable to conduct micromagnetic modeling of tetrataenite grains with least energy states containing more than two domains, a combination of analytical calculations and observations from imaging of tetrataenite domains in the tetrataenite rim abutting cloudy zones and large tetrataenite grains in plessite can constrain the extent of the two-domain region (Figure 2). The tetrataenite rim ranges in width from 200 nm to 10 μm depending on the cooling rate Uehara et al., 2011). Faster cooled meteorites, such as the Steinbach and Chinautla IVA irons, are reported to have 200-nm rims and possess magnetic domains <2 μm long (Bryson et al., 2017;Nichols et al., 2020). In comparison, the slow-cooled Estherville mesosiderite has a 5 μm-wide clear tetrataenite rim with magnetic domains that are ∼200-300 nm wide and up to 3 μm long (Nichols et al., 2020). Analytical calculations predict that the transition from a two-domain to a three-domain state occurs at ∼360 nm for cubic tetrataenite grains, which could be compatible with the 200-300 nm width of the observed domains in the tetrataenite rim, though the axial ratio of the tetrataenite rim is difficult to define as the rim runs along the full length of the cloudy zone.
Akin to the tetrataenite rim, large tetrataenite islands and individual tetrataenite grains outside of the cloudy zone microstructure can be MD (Nichols et al., 2020;Uehara et al., 2011, this study). For example, the Estherville mesosiderite has islands as large as 463 ± 32 nm (Goldstein et al., 2014). Such sizes fall in the two-or three-domain stability field depending on the axial ratio ( Figure 2) and are reported to be MD (Nichols et al., 2020). Ordinary chondrites contain zoneless plessite grains with μm-sized tetrataenite grains and are interpreted to be MD as well based on the low coercivity of remanence of the grains compared to the cloudy zone (Uehara et al., 2011).

Role of Island-Island Interactions in Cloudy Zone Tetrataenite
In our micromagnetic modeling, we focused on individual tetrataenite grains and ignored interactions between grains to elucidate the natural domain states of tetrataenite. However, the proximity of tetrataenite islands to one another (<1-2 island diameters) in the cloudy zone suggests that island-island magnetostatic interactions should be considered. Micromagnetic modeling of cloudy zone islands from the Tazewell IAB iron meteorite (Einsle et al., 2018) showed that the transition of taenite to tetrataenite through an increase in the magnetocrystalline anisotropy constant forced a taenite SV state into a tetrataenite grain with two oppositely oriented magnetization domains and a domain wall. The two-domain structure is similar to that seen in our simulations. Magnetostatic interactions between the islands act to displace the grain wall until it is destroyed, creating SD islands. This process is akin to the application of an external field in the hysteresis loops shown in Figure 5, which forces out the domain wall at H W . Treating a cubic tetrataenite island of volume v like a simple dipole with a magnetization m equal to its saturation magnetization, the radial component of the external field B created by that grain at a distance r along the dipole axis follows At the center of a neighboring island of the same size, that field produce by the tetrataenite cube is ∼275 mT, though the field could be larger or smaller depending on the geometry of the island and its neighbor. This field is the same order of magnitude as H W for the grains we modeled, indicating that magnetostatic interactions could indeed provide a method of destroying the domain wall in a two-domain grain. Initial modeling of interacting 70 nm tetrataenite cubes separated by 1 nm shows that if one grain is initialized as SD and the other begins in the two-domain state, the final equilibrium assemblage consists of the SD grain plus a two domain grain with the domain wall shifted such that the domain parallel to the SD domain direction grows ( Figure S6 in Supporting Information S1). Notably, the domain wall is not completely destroyed. As the distance between the grains grow, the degree to which the domain wall shifts is lessened. Generally, the strength of the island-island interactions remains to be quantified. The maximum field that we applied to the non-interacting tetrataenite grains in the Bacubirito plessite is the same order of magnitude as H W .
Our results support the idea that magnetostatic interactions between islands are of fundamental importance in the ability of cloudy zones to preserve the record of ancient magnetic fields. As the external field intensity necessary to completely destroy the domain wall is on the order of hundreds of mT, which is orders of magnitude stronger than any reported planetesimal dynamo (less than a few hundred μT; Maurel et al., 2021), the identification of uniform magnetization across cloudy zone islands is likely solely due to magnetostatic interactions. If no magnetostatic interactions were present, then grains above the SD-two domain threshold would be expected to be two-domain. However, the stability of the SD state in grains above the SD-two domain threshold suggests that if tetrataenite falls into a metastable SD state during ordering from taenite, it will likely remain in the SD state. While the two-domain state is stable, additional work must be conducted to determine how to interpret the magnetic record of this domain state. If an external field is present, it could bias the destruction of domain walls via island-island interactions in a certain direction, creating the preference of islands toward one easy axis as seen in XPEEM images (Bryson, Herrero-Albillos, et al., 2014).
We propose that magnetostatic interactions have a direct influence on the ability of cloudy zones to be suitable for paleomagnetic analyses and allow for slower cooling rates that would otherwise lead to islands sizes above the SD-two domain threshold. This proposal is in contrast to the conclusions of Nichols et al. (2020) that magnetostatic interactions do not have an influence on determining which cloudy zones can be the subject of a paleomagnetic study since our results show that magnetostatic interactions are fundamental in placing the islands in an SD state. However, like Nichols et al. (2020), we agree that magnetostatic interactions do not play a role in affecting whether the cloudy zone shows evidence of a past magnetic field since the destruction of the domain wall is likely due to island-island interactions based on the strength of the field necessary to destroy the domain wall.

Comparison of Predicted and Observed Tetrataenite Coercivities
As detailed in Section 3.2.2, the coercivities of individual tetrataenite grains based on our simulated hysteresis loops range from 1 to 2.6 T depending on the size and axial ratio ( Figure 6). The trend of decreasing coercivity with increasing grain size is consistent with observations of tetrataenite macrocercivities from Bryson, Church, et al. (2014) and Uehara et al. (2011). However, the reported coercivities of tetrataenite islands in the cloudy zone (∼1 T or less) that are within the size range encompassed by our models are lower than predicted from our hysteresis loops ( Figure 6). This suggests that either our model is not predicting the correct coercivities or that interactions between islands or between islands and the surface of the antitaenite matrix in cloudy zones have an effect on the island coercivities that is not included in our micromagnetic modeling. The proposal that the matrix affects island coercivities is supported by micromagnetic simulations conducted on cloudy zone tetrataenite in a "soft" magnetic matrix which leads to a decrease in coercivity as a results of exchange coupling between the islands and the matrix (Bryson, Church, et al., 2014). While antitaenite has unambiguously been identified as the matrix phase between tetrataenite island through Mössbauer spectroscopy (Blukis et al., 2017), the paramagnetic nature of the phase is inconsistent with the magnetically-soft ferromagnetic phase seen in XPEEM images (for example, Nichols et al., 2020). One potential explanation for this discrepancy is that the stress at the surface of the cut face of samples could allow a magnetic structure to stabilize in a phase that is otherwise non-magnetic (Blukis et al., 2017). If this is the case, it is uncertain whether the tetrataenite islands would have a larger effect on the stable magnetic structure of the matrix phase, or if the matrix phase impacts the islands to a greater extent.
The large grain sizes of the tetrataenite in the Bacubirito plessite do not allow for direct comparison between our modeled and observed coercivities. Assuming the curves can be propagated at larger sizes using a Lorentzian fit, we can provide a tentative estimate of the coercivity provided for larger grains by our model ( Figure S1 in Supporting Information S1). At 400 nm, predicted coercivities range between ∼850 and ∼1,050 mT depending on the axial ratio. These are higher than the observed coercivities of the uniformly magnetized grains, the majority of which fall between 250 and 324 mT. The difference between the observed and predicted coercivities could indicate some interaction between the tetrataenite and the matrix kamacite as suggested for cloudy zone tetrataen ite, but we caution that there is no guarantee that the size-coercivity relationship seen in our models is consistent out to such large grain sizes. Our modeled coercivities only span the SD stability field and the lower grain sizes of two-domain stability field. These trends might not be applicable to larger two-domain grains and higher-domain grains, especially as the coercivities appear to plateau at larger grain sizes in our model.

Using Individual Tetrataenite for Paleomagnetic Studies
Obtaining accurate records of the intensity and direction of past magnetic fields from isolated, non-interacting tetrataenite grains (e.g., in plessite) remains to be solved. One limitation is that while XPEEM experiments can image the magnetization direction of tetrataenite with high-spatial resolution, the large spacing between grains means that the number of grains necessary to build a statistically meaningful data set (∼10 4 for a 1° uncertainty in the recovered field paleodirection; Berndt et al., 2016) can make XPEEM time consuming and labor intensive. For Bacubirito duplex plessite, we can capture about ∼100 grains per XPEEM image (not all of which are SD), compared to the many hundreds of grains in cloudy zone XPEEM images and several thousands of grains imaged over the course of a cloudy zone study.
If the tetrataenite grains are dispersed amongst non-magnetic materials such as silicate grains (Figure 1c), the magnetization could be studied using a superconducting moment magnetometer, a superconducting quantum interference device (SQUID) microscope (Lima & Weiss, 2016), or a quantum diamond microscope (Glenn et al., 2017) depending on the strength of the magnetic moment. However, additional work needs to be conducted to determine how to interpret tetrataenite paleomagnetism in a bulk demagnetization study, which is complicated by the uniaxial anisotropy of tetrataenite and the importance of crystallographic orientation in magnetization direction. This is different from cloudy zone paleomagnetic studies that take advantage of the effect of external fields on island-island interactions to bias the similarly crystallographic-oriented cloudy zone in the directly of the external field (Nichols et al., 2020).

Conclusions
Tetrataenite is a potent magnetic recorder that to date has only been utilized in paleomagnetic studies through investigations of the cloudy zone microstructure. However, tetrataenite can occur as individual, non-interacting grains in a variety of meteorite groups including irons, ordinary chondrites, and primitive achondrites. Understanding the fundamental magnetic properties and domain behavior of tetrataenite is not only important in order to further our understanding of the cloudy zone, but also to take a step toward expanding our ability to conduct paleomagnetic investigations of meteorites using individual tetrataenite grains.
Here, we describe the domain states of non-interacting tetrataenite as a function of grain geometry based on fundamental magnetic principles and micromagnetic modeling. Tetrataenite in the size range of ∼10-∼160 nm occupies a uniform, SD state dependent on the grain's axial ratio. Elongated grains retain the SD state at larger sizes compared to equant grains. SD tetrataenite grains possess relaxation times greater than the lifetime of the solar system. Below 6 nm, tetrataenite occupies a SP state. Once tetrataenite is too large to occupy a uniform state, the least-energy domain configuration is a two-domain state with the domains and domain wall parallel to the direction of uniaxial anisotropy. This is observed in XMCD images of tetrataenite in the ungrouped Bacubirito iron meteorite. Importantly, tetrataenite cannot occupy a vortex state due to the high uniaxial anisotropy of the mineral. The high uniaxial anisotropy also allows the two-domain state to be stable over the lifetime of the solar system.
Through the application of an external field, the two-domain state can be forced into a uniform state through a destruction of the domain wall. The subsequent uniform state can be retained even after the external field is removed. The majority of paleomagnetic cloudy zone studies rely of tetrataenite islands, whose sizes largely fall in the two-domain state range. However, magnetostatic interactions between islands provide a method of displacing the domain wall and creating single domain tetrataenite. Therefore, island-island interactions play a vital role in making the cloudy zone a useful paleomagnetic recorder.