Hysteresis of Natural Magnetite Ensembles: Micromagnetics of Silicate‐Hosted Magnetite Inclusions Based on Focused‐Ion‐Beam Nanotomography

Three‐dimensional geometries of silicate‐hosted magnetic inclusions from the Harcus intrusion, South Australia have been determined using focused‐ion‐beam nanotomography (FIB‐nt). By developing an effective workflow, the geometries were reconstructed for magnetic particles in a plagioclase (162) and a pyroxene (282), respectively. For each inclusion, micromagnetic modeling using Micromagnetic Earth Related Rapid Interpreted Language Laboratory provided averaged hysteresis loops and backfield remanence curves of 20 equidistributed field directions together with average Ms, Mrs, Hc, and Hcr. The micromagnetic structures within each silicate are single‐domain (SD), single‐vortex (SV), multivortex (MV) and multidomain states. They have been analyzed using domain‐state diagnostic plots, such as the Day plot and the Néel plot. SD particles can be subdivided into groups with dominant uniaxial anisotropy (Mrs/Ms ∼0.5 and 10 < Hc < 100 mT) and mixed uniaxial/multiaxial anisotropy (Mrs/Ms ∼0.7 and 10 < Hc < 30 mT). Most SV particles lie on a trend with 0 < Mrs/Ms < 0.1 and 0 < Hc < 10 mT, while others display a broad range of intermediate Mrs/Ms and Hc values. SV and MV states do not plot on systematic grain‐size trends. Instead, the multicomponent mixture of domain states within each silicate spans the entire range of natural variability seen in bulk samples. This questions the interpretation of bulk average hysteresis parameters in terms of grain size alone. FIB‐nt combined with large‐scale micromagnetic simulations provides a more complete characterization of silicate‐hosted carriers of stable magnetic remanence. This approach will improve the understanding of single‐crystal paleomagnetism and enable primary paleomagnetic data to be extracted from ancient rocks.

exchange constant, the magnetocrystalline anisotropy constant, M s the saturation magnetization, and μ 0 the vacuum permeability. For magnetite the intermediate size range lies roughly between 0.070 and 0.200-0.300 μm, however, the formation of clearly defined homogeneously magnetized domains does not occur below grain sizes of about 2 μm. Nagy et al. (2019) modeled grains up to the micron range and observed domain wall formation in unconstrained micromagnetic models. What is missing, is a comprehensive survey of magnetization states and predicted hysteresis properties for naturally occurring magnetite particles across, and beyond, the SD state. Here, we develop a method to achieve such a survey based on focused-ionbeam (FIB) tomography combined with micromagnetic modeling.
To compare our model results to measurements on sized magnetite samples, we consider several isothermal measurements that have been developed to quantify the characteristic domain states of natural remanence carriers. The most common approach is based on the measurement of bulk average magnetic hysteresis parameters, such as saturation magnetization M s , saturation remanent magnetization M rs , coercivity H c , and coercivity of remanence H cr . Two intensely studied diagrams are the Day plot of M rs /M s versus H cr /H c (Day et al., 1977), and the Néel plot of M rs /M s versus H c (Néel, 1955), which are based on bulk hysteresis parameters that can be quickly and easily measured in most rock magnetic laboratories (D. Dunlop, 2002;Tauxe et al., 2002). It is well known that their interpretation is ambiguous without detailed knowledge of the underlying magnetic ensemble (Roberts et al., 2018). To relate domain state to geometrical measurements, the theoretical diagram of Butler and Banerjee (1975), and later revisions (Fabian et al., 1996;Muxworthy & Williams, 2006), is commonly used to characterize the domain state of magnetic particles based on particle length and axial ratio. This approach assumes that magnetic particles can be effectively approximated as prolate ellipsoids.
Here, we apply our workflow and method to magnetite particles within silicate-hosted magnetic inclusions. This setting is abundant and important for studying the Earth's magnetic field, magmatic evolution, geochemical interaction, stability, and oxygen fugacity (Feinberg et al., 2005;Tarduno et al., 2020). Magnetite is the most significant Fe-oxide for paleomagnetic studies, and magnetite is known to form as exsolved or included particles in silicates like plagioclase (Davis & Letters, 1981;Feinberg et al., 2005;Usui et al., 2015) and pyroxene (Feinberg et al., 2004;Fleet et al., 1980;Frandsen et al., 2004;Renne et al., 2002) within a wide range of geologic settings. Davis and Letters (1981) document that magnetite exsolves several micrometer long rods in 2-4 different directions in plagioclase from oceanic gabbros. The crystallographic relationship between magnetite needles and a plagioclase host was determined by Wenk et al. (2011). They observed a single set of magnetite needles oriented parallel to the (001) axis of plagioclase, with the needle elongation direction corresponding to the (110) direction of magnetite and the (111) planes of magnetite oriented subparallel to (120) and (120) planes of plagioclase. Within pyroxene (CPX) iron oxide inclusions occur parallel to (010), where all inclusions are elongated and with the long axis subparallel to both (100) and (001) (Feinberg et al., 2004;Renne et al., 2002). In orthopyroxene exsolved oxide, lamellae are known from optical and transmission electron microscopy. Ilmenite lamellae are present as coherent rods and blades where (001) ilmenite is epitaxially intergrown on (100) pyroxene, and the a-axis of ilmenite is parallel to the c-axis of pyroxene . There the ilmenite contained hematite lamellae parallel to the (001) of the ilmenite. The incorporation of iron into the silicates is not confirmed, but suggested to occur at primary subsolidus mineral formation and exsolve magnetite well above its curie temperature (Feinberg et al., 2004;Fleet et al., 1980) for pyroxene. The exact process for exsolving magnetite in plagioclase is so far not established, but once the orientation relationship is established the phase boundary between plagioclase and magnetite can be established (Wenk et al., 2011).
The advent of high-resolution FIB nanotomography (FIB-nt), combined with the development of powerful finite element micromagnetic (FEM) software optimized for rock magnetic applications (Fabian & Shcherbakov, 2018;Ó Conbhuí et al., 2018), provides a new opportunity to perform domain-state diagnosis from first principles. FIB-nt is performed on a dual-beam electron microscope that combines the high-resolution imaging capabilities of a field-emission gun scanning electron microscope (SEM) with the nm-precision slicing capabilities of a FIB (Brogden, 2015;Fagerland, 2014). This technique is widely used to analyze biological samples (Beckwith et al., 2015;Guehrs et al., 2017;Mulders et al., 2006), alloys (Cao et al., 2009;Ding & Jones, 2011), and geological materials (Einsle et al., 2016;Lascu et al., 2015;ter Maat et al., 2020;Warr et al., 2014). The process is destructive, as the slicing destroys the prepared sample area. However, it has higher spatial resolution than comparable nondestructive methods of X-ray computed tomography, enabling the morphology of particles spanning the SD to MD range to be reconstructed in 3D. Rather than attempting to solve the inverse problem of unmixing domain states from a series of magnetic hysteresis measurements, the combined FIB-nt-FEM approach enables forward models of domain states and hysteresis properties to be calculated from direct 3D observations of a representative particle ensemble. Evaluating magnetic signatures in geologic samples will always be challenging. The complexities of particle size, shape, spacing, chemical composition, stoichiometry, and stress will generate unique mineralogical settings, and it is a fair question raised by Lascu et al. (2015) whether domain-state diagnosis might be an unachievable ideal? Here, we show that FIB-nt-FEM is able to resolve domain-state ambiguities associated with the size, shape, and spacing of particles. In combination with other techniques, such as energy-dispersive X-ray spectroscopy (EDS) and high-resolution electron backscattered diffraction (EBSD), it has the potential to resolve ambiguities associated with chemical composition, stoichiometry, and stress in the future.

