Multi-scale three-dimensional characterization of iron particles in dusty olivine: Implications for paleomagnetism of chondritic meteorites

Sample

The sample is a grain of dusty olivine extracted from a chondrule of the Semarkona LL3.0 chondrite (sample DOC5 from Fu et al. 2014). A full description of all sample preparation steps performed prior to our study is given by Fu et al. (2014). The sample was mounted on a quartz disk stub with epoxy, and the magnetic reference axes of the grain were marked. This enables the extraction of nanoscale information to be correlated back to the macroscopic magnetic measurements made by Fu et al. (2014). After optical investigation, the quartz disk was mounted onto a 38 mm SEM stub and carbon coated. Much of the sample surface was initially obscured by epoxy (Fig. 1a). Small bright specks visible in the exposed region of the grain correspond to the metallic Fe particles that define dusty olivine. FIB-nT was performed on the region highlighted in Figure 1a. The chosen region was adjacent to the region where a TEM lamella had previously been extracted from the olivine grain (Fu et al. 2014). The same TEM lamella was used to perform the additional TEM measurements reported here. The imaging surfaces for both FIB-nT and TEM studies were acquired parallel to each other, allowing us to make direct comparisons across the length scales captured.

Figure 1

(a) SEM micrograph of Semarkona Dusty Olivine Grain DOC5. White box indicates region used for FIB-nT and S/TEM studies, (b) Top down FIB micrograph of FIB-nT MROI after trenches have been milled. Hourglass structure is a fiducial mark to assist in mill pattern placement. (c-e) The cyclic sequence of steps for each 20 nm slice in FIB-nT volume collected.

FIB-nT

FIB-nT was performed using a FEI Helios Nanolab Dual Beam microscope at the University of Cambridge. The sequential slicing and imaging sequence was controlled by the Auto Slice and View G3 (ASVG3) application. All FIB milling was performed using an accelerating voltage of 30 kV. The MROI was prepared by depositing a protective 10 μm by 15 μm by 1 μm tungsten pad using ion beam induced deposition with an ion beam current of 3 nA. The MROI was isolated from the bulk sample by selectively milling 20 μm deep trenches on three sides of the region defined by the W pad (Figs. 1b-d) using a 7 nA ion beam current. The front trench allows full viewing access to the cross-sectional surface and the side trenches minimize re-deposition effects associated with the sequential sectioning process. Figure 1b shows the tomographic region of interest immediately after clearing away excess material. Finally, a pair of fiducial marks was created before starting the automated sequence, using the 3 nA ion aperture for W deposition and 300 pA aperture for feature milling. Figure 1c is a FIB micrograph used for slice placement, showing the actual field of view used, taken at an arbitrary point in the automated slice-and-view routine. The fiducial mark seen in Figures 1b and 1c controls the placement of each slice of the tomographic sequence, whereas the fiducial mark seen in Figure 1d minimizes the amount of image drift in the SEM image stack. Each 20 nm thick tomographic slice was milled away using a 920 pA ion beam current. All milling was performed at 52° stage tilt, which is normal to the FIB.

Each image of the tomographic sequence was recorded using a dwell time of 20 μs with an 8 bit grayscale and a scan area of 1024 by 884 pixels. The horizontal field width of the final image was set to 10 μm. This gave a final pixel size of 9.766 nm. Imaging of the cross-sectional cut face was achieved using backscattered electron (BSE) imaging with the SEM operating in immersion mode at a low accelerating voltage of 2 kV with a beam current of 86 pA. We used the through the lens detector in backscattered electron mode for the strong material contrast mechanism between the olivine host matrix (dark gray) and the Fe nanoparticles (bright grays) (Fig. 2a). This combination of through lens detection and backscatter mode also eliminates contrast artifacts like shadowing and allows for higher spatial resolution over Everhart-Thornley detector geometries. SEM imaging conditions were further optimized to reduce noise in the BSE image as well as to ensure that only the cross-sectioned surface was imaged. Reduction of the electron interaction volume was achieved by using the low accelerating voltage of 2 kV. Detector noise in the image was minimized by reducing the working distance from 4 to 3.3 mm. Additionally, to protect the pole piece at the shorter working distance the sample was tilted to 47.8°. Movement between the SEM imaging position and the FIB patterning position was controlled automatically through ASVG3. Using this automation routine, 262 slices were recorded over 15.3 h, with a specified slice thickness of 20 nm resulting in a total milled thickness of 5.2 μm.

Image post-processing

Data post-processing and analysis were performed using a combination of the commercial package ORS Visual SI and the open-source platform FIJI, based on the ImageJ image analysis distribution (Schindelin et al. 2012). The 262 images were loaded as a single three-dimensional image stack into ORS Visual SI. The pixel aspect ratio was adjusted to correct for the sample tilt and the slice thickness of 20 nm was applied to the image stack. This resulted in voxel dimensions of 9.8 nm by 13.2 nm by 20 nm. Although the fiducial mark seen in Figure 1b helps to minimize image stack drift, erosion of this feature during the sequential milling of the MROI resulted in a progressive drift in the final image position on the cut face. Applying the Normalized Mutual Information alignment algorithm within ORS Visual SI the image stack alignment tool allowed for precise and jitter-free alignment of the reconstructed volume.