Samples
This study focuses on two gabbroic norite samples obtained from drill cores of the Mount Harcus Intrusion, which is a part of the Giles Complex, South Australia. The sample names 311.5 and 109.3 in Figure 1 correspond to their depth in m. The Giles-Event intrusions were emplaced into the Proterozoic Musgrave Province at about 1,090-1,040 Ma (Maier et al., 2015). The Harcus intrusion is of interest because of its striking remanent magnetic anomaly and the potential for magmatic nickel sulfide mineralizations (Austin et al., 2014). The high-coercivity remanent magnetization is interpreted to be carried by SD magnetite . In reflected light, the two thin sections in Figure 1 display a wide distribution of magnetite grain sizes, from large, discrete (>100 μm) MD magnetite with oxy-exsolution of ilmenite lamellae, and spinel needles and blades. Reduction exsolution of magnetite from ilmenite in the form of blades occurs in isolated ilmenite grains, and in some large ilmenite lamellae, which were first oxy-exsolved from magnetite Robinson et al., 2016). Silicate-hosted micrometer-to nanometer-sized inclusions of magnetite are found in feldspars, pyroxenes, and amphiboles (≪0. highest concentrations of iron-oxide inclusions and were targeted for 3D imaging using FIB nanotomography. Both thin sections 311.5 and 109.3 are composed of plagioclase, pyroxene, amphibole, intragranular quartz, and discrete iron-titanium oxides. Iron-oxide inclusions are present in primary plagioclase and pyroxene and secondary amphibole. Sample 311.5 has a high concentration of inclusions in plagioclase and amphibole. Metamorphic processes have possibly redistributed iron oxides in amphiboles. The intragranular quartz in association with amphibole also suggests the presence of fluids that may have mobilized Fe and removed it from the system. Sample 109.3 has a lower amount of amphibole and iron-oxide inclusions in plagioclase balanced by an increase of iron-oxide inclusions in the pyroxene. For samples 311.5 and 109.3, the bulk natural remanent magnetization values are 34 A/m and 4 A/m, and their bulk magnetic susceptibility values are 0.06 (SI) and 0.10 (SI), respectively.

Focused-Ion-Beam Slice-and-View Nanotomography
Raw data collection was performed using the Auto Slice & View 4 software on an FEI Helios G4 UX Dual-Beam FIB microscope at the NTNU NanoLab in Trondheim, Norway. Polished thin sections were coated with 10 nm of Au and the regions of interest were analyzed with EDS to confirm that the target iron oxide inclusions were magnetite with minor (<1 wt%) Ti. At the optimal imaging settings for the selected minerals, the electron beam drift was about 500 nm an hour, which limited the size of the block volume that could be successfully imaged. A block volume of material measuring 12 × 14 × 12 μm (L-W-H), sliced with 30 nm increments, took, on average, 20 h. With longer sessions, the failure rate increased substantially. Block volumes were also limited by hardness of the material and the ability to successfully prepare the trenches and clean the front face of the sliced volume. The extracted volume of pyroxene is therefore smaller than for plagioclase since the ion-beam has fewer issues with removing material.
Sample preparation for slice-and-view was performed by manual instrument operation, outside of the automated software. The selected area was first coated with 500 nm of Pt. Five 200 nm deep lines that extend the full depth of the milling area were cut into the Pt coating and then filled by carbon deposition using the ion beam. The selected area was then covered with a further 300 nm of Pt. The central three carbon rods were orientated parallel to the slice direction and used as image alignment references. The outside two carbon rods were oriented at an angle of ±60° to the slice direction and were used to calibrate the slice thickness and to enable alignment of images in the slice direction in cases where there are missing images or if a milling session had to be restarted. The surface preparation was finalized with 3-μm wide and 500-nm thick carbon rectangle, oriented parallel to the slice direction, and covered with an additional 300 nm of Pt. This carbon slab is used to keep the E-beam image aligned with the front face at all times. Trenches, 15-μm deep, 10-13-μm wide, and 10-μm long, were created to the left and right of the target volume, and a trench 20-25-μm wide and 20-μm long was created in front of the target volume using a standard cross section milling routine at 30 keV and 20 nA. The front and sides were then polished with cleaning cross sections, first at 1.2 nA and then with a final cleaning step at 0.75 nA. A Pt fiducial mark for ion beam alignment was placed on the thin section surface outside of the milled trenches as seen in I-beam and E-beam images in Figure 1.
For slicing, the Ga-ion beam was kept at 30 keV and 0.75 nA and slice thickness was chosen between 10 and 30 nm. Electron-beam imaging applies a mirror detector for the backscattered-electron signal in immersion mode at 5 kV and 0.8 nA. Tilt correction for the 52° angle between the ion-and electron-beam was performed in the FEI software with dynamic focus by activating the tilt-correction option. Each scan with 6144 × 4096 pixels at 10-μs dwell time, resulted in a 1535 × 1652 image with pixel size of 7.5 nm in the pyroxene, a 2037 × 2043 image with pixel size of 7.6 nm in plagioclase. Early trials included the automatic z-axis drift equal to the milling depth and used a fiducial mark outside the milling area for consistent E-beam alignment. However, because this instrument has a constant E-beam drift that is not accounted by the software, it was found to be much more efficient to use the carbon slab as a fiducial mark that remains inside the milling area. This approach reduced the need for additional automated focusing steps after each slice. The high contrast between platinum and carbon reduces the risk of losing the fiducial mark as the project runs and increases the success rate of the automated procedure.

Image and Object Filtering
The raster images of the milled front faces are provided in Tagged Image File Format and were stacked in Fiji (ImageJ) and aligned with the stack alignment plugin (Tseng et al., 2011), using the cross section of the three central carbon trenches as reference pattern. Using their high electron reflectivity and thus brightness, magnetite particles were extracted by gray scale thresholding and exported as binary images. Based on backscatter intensity and cross-referenced with the EDS analyses, anything that we could not confirm as magnetite was segmented out. In the case of complex structures and overlapping threshold values, each slice was manually "cleaned" by going through the binary and image stack (side by side), removing any artifacts, or segmented threshold values that were not magnetite. For the pyroxene stack, the thresholding was particularly difficult because the bottom sections of each image had lower overall backscatter intensity than the upper parts. This image stack had to be treated with a nonlocal means filter and was segmented using the Dragonfly software from ORS. After the binary segmentation, the stack was passed through 3 × 3 × 3 pixel 3D-mean and Gaussian-blur filters to smooth over the pixel boundaries and to reduce sharp edges in objects at the boundary of the originally cropped volume. Areas with objects in close proximity to each other were manually edited to prevent bridging artifacts when converted to a finite element surface mesh. We excluded incomplete particles at the edge of the milled volume and particles that are too small to be confidently distinguished from background and noise.

Mesh Generation
The stack generated in Fiji was imported to Paraview (Ahrens et al., 2005), then resized based on pixel length, width, and slice thickness to regenerate the actual size of the milled volume. A surface mesh of the magnetic particles was generated using a contour filter based on the magnetite brightness threshold. The resulting initial mesh was exported as a stereolithography (STL) file. Both, Meshmixer (Schmidt & Singh, 2010) and Meshlab (Cignoni et al., 2008) were then used to improve the meshes. The quality of each object was evaluated individually. An object was discarded if it did not intersect at least two slices. This imposes a lower-size limit of approximately 30 nm for particles to be modeled in 3D using this method. For geological time-scales, this corresponds approximately to the SP transition size for magnetite. The geometric mesh quality for each object was checked, and any holes filled to generate a closed mesh surface. Surface mesh node-density was reduced to generate a more smoothed surface that minimizes the appearance of artificial steps reflecting the finite distance between the measured slices. After eliminating the voxel steps between the slices, mesh node-density was again increased to generate a now smoothed mesh with a target edge length ≤8.9 nm, the exchange length of magnetite. Target edge length, mesh quality, and volume were checked in Meshlab. If the mesh still contained artificial sharp edges due to the object's small size, it was remeshed and smoothed in Meshmixer to fit the original size. Iso2mesh and a Matlab script (Fang & Boas, 2009) were used to generated tetrahedral volume meshes from the surface meshes. Also here a target edge length ≤8.9 nm was used. To keep the node density similar for each object, the edge length was reduced for smaller objects. The output file was converted to a Patran file using Git Bash (Windows computer) to call a shell script (Convert2Pat; Ó Conbhuí et al., 2018). The Patran file is then called by a script to be used by Micromagnetic Earth Related Rapid Interpreted Language Laboratory (MERRILL) for computing the magnetization structures. Aspect ratios of the final meshes were computed from the STL files based on the volume inertia tensor, by computing its eigenvalues and eigensystem by a Python routine. To ensure correct calculations, the proportional relationship between the eigenvalues was tested by comparing it to the measured short and long axis of the same particles in Meshlab.

Micromagnetic Modeling
Micromagnetic modeling uses the MERRILL, a micromagnetics package optimized for rock magnetism (Fabian & Shcherbakov, 2018;Ó Conbhuí et al., 2018). MERRILL applies finite element tetrahedral meshes to calculate local micromagnetic energies, and a boundary element method to calculate the demagnetizing energy and its gradient (Ó Conbhuí et al., 2018). For each particle in our study, first an initial remanent domain state was obtained by minimizing the total micromagnetic energy starting from a state of fully randomized spins. Next, for each particle, the upper branches of 20 hysteresis loops between the maximal fields of 350 and −350 mT was calculated in 5-mT steps. Each of these hysteresis loops was calculated for one of 20 directions chosen from a Fibonacci sphere to achieve an approximately homogenous axis distribution. Backfield remanence curves from 0 to −200 mT in 5-mT steps were calculated for the same 20 field orientations starting from a positive saturation remanence state. While it is possible to run multiple MERRILL operations simultaneously, with the computer power available (Linux server with 20 cores and 128 GB memory), it was still necessary to limit the input file size because of the exponentially increasing time used for minimization. For file sizes < 20 MB, an average hysteresis run took 15 min and an average backfield curve required 35 min. It was possible to complete one hysteresis loop for 250 particles in two field directions within 24 h. The total data set presented in this study equals to approximately 8,000 hysteresis loops and 8,000 backfield curves. The numerical average of the hysteresis loops and backfield curves over the 20 field directions for each mesh geometry approximately represents a random spherical ensemble of equal particles. Values of M rs , H c , and H cr were extracted from the average curves. Visual presentation of the magnetization structures of individual models was obtained using the open-source software ParaView (www.paraview.org; Ayachit, 2015).

Plagioclase
One cube measuring 12.77 × 14.50 × 12.77 μm (L × W × H) was extracted from plagioclase (PLAG [1]; Figures 1 and 2). This volume was reconstructed from 421 slices, each with a thickness of 30 nm. We extracted a total of 162 particles of magnetite and were able to run hysteresis loops for 140 of these. In this volume of plagioclase, the magnetite forms two different patterns differentiated by shape and size. The first and most noticeable are magnetite "rods" that are a few to several micrometers long and often extend beyond the milled volume. In Figure 2, we observe these rods oriented in two sets of parallel directions. This observation suggests that these rods are crystallographically oriented, however, we would need to cross reference EBSD and magnetite long axis pole plots to confirm this. Cross sections show that these rods/needles either form flattened or prismatic shapes, with irregular gaps (red arrow in Figure  reasonable time (one smaller 6.6 × 0.2 μm magnetite rod was successfully minimized in zero-field from an initial state of randomized moments and shown in Appendix A1). Appendix A1 shows the general observation for magnetite rods that the magnetization aligns with the length of the rod and generates vortices at each end. The 140 smaller magnetite particles are a mix of prolate, oblate, and spherical particles with 0.71 median aspect ratio and 0.131-μm median equivalent sphere diameter ( Figure 4). These are commonly associated with an unknown mineral with a lower backscatter intensity that was removed during electron reflectivity segmentation, as described in Section 2.3 (Appendix A2).

Pyroxene
The second volume is from pyroxene (PYR [2] in Figure 1) and measures 6.5 × 10.24 × 9.73 μm (L × W × H; Figure 3). This volume is reconstructed from 647 slices with a slice thickness of 10 nm, which contained total of 282 magnetite particles of which hysteresis loops were calculated for 273. Magnetite in pyroxene forms elongated particles with a median 0.31 aspect ratio (Figure 4b) that are observed in three orientations. This silicate-oxide relationship is described by (Fleet et al., 1980;Frandsen et al., 2004;Renne et al., 2002). For the small magnetite particles included in this pyroxene data set, the Ti content is expected to be low. The largest object does show sections that contain higher Ti values and is excluded from this data set. EDS overview maps for areas with small particles also contain trace amount (<1 wt%) Ti, but there is no visible change in the backscatter intensity, nor any Ti peak in EDS spot analysis and we can therefore not confirm where the Ti is located. Due to the size of the smaller objects and the spot size of the SEM, there is a possibility that the magnetite contain sections that are higher in Ti (ul-vöspinel or ilmenite). TEM or microprobe is needed to quantify the Ti content and confirm the mineral identification. Based on the cited papers above, there is good confidence to identify these smaller inclusions as magnetite.
NIKOLAISEN ET AL. 7 of 20 10.1029/2020GC009389  Figure 4 and Appendix A5 summarize sizes (A) and aspect ratios (B) for the modeled particles in both plagioclase and pyroxene. The "equivalent sphere diameter" is calculated from a sphere with equivalent particle volume. The aspect ratio was calculated using the ratio of min over max eigenvalues of the moment of inertia. The complete data set presented in Figure 4 is subdivided according to both the type of silicate host and the zero-field-minimized magnetic state ( Figure 5 and Table 2). For SD particles, we observe a median size difference between particles hosted by plagioclase (0.062 μm) versus pyroxene (0.107 μm). The plagioclase-hosted magnetite particles have a difference in shape, with median aspect ratios of 0.63 compared to 0.26 for pyroxene-hosted magnetite. The difference in shape is also observed for single-vortex (SV) particles, with median size and aspect ratio of 0.144 μm and 0.72, respectively, for plagioclase and 0.159 μm and 0.4, respectively, for pyroxene. The shape difference is not observed for multivortex (MV) particles, however, with median size and aspect ratio of 0.247 μm and 0.24, respectively, for plagioclase and 0.214 μm and 0.26, respectively, for pyroxene. There are three particles characterized as having swirl ("SW") domain structure (Appendix A3). This group refers to a distinct set of large grains with low aspect ratios (median size of 0.370 μm, and aspect ratio of 0.2 that represents a particle 0.500-μm long and 0.100-μm wide). These particles were observed only within pyroxene, and because of their large size, we were only able to complete 20 hysteresis runs for these three particles. We were unable to complete backfield runs due to insufficient computer memory. Figure 5 shows the micromagnetic states of representative particles and how they were categorized after minimizing their total magnetic energy from random initial moments. The domain state estimations represent one of many possible local energy minima. The minimization from random initial moments should ideally be repeated multiple times to categorize the full range domain states, and potentially identify the global energy minimum state. This is computationally expensive, but should be encouraged for future studies and as a part of the FIB-nt method in general. Insets show the corresponding average hysteresis loop for each particle. The arrows are color coded by their vertical (M z ) component where up = red and blue = down. The vortex core is represented by contour surfaces (in green) enclosing volumes with increased helicity of the magnetization field m defined by SD particles have uniform magnetic moments aligned throughout the particle, although the spins may bend along the demagnetizing field lines at edges and corners (flower state). The SV is a inhomogenous magnetization state with lowest total micromagnetic energy beyond the SD limit (Hubert & Schäfer, 1998). The net energy decrease is obtained by reducing the magnetic stray-field energy through flux closure at the expense of increasing exchange and anisotropy energy. Therefore, this transition can only occur when the decrease in demagnetizing energy overcompensates this increase. Visually, the SV core extends partially or fully through the particle. Rave et al. (1998a) and Hubert and Schäfer (1998) find that there is a transition from a simple vortex-state to a well-defined MD state, where vortices are embedded features within the soft-magnetic Bloch walls. In this paper this transition is described as a MV state from the zero-field minimization. The MV state is characterized by the presence of two or more well-defined vortex cores Green contour (vortex core) represents areas of increased helicity of the magnetization field. Hysteresis, size, and shape parameters for each individual particle are summarized in Table 1. Location of these 12 particles are marked with bold text in the subsequent  in magnetite grains with different sizes. They observed the formation of primitive bloch-like 71° domain walls that separate uniform regions of magnetization, similar to our observations in the three selected large particles in Appendix A3. In these particles, the magnetization pattern mainly consists of vortices or SWs that represent seeds of Bloch walls, but due to size constraints, cannot develop into conventional domain walls as they occur in well-developed MD patterns (Fabian et al., 1996;Hubert & Schäfer, 1998). Domain walls widths expected for MD magnetite are observed only in much larger particles (∼3 μm) and the SV structures are present until at least ∼1 μm (Nagy et al., 2019). Due to the lack of a distinctive terminology, the intermediate sized magnetite structures are defined as "SW" within this study.

Hysteresis Loops
Hysteresis results for all simulated inclusions, from both plagioclase (circles) and pyroxenes (triangles), are summarized in Figure 6, Figure 7 (Day plot), Figure 9 (Néel plot), Figure 11 (Butler-Banerjee plot), and Appendix A5. Each point represents the average hysteresis data acquired from 20 different field directions and is color coded according to the remanent domain state observed at zero field following energy minimization from a random starting configuration (purple = SD, yellow = SV, and teal = MV). Each symbol size is scaled according to particle size, defined as the diameter of a sphere with equivalent volume. The total number of particles adopting each remanent state is represented in Table 2.
The two red markers (one for each silicate) in Figures 7 and 9 represent the average hysteresis properties for all of the modeled particles. Assuming the two volumes studied provide a representative sample of the silicate-hosted component of magnetic mineralogy, these two points provide an estimate of the hysteresis properties that would be obtained if we measured the bulk properties of plagioclase and pyroxene extracts experimentally. These "bulk" hysteresis properties would normally be the only two points available to perform domain state diagnosis. However, when deconstructed to individual particles, we observe the data for both volumes extends across the entire range of the common diagnostic diagrams and all regions designated to specific domain states. The location of these bulk average points also suggest the domain states are unevenly distributed between the silicates. The bulk average for pyroxene is dominated by SD particles and for plagioclase is dominated by SV particles (Table 2). Figure 6 shows the distribution of particles between the domain states, host silicate, and M rs /M s ranges. For SD, there is a distinct separation, with tight clusters of M rs /M s ranges with pyroxene hosted particles at 0.5 and plagioclase hosted particles at 0.7. For SV particles, there are M rs /M s values where particles form clusters. Within plagioclase SV particles cluster at M rs /M s values ≈6.1 × 10 −2 , but there are a few points that extend the range up to M rs /M s = 0.7. SV particles in pyroxene show a bimodal distribution (red boxes) with clusters at M rs /M s values >0.3 and <0.3. MV particles show a similar tendency to form two separate clusters at M rs /M s values ≈ 0.1 and 0.4. This is most evident for particles hosted by pyroxene, but with only three particles extracted from the plagioclase is an insufficient sample size to indicate anything about specific M rs /M s ranges. What we have for now indicates that these particles cover the same M rs /M s ranges as MV magnetite in pyroxene. Three SW particles hosted by pyroxene did successfully complete minimization of the hysteresis loop. The spread in M rs /M s is from 0.12 to 3.44 × 10 −4 .

Day Plot
In Figure Figure 5 and Highlighted With Bold Text in Figures 7-11 The spread of points shows little regard for the particle size, highlighted by Figure 8, where particles with varying sizes and domain states coexist and even overlap. This observation is also evident throughout Figure 7. At M rs /M s values ≈0.5, the cluster of prolate particles is mainly hosted by pyroxene. SD particles hosted by plagioclase form a cluster at M rs /M s values ≈0.7.

Néel Plot
In a Néel plot (Figure 9), we observe a dominating SD "central ridge" that extends from 10 mT < H c < 120 mT. The SD particles mainly hosted by pyroxene at M rs /M s ≈ 0.5 align with the uniaxial SD (USD) theoretical line. Axis ratios show the majority of these particles are prolate in shape (Appendix A4) and correspond to the median aspect ratio cluster of 0.26 (Figure 4b). The feldspar-hosted SD particles (higher median aspect ratio of 0.62 in Figure 4b) form a cluster at M rs /M s values ≈0.7 and 10 mT < H c < 30 mT. These particles are oblate spheroids and the increase in H c for these particles represents a general increase in particle oblateness/flatness (Figure 9a). The low-coercivity (<10 mT) region is dominated by SV and MV particles. With decreasing coercivities below 10 mT, there is also a decrease in M rs /M s . Particles with M rs /M s values < 0.3 form a distinct cluster that is dominated by SV states (Figure 10). This cluster has a well-defined trend for particles from both plagiolase and pyroxene. While pyroxene-hosted particles are equally sized and evenly distributed in Figure 10, there is a decrease in particle size in plagioclase above H c ≈ 9 mT. Particles below 9 mT and M rs /M s = 0.07 align with the theoretical "USD + SP". Those above 9 mT demonstrate an increased upwards slope that follows the USD theoretical lines from Tauxe et al. (2002).

Butler-Banerjee Plot
The Butler-Banerjee plot in Figure 11 delineates a two-dimensional size-shape range for SD particles of magnetite and titanomagnetite (Butler & Banerjee, 1975). Contrary to the Day and Néel plots, which interpret domain state from hysteresis parameters, the Butler-Banerjee plot predicts domain states based on the length of a particle's long axis compared to its axis ratio. One is therefore able to characterize a particle's domain state purely based on shape and size. The theoretical boundaries in the background were calculated either by a domain model (Butler & Banerjee, 1975;dashed line) or by 3D micromagnetic modeling (Muxworthy & Williams, 2008;solid lines). In Butler and Banerjee (1975), the particles have cuboidal shapes with 180° domain walls, while in Muxworthy and Williams (2008), the magnetic particles are SD or SV parallelepipeds arranged in chains to assess magnetostatic stabilization in chains of magnetosomes from magnetotactic bacteria. Our data set defines a clear transition between SD and SV states that is in good agreement with the "no interaction" boundary identified by Muxworthy and Williams (2008) for isolated parallelepipeds. Our data in Figure 11 clearly reflect the influence of the silicate host mineral on the inclusion shape and subsequently its magnetic properties. While plagioclase-hosted magnetites (circles) plot predominantly at aspect ratios (AR) above 0.6, pyroxene-hosted magnetites (triangles) have AR less than 0.6. Particles that are able to adopt MV remanence states typically plot in the region with AR <0.4 and long axis >0.300 μm. The phase boundary between SD and SV/MV states in this region deviates from that of Muxworthy and Williams (2008), most likely reflecting the difference between the real particle geometries used here versus the ideal parallelepipeds used by Muxworthy and Williams (2008). Figure 11 demonstrates that MV and SW particles start to dominate over SV states with decreasing AR and increasing long axis length.

Discussion
As micromagnetics takes us beyond the SD theories that have been the foundation of rock magnetism since the pioneering work of Néel, the ability to obtain an accurate breakdown of the different classes of magnetic behavior that exist within a sample becomes increasingly important. One of the most striking results of this study is the shear range of magnetic domain states and hysteresis properties that are observed within just two small ∼2,000-600 μm 3 volumes of silicate. The observed behavior spans the entire range of popular domain-state diagnostic plots (Figures 7, 9, and 11). In contrast, bulk hysteresis measurements reduce this complexity to a single average data point (e.g., the red points in Figures 7, 9), discarding valuable information about the true nature of the underlying magnetic ensemble. Classifying a bulk sample as either "SD," "PSD," or "MD" on the basis of the Day diagram is clearly meaningless when particles in the SV state plot across the entire "SD/PSD/MD" range, and the largest MV particles plot much closer to the SD region than the MD region (Figure 7). Of the three diagnostic plots explored here, the most successful is the Butler-Banerjee plot ( Figure 11). This diagram works well because it recognizes the fundamental importance of particle size and shape as the dominant factors controlling the domain state of magnetite. However, our results also demonstrate that size and shape alone are often poor predictors of magnetic hysteresis properties. Generalized geometric information, such as equivalent spherical volume, is not sufficient to predict the magnetic properties-the specifics of individual particle geometries really do matter. Although time consuming and computationally intensive, the FIB-nt-FEM workflow presented here provides a route to rock magnetic characterization that not only yields a quantitative breakdown of the domain states present, but a realistic estimate of the range of magnetic properties associated with the particle ensemble. As the database of magnetic properties linked to specific particle geometries expands, it will eventually become possible to use look-up tables and/or machine-learning approaches to translate 3D particle geometries into useful magnetic parameters (including, e.g., energy barriers for thermal switching) without the need to perform ad-  Dunlop (2002) and Parry (1982). Regions for stable single-domain, pseudo-single-domain, and multidomain from Dunlop (2002). Bold text indicates the positions of particles from Figure 5 and Table 1. ditional micromagnetic simulations. With increased automation, the experimental methods outlined here could indeed become a routine feature of future rock magnetic characterization workflows.
The Néel plot (Figure 9) shows a distinct "central ridge" of SD particles along the USD line at M rs /M s ≈ 0.5. Tauxe et al. (2002) argue that the H c increase for SD particles is proportional to an increase in the length to width (L/W) ratio. Figure 4a in Tauxe et al. (2002) shows that the increase from H c = 30 to 70 mT for USD particles is caused by an increase in L/W from 1.3:1 to 2:1 While we do observe that particles with H c >100 mT have high L/W ratio (>10:1), for coercivities <100 mT, we do not find a clear correlation between coercivity and aspect ratio, which may be due to the large spread of USD particle sizes (0.022-0.253 μm). Our data imply that the size of a USD particle is at least as important as its aspect ratio in determining the switching field, highlighting again the need for detailed geometrical information to produce an accurate prediction of the magnetic properties of an ensemble. Another interesting modeling result is that the characteristic coercivity range for SD magnetite particles predominantly covers the range 10 mT < H c < 100 mT, which is in line with experimental results . Without external stress, higher H c values appear to require very unusual particle geometry, while SD structures with H c < 10 mT can occur only in a very narrow grain size and shape region.
A cluster of feldspar-hosted SD particles with M rs /M s values of 0.7 (Figures 8 and 9) falls in between the theoretical values for cubic SD (M rs /M s = 0.83-0.87; Joffe & Heuberger, 1974;Tauxe et al., 2002) and uniaxial SD (M rs /M s = 0.5). This remanence ratio is similar to the value of M rs /M s = 0.75 calculated by Harrison et al. (2019) for particles with a combination of uniaxial anisotropy (restricting moments to a basal plane) and hexagonal anisotropy (defining multiple easy axes within that basal plane). This cluster contains mainly oblate particles (Appendix A4), and therefore the intermediate value of M rs /M s is entirely consistent with these particles having a mixed uniaxial/multiaxial anisotropy. A positive correlation between coercivity and oblateness is observed (inset to Figure 9). This improved correlation may be a result of the smaller size variations within this oblate cluster (0.023-0.090 μm) compared to that of the pyroxene-hosted USD prolate cluster.
The SV domain state is observed across a wide range of particle sizes and shapes, and produces correspond- Another key result of this study is the observation of SV and MV particles that display hysteresis properties similar to those expected for much smaller SD particles (Table 3). The MV state is stable in large grains (median 0.225 μm) with low aspect ratios <0.5 (Figures 4a and 4b). Despite their larger sizes, these particles plot more toward the upper left of a Day diagram (SD and PSD regions), rather than toward the MD region. Large, slightly flattened needle (a ≫ b > c) particles with low aspect ratio readily adopt MV states during zero-field relaxation from a random starting configuration (e.g., PYR-23 in Figure 5), but are forced into a uniformly magnetized remanence state after exposure to a saturating field. This produces an apparent disconnect between the domain-state classification used in this study and the observed hysteresis properties. This disconnect becomes significant if, as intended, the zero-field relaxation state is more representative of a weak-field thermoremanent magnetization (TRM) state than the corresponding saturation remanence state. Such disconnects have been observed in dusty olivine by Lappe et al. (2013), where stable SV states were observed after acquisition of weak-field TRM but were converted to metastable SD states after application of a room-temperature saturating field. Our study raises the question of whether high-field hysteresis measurements are always a useful means to classify the potential weak-field remanence states adopted by a given particle geometry, especially for SV and MV states, which may well be the dominant carriers of remanent moment. Such ambiguities are avoided using the FIB-nt-FEM approach, ensuring that NIKOLAISEN ET AL. 14 of 20 10.1029/2020GC009389  Table 1. Theoretical lines of USD, USD + SP and CSD from Tauxe et al. (2002). Inset (a) shows the increasing oblateness for the SD particles at M rs /M s = 0.7 with increasing H c . SD, stable single-domain; SP, superparamagnetic; USD, uniaxial SD an appropriate combination of physical models is always used for a given ensemble. A full theory of the remanence state adopted by such particles requires us to calculate the multitude of energy barriers between all possible remanent states. Although this is an active area of research, it is outside the scope of the current study. As discussed above, the state adopted will be highly dependent on the external field history. By introducing the energy barrier calculations, we would expect to see a better separation line between the remanence states of all the particles, resulting in adjustments to the theoretical lines in Figure 11 from Butler and Banerjee (1975) and Muxworthy and Williams (2008).
The blue particles seen in Figures 9-11 are large (0.370-μm median equivalent sphere diameter) oblate to prolate shapes (Figures 4a and 4b and Appendix A4) which show a tendency to form domain walls when minimized in zero-field (Appendix A3). The backfield curves were not able to be calculated on a reasonable timescale, and therefore are not plotted on a Day diagram (Figure 7). On a Néel diagram (Figure 9), two SW particles plot close to the origin, as would be expected for MD behavior. One, however, displays similar properties to MV particles, suggesting that they share many characteristics, and that the transition from MV to SW to MD is likely gradual one for particles with this shape.
In Figure 7, the theoretical Day plot SD-MD and SD-SP mixing lines show the modified calculated binary curves Dunlop (2002) and the SD-MD mixing curve derived from natural magnetite mixtures in Parry (1982). Because effects like thermal variations, stress, interactions, impurities, or inclusions are not included in the models of this study, the hysteresis properties of our modeled particles represent ideal magnetite crystals. When the entire data set is compared to theoretical mixing lines, the closest agreement occurs for the SD-MD mixing line from Parry (1982), that represents a grain-size trend of measured mixtures of SD and MD particles. The samples used by Parry (1982) were generated by crushing natural MD magnetite and sieving for various grain sizes. Because crushing introduces internal stress anisotropy that may be larger than the magnetocrystalline anisotropy, the correlation with our data set appears to be spurious, but based on the fact that our data consistently lie above the binary mixing lines from Day et al. (1977) and Dunlop (2002).
It is unavoidable that the magnetic models presented here contain several uncertainties related to reconstruction artifacts, model assumptions, and unknown physical influence factors. Shape uncertainty is directly connected to the extraction of selected mineralogy from the binary stack and the accuracy of the  Table 1, and the bulk average for plagioclase. Particle PYR-07 is displayed in Appendix A3.
resulting finite element meshes. The alignment of fiducial marks in the image stack to reconstruct the volume contributes to the shape uncertainty. The carbon rods are deposited from carbon gas reacting with the ion beam. As observed in the reassembled stack, there are irregular shifts in image placement that also impact the reconstructed shape, especially of the smallest particles. Evidence for this is seen in the meshes as irregular pinching in elongated particles, where images are shifted. Thresholding in FIJI is another uncertainty in reconstructing the true shape of particles. With the variation in size, the extraction of the smaller particles will become suppressed when selecting the threshold based on the median size. A simple calculation of the volume change when increasing the particle volume by adding a pixel in each direction leads to an estimate of up to 10%-25% for the larger particles and 50%-130% for the smaller particles (based on an increase in the representing sphere diameter of 7-nm pixel width). The true change in volume is dependent on more complicated measures such as the particle orientation and the milling thickness. We were also concerned that variable amount of mesh smoothing would possibly change the hysteresis parameters and therefore the domain state. To test this, we calculated average hysteresis properties for three particles of different sizes at three different stages of smoothness. First stage is where the particle shape is strongly voxelated (directly from the binary image), second is at intermediate smoothing, and the last where all sides are completely rounded (axis ratio is preserved in all stages). From the resulting hysteresis properties, we did not observe any significant changes beyond the deviation expected for models starting with randomized moments. This is in agreement with a study by Rave et al. (1998b) Figure 11. Butler-Banerjee diagram that shows the relationship between a particle's long axis and its axis ratio. Axis ratio = 0 references an infinitely long needle and axis ratio = 1 references an equidimensional particl (e.g., cube or sphere). Bold text indicates the position of particles visible in Figure 5 and Table 1, except for particle PYR-07 which is displayed in Appendix A3.
resolutions below the exchange length, particle surface discretization does not influence micromagnetic model results. Model assumptions that contribute to the overall uncertainty are the material parameters for magnetite that in the natural minerals can vary for example due to impurities. Physical unknowns are lattice orientations, related to the direction of the easy anisotropy axes, or internal stress (Hodych, 1990;ter Maat et al., 2020).
This study presents a workflow for extracting the 3D geometry of a nanoscale magnetic particle ensemble and calculating its micromagnetic properties. To date, it is the most comprehensive micromagnetic data set of natural remanence carriers and acts as a proof-of-concept study for obtaining statistical data on the magnetic domain state and hysteresis properties of silicate-hosted magnetic inclusions. This work is of particular relevence to single-crystal paleomagnetic studies, which offer the best possibility of extracting reliable primary paleomagnetic remanence data from ancient rocks with complex geological histories. The key challenge in such studies is targeting primary remanence carriers that are most likely to have escaped thermal, chemical, or viscous remagnetization. In the case of the Harcus intrusion, stable remanence carriers hosted in primary feldspar and pyroxene grains could be targeted directly, avoiding both remanence carriers hosted by secondary amphibole as well as interstitial MD magnetite that could carry a viscous overprint. The combination of FIB-nt and micromagnetic simulations provides a full statistical breakdown of the domain states present in a representative single crystal. Further studies will enable the range of local energy minimum states adopted by the SD, SV, MV, SW, and MD states, as well as the energy barriers that separate them, to be calculated, and the corresponding blocking/unblocking behavior during natural and laboratory cooling/heating cycles to be modeled. With that we would be one step closer to a comprehensive model of remanence acquisition that takes us beyond Néel's SD theory. Ultimately this will increase both the reliability of single-crystal paleomagnetic studies and our confidence to interpret them, paving the way to unlocking hitherto inaccessible parts of the geomagnetic record.

Conclusions
Starting from thin sections of geological samples, selected silicates were milled with SEM-FIB slice-andview. This resulted in two successfully reconstructed volumes of plagioclase (2,366.14 μm 3 ) and pyroxene (647.35 μm 3 ). Thresholding of binary image stacks extracted a total of 413 particles of magnetite measuring between 0.022 and 0.420 μm, which were then individually meshed. Separate micromagnetic hysteresis modeling was performed on each mesh in MERRILL with 20 different external field direction. Magnetic hysteresis properties for the total data set showed M rs /M s span from 8 × 10 −4 to 0.7, H cr /H c from 1.14 to 127 and H c from 3 × 10 −4 to 107.7 mT. Because the geometries of the individual magnetites are exactly known, the results were used to probe three of the commonly used domain-state diagnostic plots. Our results highlight the following points: 1. FIB-nt protocols have been developed for 3D imaging of silicate-hosted magnetite inclusions down to the stable SD limit of 30 nm 2. Silicate-hosted magnetite inclusions in plagioclase and pyroxene span a wide range of domain state behaviors, including SD, SV, MV, SW, and MD 3. In plagioclase, a dominant cluster of oblate SD particles is found with mixed uniaxial/multiaxial anisotropy and M rs /M s = 0.7. In pyroxene, SD particles are predominantly prolate with uniaxial anisotropy and M rs /M s = 0.5 4. The range of hysteresis properties observed in just two ∼650 − 2300 μm 3 volumes of pyroxene and plagioclase, span almost the entire range of common domain-state diagnostic plots. This indicates that bulk average hysteresis parameters cannot predict a unique domain state, even for the particle ensemble inside a single-exsolved silicate crystal 5. Because shape and size equally determine the hysteresis properties of natural magnetite inclusions, the Butler-Banerjee plot performs best in terms of domain-state diagnosis, despite failing to account for MV states 6. The morphology of the magnetite inclusions highlights the need to extend Néel's SD theory in the context of single-crystal paleomagnetism toward SV and MV magnetization states, as these may be dominant carriers of stable paleomagnetic remanence in silicates