The image stack was segmented to define the Fe particles and olivine matrix (Fig. 2a). Several approaches for this have been documented in the literature (Holzer et al. 2004; Bushby et al. 2011; Keller et al. 2011a). We found that the best results were obtained using simple grayscale thresholding followed by noise reduction and manual artifact removal. Middle gray levels (114-212) associated with Fe particles were assigned a saturating value of 255. All other pixels were set to 0. After applying this binary segmentation, residual detector noise was eliminated by one pass of the FIJI/ImageJ despeckle filter. Next, residual W deposition, milling curtain artifacts and particles partially lying on the edges of the reconstructed volume were manually removed from the image stack. Close inspection of the original images revealed saturated black and saturated white pixels around some Fe particles. It is thought that these are pockets of silica glass formed as a by-product of the solid-state reduction reaction that created the Fe particles (Leroux et al. 2003). As these regions are non-conductive, they tend to induce charging effects in the images registering as saturated white. Related to these are saturated black voxels that upon inspection of the original image stack can be seen to be either non-charging glass (which charge up in subsequent images) or holes in the olivine crystal. The holes and charging artifacts can be seen in Figure 2a. Grayscale binary segmentation does not fully isolate the Fe particles from these parasitic features and so these regions were manually removed from the segmented images. Figure 2b shows the same image after binary segmentation and manual artifact removal.

Figure 2

(a) BSE SEM micrograph from slice 151 where medium gray is olivine and the light gray are the Fe nanoparticles. Also visible are non-crystalline materials such as charging of amorphous silica (saturated white) and holes and/or uncharged silica (saturated black). (b) Binary segmented image of same BSE image after denoising and manual artifact removal this leaves particles as white voxels and all other materials as black.

Quantitative three-dimensional analysis of the resulting scaled aligned and segmented image stack was performed using the ImageJ/FIJI plug-in BoneJ version 1.3.15 (Doube et al. 2010; Carriero et al. 2014). The BoneJ Particle Analyzer generates a surface mesh for each particle and determines a best-fitting ellipsoid to that surface. Errors in the ellipsoid fitting are of the order of the voxel size and therefore can become significant for particles that are defined by a small number of voxels. Additionally non-physical results were obtained for groups of voxels that appeared only on one image of the stack. These voxel groups consisted of either segmentation noise that was not fully removed, or very small particles with radii less than 20 nm in one dimension (i.e., these are particles that are less than the thickness of a tomography slice). This establishes our 3D resolution limits and means we do not observe any particles with a diameter of less than 40 nm. However, TEM studies did not reveal any particles smaller than this. After artifact removal we were left with 246 particles in the 710 μm³ volume analyzed by FIB-nT.

TEM

To complement the mesoscale tomographic information provided by FIB-nT, we also performed high-resolution TEM studies of individual particles using scanning transmission electron microscopy (STEM). The TEM lamella was fabricated using FIB milling (Fu et al. 2014). All studies were performed at the Technical University of Denmark (DTU) on an FEI Titan 80-300 TEM equipped with a field-emission electron source, monochromator, spherical aberration probe corrector, Lorentz lens, and electrostatic biprism. All measurements were made using 300 kV accelerating voltage. Preliminary Lorentz microscopy and electron holography observations were reported by Fu et al. (2014).

Here we present bright-field TEM and STEM imaging, electron diffraction data as well as dark-field STEM tomography results. Both TEM and bright-field STEM modes allow us to image crystallographic features such as dislocations and sub-grain boundaries (Williams and Carter 2002; Crewe and Nellist 2009). Dark-field STEM tomography of a single particle was performed using a camera length of 130 mm to produce a strong material contrast between the Fe particles and the olivine crystal. The tomographic series was collected at a magnification of 28 500x (giving a pixel size of 3.26 nm) using the high-angle annular dark-field (HAADF) detector. The tilt series consisted of an image taken every 2° for tilts from −76° to +76°. Alignment and reconstruction of the tilt series was achieved using Inspect 3D with the SIRT algorithm (Gilbert 1972). Visualization was performed using Avizo Fire.

Micromagnetic modeling

A selection of Fe metal particles, representing the range of sizes and shapes within the ensemble, were chosen to perform detailed micromagnetic simulations. Each particle was cropped from the segmented FIB-nT stack and converted to a tetrahedral finite-element mesh in a multi-step process. An example of the initial geometry of a particle defined by the FIB-nT is shown in Figure 3a. Each rectangular block represents a single 9.8 × 13.2 × 20 nm voxel. This representation of the particle was used to generate a bounding polyhedron that best approximates the actual particle surface, where each point of the resulting mesh must solve a Poisson boundary condition (Fig. 3b). We further refine the surface mesh by passing it through a surface smoothing routine (Fig. 3c). We coarsen the smoothed surface using a Delaunay triangulation routine to produce a surface mesh at the desired resolution of 5 nm (Fig. 3d). The final triangular surface mesh was imported into the software package CuBit (KitWare), where it was turned into a tetrahedral volume mesh. We used an initial surface mesh with average node spacing 5 nm, which was then used to generate tetrahedral nodes on average every 5 nm throughout the volume. Although this resulted in a mesh size slightly bigger than the 3.4 nm exchange length for iron, it enabled the number of elements in the model to be kept below ~300 000 and provides acceptable resolution for modeling simple vortex micromagnetic structures.

Figure 3

Transformation of tomographic particle voxels into tetrahedral meshes for FEM modeling for particle number 8, the smallest prolate particle captured in the FIB-nT volume. (a) Initially the individual particles form a rough box-like volume, described by the voxels of the original FIB-nT volume. (b) By fitting a bounding polyhedron subject to the Poisson equation, a rough triangular surface mesh is generated. (c) Applying a smoothing filter to the Poisson surface results in a mesh with triangular nodes spaced every <2.5 nm. This is too high a mesh density for efficient micromagnetic modeling. (c) Using a Delaunay triangulation function, the surface mesh is coarsened to the desired 5 nm resolution triangular surface mesh. This surface is then used to seed the interior tetrahedral elements of the same size.

Micromagnetic modeling was performed using Micromagnetic Earth Related Rapid Interpreted Language Laboratory (MERRILL), a micromagnetics package optimized for the rock magnetic community developed by K. Fabian and W. Williams (Williams and Fabian, personal correspondence 2016). MERRILL uses a Finite Element Method/Boundary Element Method (FEM-BEM) to solve for the magnetic scalar potential inside the particle and thereby calculate the demagnetizing energy of the system. The use of FEM-BEM avoids the need to discretize the non-magnetic volume outside the particle. Simulations were performed by minimizing the total micromagnetic energy. This consists of summing the exchange, cubic anisotropy, magnetostatic, and demagnetizing energies. Energy minimization was performed using a conjugate gradient method, specially adapted to micromagnetic problems. MERRILL has been successfully tested against mMAG Standard Problem 3 (http://www.ctcms.nist.gov/~rdm/mumag.org.html).

Material parameters used were appropriate for pure iron at room temperature: saturation magnetization M_s = 1715 kA/m, exchange constant A = 2 × 10^-11 J/m, and cubic anisotropy with K₁ = 48 kJ/m³ (Muxworthy and Williams 2015). We arbitrarily set the cubic <100> axes parallel to the X, Y, and Z axes of the volume reconstruction. We will show that this does not greatly influence our analysis, as the magnetic behavior of the particles studied is dominated by shape rather than magnetocrystalline anisotropy.

Using the method of Hubert (1967), it can be estimated that magnetostrictive effects become important for iron only when the magnetostrictive energy density 9/2 (c₁₁-c₁₂) λ₁₀₀² is of similar size as the energy density ~2(A K₁)^1/2/d generated by a 180° domain wall, where d is the dimension of the particle (Hubert 1967; Fabian et al. 1996; Hubert and Schäfer 1998). For iron with elastic constants c₁₁ = 241 GPa, c₁₂ = 146 GPa (Lee 1955), and λ₁₀₀ = 22 × 10^-6 (Radeloff 1964), this occurs only for d > 9 mm such that magnetostriction can safely be neglected for all modeled particles. Also magnetoelastic interaction is neglected because it is assumed that any strain involved in the formation of the particles has been relaxed by plastic deformation such that no noticeable internal stress field is present.

Simulations were performed using an Apple iMac with a 3.4 GHz Intel i7 processor and 24 GB of RAM. Each particle was initialized with uniform magnetization along either the X, Y, and Z axes of the reconstructed volume. Fields varying from 1000 mT to −1000 mT in steps of 10 mT were applied along X, Y, and Z. The converged set of magnetic moments obtained after each field step was subjected to small random rotations (maximum angle 20°) and then used as the basis for the starting condition for the next field step. This step is to ensure that the energy minimization is not trapped in a local energy minima. The average magnetization projected on each of the X, Y, and Z axes were calculated at each step to generate the upper branch of the hysteresis loop. Lower branches were not calculated directly using micromagnetics, but are presented for visualization purposes under the assumption that these are symmetrically equivalent to the upper branch.

FIB-nT

The reconstructed dusty olivine volume is shown in Figures 4a-4c. A qualitative analysis (see movie in supplemental^[1] information) reveals that: (1) the particles are loosely arranged in planar sheets, (2) the particles tend to be flattened in the direction perpendicular to the sheets, and (3) there is a preferred orientation of particle elongation along a direction within the sheets. We will present the crystallographic analysis below when we present the TEM and STEM results. Particles are widely distributed in terms of their size and aspect ratio. Figure 5 summarizes this distribution by plotting a histogram of the best-fit ellipsoid diameters for the major (Fig. 5a), intermediate (Fig. 5b), and minor (Fig. 5c) axes. The average particle radii are 242 ± 94 nm by 199 ± 80 nm by 123 ± 58 nm. To classify the particles in terms of their tendency toward either uniaxial prolate or uniaxial oblate symmetry, we plot the aspect ratio of the major to intermediate radii against the aspect ratio of the intermediate to minor radii in Figure 6a (Flinn 1962). This “Flinn” plot also scales the size of each data point to the size of the major radius. The line y = x separates prolate particles (above the line) from oblate particles (below the line). There are significant populations of small, flattened particles that plot close to the horizontal (uniaxial oblate) and elongated particles that plot close to the vertical (uniaxial prolate) axes. As seen qualitatively in the Flinn plot and quantitatively in the histogram in Figure 6b, 95% of the particles are classified as oblate in aspect ratio with a mean Flinn ratio of 0.76. This gives the majority of particles a triaxial symmetry. This non-spherical and non-uniaxial aspect ratio has profound implications for the magnetic anisotropy of each particle, which we will explore in more detail using FEM models. Of the reconstructed particles, only 11 possess a prolate aspect ratio. The three smallest of these are close to the limitations of our reconstruction resolution. For these particles, the radii are between 2 to 6 voxels long in any one direction, which means that the uncertainty in the ellipsoid-fitting algorithm can be on the order of the size of the particle. We found that extracting these individual particles from the larger tomographic volume and rerunning the BoneJ analysis lead to small changes in the fitted ellipsoid parameters, causing them to be reclassified as oblate. Reported population statistics refer to the results of the whole ensemble fitting, which are accurate for all but these very small particles.

Figure 4

Orthogonal views of the reconstructed volume on the (a) X-Y plane, (b) X-Z plane, and (c) Y-Z plane. (d) Close-up detail of several larger particles where faceting is observable.

Figure 5

Histograms of best-fitting ellipsoid diameters for the (a) major, (b) intermediate, and (c) minor axes of each Fe particle in the FIB-nT reconstruction.

Figure 6. (a)

Flinn plot of the particle aspect ratios. The 45° line describes particles with a Flinn ratio of 1. The particle aspect ratios are plotted with color and point size scaling as a function of the major radius length in nanometers. (b) Histogram of Flinn ratios (<1 for oblate, >1 for prolate) from the FIB-nT volume.

The orientations of major, intermediate, and minor ellipsoid axes are shown in Figure 7. All three axes show pronounced clustering of their principal axes, confirming the qualitative assessment above (note that the ellipsoid fitting does not distinguish between positive and negative vectors). The minor ellipsoid axes (blue circles) are tightly clustered in the direction normal to the sheets, whereas there is a much broader spread of major and intermediate axes within the sheets (red and green circles, respectively). Solid triangles in Figure 7 show the principal axes of the anisotropy of susceptibility of anhysteretic remanent magnetization (ARM) for the whole dusty olivine chondrule, as determined by scanning SQUID microscopy (Fu et al. 2014). The coordinate system used for reporting results is based around the orientations of the FIB-SEM microscope. We have transformed the coordinates Fu et al. (2014) to agree with the microscope system. The average orientation of the major ellipsoid axes coincides with the direction of highest ARM susceptibility, whereas the intermediate and minor susceptibilities lie at angles of ~40° from the average orientations of the intermediate and minor ellipsoid axes, respectively.

Figure 7

Stereographic projection showing the orientation of major (red circles), intermediate (green circles), and minor (blue circles) axes of the best-fitting ellipsoids. Open circles are upper hemisphere, closed circles are lower hemisphere. The three solid triangles plot the measured the anisotropy of ARM susceptibility for the entire dusty olivine grain studied (Maximum susceptibility is red, intermediate is green, and minor is blue).

Tem and Stem analysis

Figure 8a shows a bright-field STEM mosaic of the entire lamella. The prominent dark feature running horizontally along the top of the image is the remains of the Pt capping layer that was deposited on the surface of the sample to protect it during FIB milling. The Pt layer became partially detached from the surface after the lamella was plasma cleaned prior to insertion in the TEM. The Fe particles appear dark against the mid gray olivine background. Silica glass regions appear as light gray blebs parasitic to the Fe particles. Electron diffraction patterns were obtained for the olivine and Fe particles inside the red-boxed region shown in Figure 8a. Figure 8b has been oriented with respect to the region highlighted in larger montage image. The upper right-hand diffraction pattern in Figure 8b was collected from olivine in this region with an α-tilt of −3.5° and a β-tilt of −1.1°, and corresponds to a [130] zone axis (~9° from [010]). The diffraction pattern in the lower left corner of Figure 8b was obtained from the smaller of the two Fe particles with an α-tilt of 8.8° and a β-tilt of −6.8°, and corresponds to a [001] zone axis. The short axis of the particle corresponds to the [100] direction of Fe and lies normal to the (001) plane of olivine. This is in agreement with the orientation data collected using the FIB-nT (Fig. 8, stereonet inset upper left), and highlights how the underlying texture of the dusty olivine is crystallographically constrained by the olivine host.

Figure 8. (a)

STEM bright-field montage showing Fe particles (dark), olivine (gray), and silica blebs (bright). Sub-grains are visible as small changes in grayscale in the olivine host. (b) Red inset: upper left panel shows bright-field TEM of the two highlighted Fe nanoparticles; upper right panel shows the diffraction pattern of olivine matrix collected at an a-tilt of −3.5° and a b-tilt of −1.1°; lower left panel shows the diffraction pattern for the smaller Fe particle; lower right panel shows the STEM tomographic reconstruction of the smaller Fe nanoparticle. (c) Blue inset shows a bright-field TEM image of the highlighted Fe nanoparticle. Faceting parallel to sub-grain boundaries, defined by multiple dislocations, can be observed. Stereonet inset upper left reproduces FIB-nT orientation data rotated to be aligned with STEM bright-field image.

The lower right-hand image in Figure 8b is a 3D visualization of the small Fe particle obtained using STEM tomography (movie of particle in supplemental¹ materials). The voxel size is 3.3 nm, at least a factor of 3 smaller than the voxel size obtained via FIB-nT. The particle dimensions are 232 × 205 × 232 nm (X, Y, Z with respect to the image plane), giving it a slightly oblate profile. These dimensions mean that the particle is one of the smaller ones in the population. It is representative of the largest complete particle that could be imaged using STEM tomography (larger particles were truncated by the surfaces of the TEM foil). Particles of this size (and smaller) are observable using the FIB-nT approach (albeit at considerably lower spatial resolution), so there is some overlap between the size ranges accessible by the two techniques. 2D particle analysis of this lamella does not reveal any particles smaller than those observed in the FIB-nT volume.

Figure 8a shows sub-grain boundaries in the olivine running parallel and perpendicular to [001]. As the diffraction pattern for the olivine in Figure 8b is oriented with respect to Figure 8a, we see that the olivine’s [001] is perpendicular to the sharp sub-grain boundaries. In contrast, the sub-grain boundaries running parallel to the [001] are slightly blurred. This broadening we interpret to be the result of the trace of the (010) intersecting the 100 nm thick lamellae surface with an angle of around 9° as noted above. These observations about the arrangement of the sub-grain boundaries are in line with the previous study by Kirby and Wegner (1978), which demonstrated that the dislocation arrays in olivine concentrate along the {100} lattice planes. The dislocations defining the sub-grain boundaries are more clearly visible in the bright-field TEM image (Fig. 8c), which was taken from the region outlined in blue. The sub-grain boundary seen on the right of this image is parallel to (001) of the olivine and lies parallel to a prominent (100) facet of the Fe crystal. This subgrain boundary is also parallel to the plane containing the major and intermediate axes of the ensemble (inset upper left). Previous studies of natural and synthetic dusty olivines (Leroux et al. 2003; Lappe et al. 2013) have suggested that the Fe nanoparticles arrange along dislocation arrays associated with the sub-grain boundaries. By measuring the crystallographic information of the lamella from the same region, we are able to demonstrate that the particle ensemble does indeed arrange in sheets related to crystallographic planes of the olivine host crystal.

Micromagnetic simulations

Remanence states and magnetic moments

The results of micromagnetic simulations for nine selected particles (eight particles extracted from the FIB-nT stack plus the STEM tomography particle shown in the red inset of Fig. 8) are summarized in Table 1. A range of oblate (Flinn ratio less than 1) and prolate (Flinn ratio greater than 1) ellipsoids are represented, with volumes ranging from the smallest in the ensemble (Particle 48) to those closer to the average (Particle 165). Remanence states were obtained after applying a saturating field of 1 T along the X, Y, and Z directions of the reconstructed volume, and then stepping the field down to zero in steps of 10 mT. The squareness M_rs/M_s is the magnitude of the total remanence vector normalized to the saturation moment of the particle. To aid comparison and to give a better sense of the magnitude of the magnetic moment of each particle, we define “relative Mr” as the total remanent moment of the particle divided by the total remanent moment of a uniformly magnetized 25 nm diameter sphere. This size corresponds to the upper threshold for single-domain (SD) Fe (Muxworthy and Williams 2015).

Table 1

Summary of micromagnetic simulations

Particle ID	Major radius (nm)	Inter. radius (nm)	Minor radius (nm)	Flinn ratio	Mesh vol. (nm3)	Mrs/Ms X Y Z	Hc (mT) X Y Z	Rel. M3	Min. irrev. field (mT)	Domain state
STEM	119	117	98	0.85	5.72E+06	0.0086	0.2	6	270	II
						0.0086	0.5	6	330	II
						0.0086	7	6	260	II
48	36	35	23	0.68	1.43E+05	0.066	42	1.16	260	I
						0.065	19	1.14	310	I
						-0.067	-20	1.16	300	I
321	54	44	28	0.78	2.97E+05	0.032	27	1.15	380	I
						0.036	2.1	1.32	260	I
						-0.032	-4.4	1.15	410	I
364	34	26	16	0.8	2.46E+05	0.042	37	1.27	290	I
						-0.03	-5	0.9	300	I
						-0.042	-4	1.27	>1000	I
8	138	79	54	1.2	2.24E+06	0.012	11	3.1	440	III
						0.022	8	6.1	260	III
						0.022	5	6.1	240	III
233	536	163	71	1.43	2.27E+07	0.013	8	36.7	230	III
						0.016	4	44.8	60	IVa
						0.0035	0.5	9.7	100	IVb
75	259	133	84	1.23	1.13E+07	0.0059	6	8.1	430	III
						0.016	4	21.7	260	III
						0.018	6	25.8	40	IVa
165	262	228	151	0.76	3.73E+07	0.0027	2.5	12.2	250	II
						0.0035	0.64	16	330	II
						0.0035	0.04	16	490	II
155	213	198	73	0.4	1.33E+07	0.004	3.6	7	500	I
						-0.003	-0.09	5.6	260	I
						0.004	1.1	7	300	I

No particles were small enough to adopt an SD state. Instead all particles adopt either pseudo-single domain (PSD) or emerging multi-domain (MD) states, consisting of either a single vortex or multiple vortex/wall-like structures (Fig. 9). To highlight the orientation and nature of vortex cores, the magnitude of the vorticity of the moment vector field (| ∇ × M|) is plotted as an isosurface (green in Fig. 9). The choice of isosurface magnitude is somewhat arbitrary (too large and only the ends of the core are highlighted; too small and the surface extends too far from the core region). We chose the largest value that would produce a continuous trace of the core from one surface termination to another. The remanent states of the particles studied can be divided into four general categories: (1) single vortex with core oriented close to the minor axis of the best-fitting ellipsoid (Fig. 9a); (2) single vortex with core oriented within the plane defined by the minor and major axes of the best-fitting ellipsoid (Fig. 9b); (3) single vortex with core oriented parallel to the major axis of the best-fitting ellipsoid (Fig. 9c); and (4) multiple vortex/wall-like structures containing two or more cores (Fig. 9d).

Figure 9

The four types of pseudo-single domain state observed in the FEM simulations. Color bar for arrows shows the direction and magnitude of the normalized M_z component of the particle magnetization. Green isosurfaces are selected to highlight the magnitude of the vorticity of the moment vector field. These vary particle to particle and are selected to present the minimum continuous surface. (a) Type I: single vortex with core oriented parallel to the minor axis of the best-fitting ellipsoid for particle 155. (b) Type II: single vortex with core oriented within the plane defined by the minor and major axes of the best-fitting ellipsoid for particle 165. (c) Type III: single vortex with core oriented parallel to the major axis of the best-fitting ellipsoid for particle 75. (d) Type IV: multiple core/wall-like structures for particle 233 with the magnetic field applied along the Y and Z-axis, respectively.

Type I behavior was typically observed in oblate particles with lower Flinn ratios (flattened ellipsoids). Type I particles display straight cores located at the center of the largest face. The relatively strong demagnetizing field of the vortex core in this case can be shielded by antiparallel spin tilting at the outer rim of the oblate particle. Type II behavior was typically observed in oblate particles with higher Flinn ratios (triaxial ellipsoids). Type II particles display curved cores that adopt a sigmoidal trajectory through the center of the particle. The ends of the core lie normal to their surface terminations, as requested by the micromagnetic boundary conditions (e.g., Hubert and Schäfer 1998). Type III behavior was observed in the three prolate particles. Type III particles display cores that track the major ellipsoid axis in the central section of the particle. The ends of cores again lie normal to their surface terminations, causing deviations and distortions of their trajectory. Type IV behavior was observed in the two large prolate particles and already represents a diamond-shape domain pattern with preference for 90° walls, as previously observed for iron thin films and whiskers (e.g., Fig. 5.94 in Hubert and Schäfer 1998). The remanent state of these particles was more sensitive to the direction of the saturating field than oblate particles with similar volume. In the case of Particle 233 (the largest and most elongated prolate particle), states III, IVa, and IVb were adopted for saturating fields applied along X, Y, and Z, respectively (Table 1). State IVa contains two curved cores corresponding to Bloch walls that split the particle into magnetic domains; IVb is an efficient flux closure structure containing 10 magnetic domains (Fig. 9d).

The morphology of the vortex core evolves with increasing volume of particle (Fig. 10). Small particles contain well-defined cylindrical cores (Fig. 10a). Larger particles develop cores with a “winged” structure, with wings protruding along the directions of emerging domain walls (Fig. 10b). In larger prolate particles, the core is poorly defined, becoming flattened and developing off-shoots (Fig. 10c) and loops (Fig. 10d), suggestive of emerging MD behavior.

Figure 10

Evolution of vortex core morphology. (a) The vortex core of Particle 48 plotted for a helicity level of 180. (b) The vortex core of Particle 165 plotted for a helicity level of 40. (c) The vortex core of Particle 8 plotted for a helicity level of 80. (d) The vortex core of Particle 75 plotted for a helicity level of 75.

The emerging walls can be better defined and visualized as isosurfaces of the magnetocrystalline anisotropy energy, which highlights those regions where the magnetization points away from the <100> magnetocrystalline easy axes (Fig. 11). As before, the choice of isosurface magnitude is somewhat arbitrary (too large and the domain walls develop holes; too small and the walls become unreasonably wide). We chose a value that generated the thinnest continuous wall structures. Examining the remanent state for each of the orthogonal magnetization directions, a variety of domain wall-like behaviors can be identified. Figure 11a depicts the Type III remanent state observed after applying saturating fields along X. A cross section through the anisotropy surface of this state is shown in Figure 11c. The cross section reveals the presence of four wall-like structures associated with the rotation of spins around a central vortex core. Due to its alignment with a perpendicular easy axis, the vortex core appears as the cylindrical hole in the center of the anisotropy surface in Figure 11c. Near the particle center, the domain walls tend to be thinner and better defined. Walls broaden toward the particle surface, as the magnetization adapts to the surface morphology, driven by the need to avoid high magnetostatic energy. Figures 11b and 11d show Type IVa and IVb behavior, revealing the presence of well-defined 90° domain walls. Again, there is a tendency for domain walls to be thinner in the grain interior, broadening and adopting a more Néel-like character as the spins adapt to the surface morphology (Hubert and Rave 1999).

Figure 11

Thresholded anisotropy energy surfaces for Particle 233, showing wall regions corresponding to the highest magnetocrystalline anisotropy energy. The color of each surface represents the projection of wall magnetization along the minor axis of the best-fitting ellipsoid (~ parallel to X). Red is positive, blue is negative. (a) Type III remanence state resulting from magnetic field applied along the X-axis. Dashed green line shows approximate position of cross section in c. (b) Type IVa remanence state resulting from magnetic field applied along the Y-axis. (c) Cross section through the anisotropy surface in a. (d) Type IVb remanence state displaying well-defined 90° walls resulting from magnetic field applied along the Z-axis.

Remanence vectors with respect to particle shape

Stereograms showing the predicted remanence vectors obtained after applying fields along X (black circle), Y (black square), and Z (black triangle) are shown in Figure 12 for selected particles. Also shown are the orientations of the minor, intermediate, and major axes of the best fitting ellipsoids (blue, green, and red diamonds, respectively), and the corresponding traces of the planes normal to these directions. For single-vortex states, the remanence vector is dictated primarily by the magnetization of the vortex core, but the precise remanence direction is modified significantly by switching the sense of rotation of the bulk spins. Denoting the core direction as either up or down and the sense of bulk spin rotation as being either left or right, this results in two pairs of distinct remanence directions for each particle: (up left and down right) and (up right and down left) (Fig. 12a). The remanence states in each pair are antiparallel to each other, but lie at some angle to the other pair. Remanence directions for two Type I particles are shown in Figures 12b and 12c. The angle between remanence pairs is 32° for Particle 364 (Fig. 12b) and is 80° for Particle 155 (Fig. 12c). In Particle 364 (a small oblate particle), the remanence lies close to the core direction, i.e., close to the minor axis of the best-fitting ellipsoid (Fig. 12b). Note also that the remanence states obtained for this particle after applying a saturating fields along X and Z are identical (only the X state is plotted in Fig. 12b). In Particle 155 (the most extreme oblate particle), the remanence obtained after magnetizing along Y lies much closer to the major axis of the best-fitting ellipsoid than the minor axis (Fig. 12c). A Type II particle is shown in Figure 12d. The angle between remanence pairs is small (16°) and the remanence lies at an intermediate angle between the minor and major axes of the best-fitting ellipsoid, parallel to the average core orientation. A Type III particle is shown in Figure 12e. The angle between remanence pairs is 56°, with the remanence lying close to the minor axis for fields applied along X and close to the major axis for fields applied along Z. In Type IV particles (not shown), the remanence generally lies closest to the major axis of the best-fitting ellipsoid. Those particles that show a change of domain type with field direction (Table 1) display a correspondingly wider range of possible remanence directions.

Figure 12. (a)

Schematic drawing depicting the four possible combinations of bulk spin rotation and vortex core direction. (b-e) Stereograms showing the remanence vectors obtained after applying fields along X (black circle), Y (black square), and Z (black triangle). Diamonds show the orientation of major (red), intermediate (green), and minor (blue) axes of the best-fitting ellipsoid. Lines show the traces of planes normal to the corresponding axes of the best-fitting ellipsoid. Closed symbols and solid lines are lower hemisphere; open symbols and dashed lines are upper hemisphere. Results are shown for two Type I particles (b and c); (d) a Type II particle; and (e) a Type III particle. Note also that the remanence states obtained for this particle after applying a saturating fields along X and Z are identical (only the X state is plotted in b).

It is possible to access each of the four states by applying a suitably oriented saturating field. For example, applying saturating fields to Particle 155 along X and Y switches the sense of vortex rotation while retaining the direction of core magnetization (Fig. 12c). In small fields, the four states are separated by energy barriers that could, in principle, be overcome by thermal fluctuations. However, the energy barriers associated with switching the sense of vortex rotation while retaining the core direction (or switching the core direction while retaining the sense of vortex rotation) are likely to be very high compared to the barriers associated with switching both together. This is because the former process would require considerable internal disruption to the micromagnetic state and a correspondingly high exchange energy penalty, while the latter can be achieved simply by 180° rotation of the micromagnetic state against the shape anisotropy of the particle. Calculating these energy barriers is the next computational challenge, and will ultimately enable the acquisition of remanence during cooling of these particles to be modeled.

Hysteresis loops

Hysteresis loops are shown in Figure 13 for selected particles and applied field directions. The magnetic response of all particles is dominated by reversible magnetization processes (e.g., the rotations of bulk spins toward the field). The reversible component of magnetic susceptibility is highest (lowest) for fields applied parallel to the major (minor) axis of the best-fitting ellipsoid. Irreversible magnetization processes (e.g., nucleation of vortices, irreversible switching of vortex core position, core orientation or core magnetization, changing sense of bulk spin rotation, denucleation of vortices) produce small steps in magnetization superimposed on the large reversible component. This leads to loops characterized by very low values of coercivity (H_c) and squareness (M_rs/M_s). Highest coercivities (30-40 mT; Table 1) are observed in the small Type I particles when fields are applied along X (e.g., antiparallel to the core direction, close to the minor axis of the best-fitting ellipsoid; Fig. 13a). Typical coercivities are of the order of a few milliteslas or less (Fig. 13d). Negative values of coercivity listed in Table 1 highlight an unusual behavior (e.g., Fig. 13c), whereby the upper branch of the hysteresis curve reaches the M = 0 axis at a positive applied field. This behavior leads to a self-reversal of saturation isothermal remanent magnetization (SR-SIRM), in which a component of saturation remanent magnetization is antiparallel to the saturating field direction.

Figure 13

Calculated components of magnetization parallel to X (blue), Y (green), and Z (red) as a function of applied field for the upper (solid) and lower (dashed) branches of the major hysteresis loop. (a) Type I particle (364) with field applied along X. (b) Type I particle (48) with field applied along Y. (c) Type I particle (48) with field applied along Z. (d) Type III particle (233) with field applied along X.

Despite the low M_rs/M_s values, the large volume of the particles means that their total moments are at least equivalent to that of a 25 nm diameter SD particle, and in many cases significantly greater (relative M_r values vary from ~1 up to ~44; Table 1). Despite the low H_c values, the remanence states are also highly stable with respect to applied fields. We define the stability of a remanence sate in terms of the “minimum irreversible field” (Table 1). To calculate this field, the spin state obtained at each negative field of the upper hysteresis branch (from −10 to −1000 mT) was chosen as an initial configuration, and the micromagnetic energy was minimized under zero field. The minimum irreversible field is the smallest negative field required to produce a change in the remanent state of the particle. Minimum irreversible fields are typically several hundred milliteslas, with several showing remanence states stable to more than 400 mT and one simulation showing no significant change even up to 1000 mT. Note that the coercivity of remanence (H_cr) is greater than or equal to the minimum irreversible field. Based on the values in Table 1, all particles would plot in the MD region of a Day-Dunlop plot (Dunlop 2002). However, such a comparison grossly misrepresents the remanence-carrying potential of these particles. Although the ratio of H_cr/H_c is similar to MD samples, which make notoriously poor paleomagnetic recorders, the absolute values of Hcr are at least an order of magnitude larger. Furthermore, MD materials have vanishingly small values of the minimum irreversible field, since irreversible changes to the remanent state of an MD material can be achieved by movement of weakly pinned domain walls. For the vortex states studied here, no changes whatsoever are observed in the remanence state until fields of several hundred milliteslas are applied. A better comparison of the remanence carrying potential of vortex states is obtained by first-order reversal curve (FORC) analysis, which is not adversely affected by the large component of reversible magnetization associated with the rotation of bulk spins. The range of irreversible fields calculated here corresponds very well to the distribution of irreversible magnetization observed in synthetic dusty olivine samples using FORC diagrams (Lappe et al. 2011, 2013) and to the high stability of natural remanent magnetization with respect to alternating-field demagnetization observed by Fu et al. (2014). In comparison, FORC diagrams of MD materials typically show irreversible magnetization restricted to fields less than a few milliteslas (Church et al. 2011; Lindquist et al. 2015). Unsurprisingly, the lowest irreversible fields observed here are for Type IV particles, which contain more MD-like structures. Even here, though, irreversible fields of 40-100 mT are typical.

Rock magnetism of realistic ensembles

Conventional characterization of the remanence carriers in rocks typically relies on either optical or SEM imaging of polished surfaces, or TEM imaging of thin foils. Although such 2D methods are an essential part of the qualitative characterization process, our numerical simulations emphasize just how important 3D knowledge of particle geometry is for quantitative modeling. Shape and crystallographic orientation of an individual particle controls the orientation of its vortex core, or equivalently the position and type of its domain walls. Near the surface the bulk spins adapt to the faceting of the particle. In symmetric particles, we might expect that the bulk spins cancel each other out and that the total remanence would be dominated by uncompensated spins within the core. A striking example of where this assumption breaks down is shown in Figure 14 (Particle 155, a Type I particle with a pronounced oblate geometry). In this case, the core is parallel to the minor axis of the best-fitting ellipsoid, but the net remanence lies at a large angle to this, and rotates by 80° as the sense of bulk spin rotation changes (Fig. 12c). The explanation for this behavior is the combination of the short length of the core, which reduces its contribution to the net moment, and the uneven length of opposing surface facets, which creates a significant contribution from uncompensated bulk spins. In Figure 14, this unbalancing is highlighted by plotting the anisotropy energy of the domain walls. The width of the walls and the sizes of the four resulting quadrants can be seen to be of unequal sizes. In the configuration shown, there are more spins pointing along +Z due to the larger facet on right than there are compensating spins along −Z due to the small facet on the left. A similar situation is known for vortex states in sufficiently large uniaxial particles, where several metastable magnetization states may exist that are related to edge moments that can be aligned parallel or anti-parallel to the global demagnetizing field (Fig. 7 in Rave et al. 1998). The importance of bulk spins in controlling the remanence of vortex states in particles with realistic morphologies is not generally appreciated.

Figure 14

Thresholded anisotropy energy surface for Particle 155, showing wall regions corresponding to the highest magnetocrystalline anisotropy energy. The color of each surface represents the projection of wall magnetization along the minor axis of the best-fitting ellipsoid (~parallel to X). Red is positive, blue is negative. Note the larger size of the right-hand domain (red) compared to the left-hand domain (blue) is due to the different size of the corresponding surface facets, and leads to a net contribution to the remanence of the particle from the bulk spins.

The sheet-like arrangement of particles within the reconstructed volume (Fig. 4) is expected to generate significant remanence anisotropy. Such anisotropy is well documented in single-crystal paleomagnetism studies (Feinberg et al. 2004, 2005; Fu et al. 2014) and must be corrected before the measured remanence directions can be interpreted quantitatively. Fu et al. (2014) measured the anisotropy of ARM susceptibility for the same grain studied here (Fig. 7), finding normalized values of 0.45, 0.56, and 1 for the minimum, intermediate, and maximum susceptibilities, respectively. The maximum ARM susceptibility is observed for fields applied along Z, which corresponds to the average orientation of major ellipsoid axes. This observation is most consistent with remanence carried by prolate Type III and IV particles and oblate Type I particles such as those shown in Figure 14, and suggests that such particles are more prevalent in regions of the grain outside the reconstructed volume. The fact that the minimum susceptibility axis is not aligned with the average orientation of minor axes is explained by the presence of Type I and II particles, which contribute to the remanence when fields are applied normal to the sheets. A full model of the entire ensemble, taking into account the magnetostatic interaction fields between particles and the distribution of shape anisotropy, would be necessary to fully describe the anisotropic response of the system (Hargraves et al. 1991) as well as accounting for the presence of much larger particles that may occur outside the analyzed volume. Such calculations could serve to improve dramatically the interpretation of single-crystal paleomagnetic studies, and minimize the number of repeat measurements needed to reach a statistical significance equivalent to bulk paleomagnetic studies.

Implications for chondrite paleomagnetism

This study demonstrates that particles in the lower half of the size distribution adopt a single vortex state (Type I and II), with the larger particles adopting Type III and IV (i.e., emerging MD) behavior. Lappe et al. (2011) similarly identified a dominance of single vortex states, but noted also a significant number of SD particles. It appears that the size distribution of particles in synthetic dusty olivine is shifted to smaller values than the natural sample studied here. Nevertheless, given the abundance of single vortex states in both cases, calibrations of non-heating paleointensity methods using synthetic dusty olivine (Lappe et al. 2013) remain largely valid. Similarly, given the high field stabilities observed here for single-vortex states (Table 1), the conclusion that dusty olivine is capable of recording and maintaining a faithful record of pre-accretionary remanence (Lappe et al. 2011) also remains valid. Indeed, the combination of high switching fields and large volumes of single-vortex states should translate to carriers with very high thermal stability (i.e., high blocking temperatures). High thermal stability of vortex states at temperatures up to the Curie temperature has recently been demonstrated in magnetite using electron holography (Almeida et al. 2014, 2016). In principle, micromagnetic simulations can be used to calculate energy barriers between alternate remanence states, and thereby model the acquisition of thermoremanent magnetization. Such calculations are beyond the scope of the present study, but are an obvious next step in the nanopaleomagnetic approach.