Mapping Magnetic Signals of Individual Magnetite Grains to Their Internal Magnetic Configurations Using Micromagnetic Models.

Micromagnetic tomography (MMT) is a technique that combines X-ray micro computed tomography and scanning magnetometry data to obtain information about the magnetic potential of individual grains embedded in a sample. Recovering magnetic signals of individual grains in natural and synthetic samples provides a new pathway to study the remanent magnetization that carries information about the ancient geomagnetic field and is the basis of all paleomagnetic studies. MMT infers the magnetic potential of individual grains by numerical inversion of surface magnetic measurements using spherical harmonic expansions. The magnetic potential of individual particles in principle is uniquely determined by MMT, not only by the dipole approximation, but also more complex, higher order, multipole moments. Here, we show that such complex magnetic information together with both particle shape and mineral properties severely constrains the internal magnetization structure of an individual grain. To this end, we apply a three dimensional micromagnetic model to predict the multipole signal from magnetization states of different local energy minima. We show that for certain grains it is even possible to uniquely infer the magnetic configuration from the inverted magnetic multipole moments. This result is crucial to discriminate single-domain particles from grains in more complex configurations such as multi-domain or vortex states. As a consequence, our investigation proves that by MMT it is feasible to select statistical ensembles of magnetic grains based on their magnetization states, which opens new possibilities to identify and characterize stable paleomagnetic recorders in natural samples.

Introduction

Paleomagnetic studies commonly try to extract from geological samples the primary characteristic remanent magnetization, acquired during rock formation (Butler, 1991). These studies put considerable effort in selecting rocks with a narrow magnetic grain size distribution. These are considered to have optimal recording properties, and are preferred because even careful stepwise demagnetization of bulk magnetic measurements cannot always reliably separate primary from secondary, or stable from unstable remanence carriers (Dunlop & Özdemir, 1997). The recording properties of larger grain sizes, for example, of those that dominate the magnetic signal of lavas, are often problematic and poorly understood. One way to minimize the influence of unwanted remanence carriers, pioneered by Tarduno et al. (2001, 2006), is to mechanically pick from the whole rock only those fractions that contain reliable magnetic carriers. Another option, which is only now becoming technically feasible, is to recover the magnetization of all magnetic grains individually together with exact information about position and shape of the carriers. The statistical ensemble for paleomagnetic information can then be selected based on this information.

The analysis of the magnetic signal of individual magnetic grains in rock samples has progressed significantly with modern micro‐to nano‐scale scanning magnetometry techniques such as superconducting quantum interference device (e.g., Egli & Heller, 2000; Lima & Weiss, 2016; Weiss et al., 2007) microscopy or quantum diamond microscopy (QDM; e.g., Farchi et al., 2017; Fu et al., 2020; Glenn et al., 2017; Levine et al., 2019). In particular, by combining the surface measurements of the magnetic field of a rock sample with the spatial location and geometric profile of individual grains, which, for example, can be obtained via X‐ray computed tomography (microCT), it is possible to recover dipole and higher order magnetic moments from individual grains. This method has been demonstrated in practice by de Groot et al. (2018) with a technique known as micromagnetic tomography (MMT; de Groot et al., 2021), in which the magnetic dipole moments of individual magnetic particles in a synthetic sample were successfully recovered by means of numerical inversion of both scanning and microCT data. Furthermore, it has been proved by Cortés‐Ortuño et al. (2021) that higher order moments can also be retrieved, especially for shallow grains, which usually have a strong signal. These higher order solutions are unique for topologically separated grains (Fabian & de Groot, 2019) and are recovered by approximating a grain by a point magnetic source (Cortés‐Ortuño et al., 2021). By taking into account higher order moments in the numerical inversions, the error between the numerically inverted field and the observed (measured) field decreases substantially when grains that produce complex magnetic signals with multipolar character are present in the sample. In consequence, the magnetic moments, and thus the magnetization, of these higher order signal grains are recovered with more accurate and reliable solutions. If it is the case that grains are mostly homogeneously magnetized, the alternative method of de Groot et al. (2018), where grain volumes are decomposed into cuboids to obtain the dipole moments of the grains, can produce correct solutions with fewer free parameters to be fixed by the observed data. Although this is difficult to occur in natural samples due to the large size distribution of the grains, reasonable dipolar moment solutions can still be retrieved by using a selection criteria to accept adequate grain solutions, as shown by de Groot et al. (2021).

A fundamental limitation of the MMT numerical inversion is that it cannot resolve the exact magnetization structure within a grain, because the magnetic potential outside any spherical volume can be created by equivalent charges on the sphere surface. Therefore, further information about the internal magnetic configuration of individual grains is not present in the surface signal. On the other hand, the inverted dipole and higher order magnetic moments per grain are in principle unique and this information effectively constrains the physically possible magnetic state in the grain. This constraint can be sufficient if all local energy minima (LEM) of the grain system are known and only one of them corresponds to the observed multipole moments. The main motivation for unveiling the internal magnetic configuration of individual magnetic grains is the possibility to select statistical ensembles of stable paleomagnetic recorders by classifying particles with similar properties. This is applicable to natural samples because they contain a wide range of grain sizes and geometries where a large number of different magnetic domain states can be hosted as stable configurations (Nikolaisen et al., 2020).

If the shape and mineral properties of a magnetic particle is known, its physically possible stable magnetization structures can be modeled by numerically solving the micromagnetic equations that describe both the magnetic interactions experienced by the magnetic moments in the particle and the equilibrium conditions for a stable state (Ó Conbhuí et al., 2018). The magnetization structure of the grain in this approach is represented by a continuous unit vector field, from which all relevant magnetic energy contributions within the particle, such as exchange, magnetocrystalline anisotropy, or demagnetizing energy can be calculated. Physically stable configurations correspond to local minima of the total energy. In this work, grains are modeled numerically by means of the Micromagnetic Earth Related Robust Interpreted Language Laboratory (MERRILL) finite element (FE) code (Ó Conbhuí et al., 2018) in order to: (a) obtain different LEM states in magnetite grains of different size and geometries. (b) Model the stray field signal of the grain which is directly linked to its internal magnetic configuration. (c) Infer the magnetic state of the grain through numerical inversion of the modeled stray field signal, which allows to recover the multipole moments of the grain system (see Figure 1). In the case of (a), the simplest geometry considered here is a 40 nm radius sphere where the LEMs are vortex and single‐domain states that can be completely classified. More complex geometries are then studied by modeling particles from synthetic (de Groot et al., 2018) and natural (Nikolaisen et al., 2020) samples. For the size ranges considered in this study, three‐dimensional vortex configurations are generally observed as LEM. Such states are known to be good paleomagnetic recorders because of their high magnetic stability (Nagy et al., 2017). Regarding steps (b) and (c), these are directly comparable to experimental evidence since microscopy data measures surface magnetic flux data from the stray field signal of the grains. By only using these magnetic surface data together with magnetic mineralogy, location, and shape of a grain, this study demonstrates that in some cases it is possible to uniquely deduce the internal magnetic state. Thus, the information contained in the inverted dipole and higher order moments, together with the physical constraints for micromagnetically stable states, can completely determine the magnetization structure.

Figure 1

Summary of the technique of micromagnetic modeling of grain data from micromagnetic tomography (MMT) developed in this study. Tomography data is applied to reconstruct grain models and simulate the grain signal based on the internal magnetic states. Scan data from MMT is numerically inverted into the grain locations from micro‐CT. Magnetic moments from both approaches are compared to constrain the possible magnetic states in the MMT data.

A summary of the methodology developed in this investigation is shown in Figure 1. The comparison of magnetic moment values of experimental data with the simulations is restricted by the resolution of the tomographic data that allows to make realistic models of the sample grains. Grain sizes in rock samples are above the micrometer scale which are still challenging to numerically simulate. Development in tomographic techniques, such as nano‐CT, high‐resolution magnetometry, such as scanning QDM, and code optimization, might allow to validate the proposed micromagnetic modeling framework for large scale natural grains.

Methods

Micromagnetic Modeling

To numerically model magnetic grains we apply micromagnetic theory to describe both the material and the magnetic interactions in the continuum. Magnetic moments in the material particle volume are treated as a continuous unit vector field: where S 2 denotes the surface of the three‐dimensional unit ball. In real space, a closed unit ball centered at the origin is defined as , with being the position vector and denoting the Euclidean norm. Local energy minimum (LEM) states are obtained by minimizing the discretized total magnetic energy functional using the FE micromagnetic code MERRILL (Ó Conbhuí et al., 2018). Within the FE technique, the grain model is discretized into a mesh using tetrahedral elements and a magnetization unit vector m with ‖m ‖ = 1 is specified at every node of the FE mesh, as shown in Figure 2.

Figure 2

Overview of the micromagnetic modeling of an individual magnetic grain. The grain geometry is represented by a finite element mesh that consists of a volume V partitioned into tetrahedrons. At every node i of the tetrahedral mesh a magnetization unit vector m is defined, which is determined by two spherical orientation angles (θ, φ). At every measurement point at the scan surface (sensor), the magnetic signal is approximated by the sum of the normal‐to‐the‐surface magnetic field component B from all node dipole moments , with M s as the saturation magnetization and V as the node volume. This simulated magnetic signal is numerically inverted for the multipole expansion coefficients of a point source located at the centroid of the grain geometry. The coordinate system is shown at the bottom left.

After a LEM is found, the stray field signal at a scan surface sufficiently far away from the magnetic grain is approximated by the sum of the signals of all magnetic dipole moments (the following notation is based on the rank‐n multipole tensor Θ ( which has 2n + 1 independent components (Cortés‐Ortuño et al., 2021)

belonging to the nodes of the FE mesh, where M s is the material saturation magnetization, and V the volume represented by node i. This volume is specified as the sum of a quarter of the tetrahedron volume of every tetrahedron adjacent to the FE node where the dipole is located (Ó Conbhuí et al., 2018). According to this definition, the total magnetic dipole moment from a micromagnetic simulation can be computed by summing all node dipole moments or, equivalently, using the weighted mean of dipole moments as with the total volume of the FE mesh The magnetization at every node is defined as

and the mean magnetization value of the sample is specified as:

The norm of M sim, that is, the grain magnetization from the simulation, will be used as a benchmark against the results of inversions.

The magnetic field B(s) at some scan locations s is then approximated as the sum over all L nodes of the mesh:

This surface signal is discretized into N × N sensor points (see top of Figure 2) which each contains the field contribution from all the node magnetization vectors. For this study, the field component B normal to the scan surface is chosen to make the results comparable to experimental results (de Groot et al., 2018). According to this, the surface signal is defined in the matrix . In the MMT framework, the surface data matrix , the grain location, and the grain shape are the available information to analyze the magnetic properties of the grain via a numerical inversion of the magnetic signal. Multipole inversion is performed by representing the grain as a magnetic point source, which is the center of a multipole expansion of the magnetic potential of the particle, as described by Cortés‐Ortuño et al. (2021). The shape information in this case only ensures that the modeled particle is topologically separated from the other sources, such that theoretical uniqueness of its potential is ensured by the unique‐source assignment theorem (Fabian & de Groot, 2019). The inversion returns estimates of dipolar and higher order magnetic moments of the grain. By considering higher order moments, like quadrupolar or octupolar moments, the potentials of complex inhomogeneous magnetization structures are represented more accurately. This is observed as a decrease of the residual between the directly modeled and the inverted magnetic field signal (Cortés‐Ortuño et al., 2021). In the multipole inversion, the magnetic grain is considered as a point source because the non‐uniqueness of the potential field inversion inhibits the method to directly solve for the internal magnetic structure of the particle. Hence, the obtained information about the magnetic states of the particle can only be classified through magnitude and orientation of their multipole moments. The dipole moment together with the total volume V of the particle is used to calculate the average magnetization of the magnetic grain. Higher order moments cannot be used directly to further interpret the signal of complex magnetic states because they also depend on the choice of the expansion center. However, these terms are fundamental to improve the fit of the measurement signal, and thus contribute significant information about the deviation from a pure dipole signal.

The following sections discuss the application of the micromagnetic modeling and multipole inversion to different grain systems. In order to verify that the calculations made by the code for the grain systems are correct a test model based on a cuboid grain has been developed in Section S1 in the Supporting Information S1, where analytic and numerical solutions are compared.

Multipole Inversion of Magnetite Spheres

For a proof of concept, we use the micromagnetic states of a magnetite sphere with 40 nm radius at −100°C. At this low temperature, the magnetic material parameters, which are summarized in Table 1 and obtained from the MERRILL code (Ó Conbhuí et al., 2018), lead to only a small number of LEM states allowed by the system. Moreover, at this temperature, magnetite still possesses cubic magnetocrystalline anisotropy with easy axes oriented in the [111] directions. For simplicity, stress effects in the material are not considered. The FE model of the mesh is constructed such that the edge lengths of the tetrahedra are smaller than the exchange length (see Table 1). The exchange length parameter determines the scale at which atomic exchange interactions dominate over the stray field interactions, which means the scale where the magnetization field changes its direction. Hence, it is a characteristic length of magnetic inhomogeneities such as domain walls. Using edge lengths larger than ℓ ex (see Section 2.4) is suited for large samples where magnetic domains such as vortices, arise from the effects of the demagnetizing field, and where inhomogeneities are only approximated. Furthermore, this length parameter influences switching fields in dynamic processes, which are of no interest here. Other exchange lengths are defined for systems with uniaxial anisotropy or interacting with a field. For more details see Andreas et al. (2014) and Kronmüller and Fähnle (2003).

Table 1

Material Constants of Magnetite at −100°C Obtained From MERRILL (Ó Conbhuí et al., 2018)

Magnetic constant
Exchange	A	14.343	pJ m−1
Cubic anisotropy	K 1	−9.490	kJ m−3
	K 2	−1.921	kJ m−3
Saturation magnetization	M s(−100°C)	0.502	MA m−1
Exchange length	ℓ exch	9.517	nm

For the magnetite sphere mesh model the 99% of tetrahedron edge lengths are smaller than 6.52 nm with a mean value of 4.61 nm (see Figure S3 in Supporting Information S1 for details). The mesh was modeled using the Gmsh (Geuzaine & Remacle, 2009) code and translated to a MERRILL input format using the meshio (Schlömer, 2021) library. To systematically find a large number of LEM states, the system is initialized with three sets of initial conditions:

Two hundred magnetic configurations with magnetization vectors randomly oriented.

Using 26 uniformly magnetized states oriented in the directions of the Cartesian axes and their diagonals for each plane (which are eight directions per each of the xy, yz, and xz planes), and the eight 〈111〉 directions. A small random component is added to these states.

The same 26 uniform states but with a larger random perturbation to the magnetization vectors (see Figure S4 in Supporting Information S1).

Cases (2) and (3) are run eight times each. From the initial state, the energy of the system is then minimized with a conjugate gradient method in Cartesian coordinates which finds the LEM that is closest to the given initial state. As a result, we obtain a set of LEM configurations that can be classified in two categories: uniform or single‐domain states, and vortex states. Vortices are inhomogeneous configurations characterized by a flux closure pattern of the magnetization vectors around a core which tries to minimize the stray field energy. Snapshots of these configurations are shown in Figure 3, where vortex states are visualized by plotting an isosurface of the helicity of the unit vector magnetization field that is defined by hel(m) = m ⋅ ∇ × m. For both, the uniform and the vortex states, only a finite number of orientations are LEM states as the spherical symmetry is broken by the minimization of the cubic anisotropy energy. Because K 1 < 0, no uniform state in any of the 〈100〉 directions is found. In contrast, the total moment of a vortex state can be oriented along the 〈100〉 axes (e.g., state 22 in Figure 3). It occurs that, even though the magnetization vectors in the vortex core are aligned close to a hard anisotropy axis, the overall better alignment of the surrounding magnetization vectors with the easy axes reduce the total field energy. This magnetization alignment can be visualized using the local anisotropy energy density, or anisotropic deviation of moment, as defined by Nagy et al. (2019). In larger particles, the vortex core is not fully aligned with the hard anisotropy axes but forms three tube‐like maxima of the anisotropic deviation of moment along the core (Nagy et al., 2019; Witt et al., 2005).

Figure 3

Dipole classification of the local energy minima (LEM) states of a spherical grain as inferred from multipole inversions. The top diagram summarizes the inverted dipole orientations of different LEM states of a magnetite sphere of 40 nm radius. The inversion includes multipoles up to octupole order, and is based on a modeled scan signal at the surface height z = 0.6 μm, whilst the sphere center is located at z = −0.2 μm. The polar (black) and azimuthal (red) angles represent the mean inverted dipole orientations from multiple simulation runs that were grouped into 42 different LEM states. Snapshots of these 42 LEMs are plotted below the diagram. Vortices are depicted by isosurfaces of their helicity scalar field. Uniform and vortex states are colored by the z‐component of the unit magnetization vector. The orientation of a uniform state is indicated in cases where arrows are not visible.

After the energy minimization of the magnetic states, their stray field signal is computed at a symmetric scan surface above the spherical sample at the origin. The scan surface is defined by a 1.5 × 1.5 μm grid of sensor points in the xy‐plane with a spacing of 10 nm in both the x and y‐directions, while the scan height is varied between 200 nm and 1.2 μm in steps of 200 nm. Snapshots of the scan signals are shown in Figures 4a and 4b for a vortex state and a uniform state, respectively. The signals of these two states look similar when the magnetization field is, in average, oriented in the same direction, with the difference that the field of the uniform state is stronger. Multipole inversions of the scan signals into a magnetic point source defined at the geometric center of the sample are now obtained by applying the methodology of Cortés‐Ortuño et al. (2021) for each of the 564 simulation results at the different scan heights and using three different multipole inversion orders.

Figure 4

Magnetic signal and inverted magnetization for a vortex and a uniform state in a magnetite sphere. Panels (a) and (b) depict the magnetic signal of the magnetic configurations simulated with micromagnetic modeling at a height of 0.6 μm. These signals are used in the octupole order inversions of Figure 3. Panels (c) and (d) summarize the mean magnetization obtained from numerical inversions, up to different multipole order, as a function of scan height. These magnitudes are scaled by the saturation magnetization of magnetite at −100°C. The averaging is done for states with similar inverted magnetization angle and norm, which here are classified into state 11 (a) and 33 (b). The scan height specified in both (a) and (b) is marked with a dashed line. The angle deviation and relative error of the mean inverted magnetization (compared to the expected simulation result) for an octupole inversion at this height is indicated at the bottom of (a) and (b).

The results of inversions of the signal of a vortex and uniform state are shown in Figures 4c and 4d, respectively. The inversions are computed for increasing multipole order and as a function of scan height by plotting the magnetization of the configurations in normalized units. The scan surface heights are specified starting at z = 0 with the top of the magnetite sphere at z = −0.16 μm, although only its geometric center is used for the inversions. Every data point results from averaging the magnetization from simulations with similar inverted angle and magnetization norm, as explained in the next paragraph. In the case of vortex states, the magnetization, which is computed from the inverted dipole moments, converges to a magnitude that is a fraction of the saturation magnetization of magnetite, namely 0.68 M s(−100°C), as the scan height increases. This trend occurs for any multipole inversion order (see Figure 4c). In contrast, for uniform states the magnetization converges toward the saturation magnetization. This result confirms that inversions in real grains producing magnetizations smaller than the saturation magnetization of the material are an indication that the grains are likely to have an internal multi‐domain or complex magnetic configuration. The convergence of the inversions are validated by the small magnitudes of both the angle deviation and the relative error of the inverted magnetization from an octupole inversion with respect to the expected magnetization from the micromagnetic simulation M sim (also averaged from states with similar magnetization angle and norm), which are shown at the bottom of Figure 4. The inversion results for every state of Figure 3 are specified in Table S2 in Supporting Information S1.

In order to classify the magnetic configurations through the numerical inversions, both the orientation and magnitude of the inverted magnetizations are compared with each other using the results of an octupole inversion at a scan height of 0.6 μm. According to this, configurations are grouped when making an angle smaller than 0.15 radians and when the norm difference is smaller than 5% of magnetite's M s(−100°C). As a consequence, an average polar, θ, and azimuthal, φ, angle are obtained from every group of magnetic states, by applying a circular mean for the azimuthal angle. This resulted in 45 different magnetic configurations. A set of 42 states was selected by discarding 3 higher energy or metastable LEM (see Figure S5 in Supporting Information S1) and are illustrated in Figure 3 using a scan height of 0.6 μm. Energies, number of states selected of every group and inversion errors are specified for all the 42 configurations in Table S2 in Supporting Information S1. In Figure 3, the mean polar angles are sorted in increasing order. Additionally, mean azimuthal angles from groups of states with similar polar angles, for example, states 17–26, are again sorted in increasing order. From the figure, it can be seen that there are states which have equivalent spherical angles, such as states 8 and 9. These cases correspond to uniform and vortex domains with equivalent orientations, which produce similar stray field geometries. These states can be distinguished by the magnitude of the magnetization vector which is nearly 1.5 times smaller in the vortex state as compared to the uniform state, as depicted in Figures 4c and 4d. This difference is due to the flux closure in the vortex state which renders parts of the magnetization invisible in the potential field, as a completely flux‐closed loop is an annihilator of the potential. Another set of configurations that share the same inverted spherical angles but cannot be uniquely recognized are vortex configurations with different chirality. For instance, state 15 in Figure 3 has a left‐hand sense of rotation of the magnetization around the vortex core, but it could equally have a right‐hand orientation, as in state 17. This occurs because the component of the field produced by the magnetic moments, that is oriented with the vortex core, is the only one that is not averaged out in the field signal at the scan surface. Therefore, the inverted dipole moments do not contain any information to distinguish opposite vortex chiralities and it is necessary to analyze the next order, that is, quadrupole, field to identify these configurations. This will be studied in detail in Section 3. More exotic configurations were also found, such as state 14, where two vortices of opposite chiralities are connected by a singularity known as Bloch point. Or the vortex state 27 which is possibly a metastable configuration since the energy minimization did not converge fully under the numerical stopping criteria of the simulations (see Figure S5 in Supporting Information S1). The energy of states containing a singularity is substantially larger than that of the vortex and uniform configurations, and thus they are less stable and less likely to be observed experimentally. There are possibly a few other metastable configurations that were missed during energy minimization from the selected set of initial conditions. For instance, in‐plane uniform states alike to state 21 (which is aligned along the [1-10] direction) should also be observed aligned toward other directions allowed by anisotropy.

Synthetic Sample Grain

A more realistic scenario for the application of the methods developed in Section 2.1 is the modeling of magnetite grains from the synthetic sample studied by de Groot et al. (2018). These particles were prepared and described by Hartstra (1982) and their diameters range between 5 and 35 μm. Both the geometric properties and positions of the magnetite grains are recovered from the three‐dimensional density distribution within the sample, acquired by microCT (Sakellariou et al., 2004). In de Groot et al. (2018), the microCT data is used to model the magnetic particles as aggregations of neighboring equally sized voxels. In a next step, this data is coarsened by grouping the voxels in clusters of cuboids of varying size while keeping the geometric shape of the grain.

In order to model the grain geometries as FE meshes, the vertices of the grain cuboids of de Groot et al. (2018) are used to reconstruct a surface that can be partitioned into tetrahedral units. If the majority of the surface of a grain exhibits a convex curvature, an efficient surface reconstruction is achieved using a convex hull algorithm. In the case of more complex geometries, a sophisticated surface reconstruction method for unstructured point clouds is required, such as the Advancing Front method or the Poisson Surface Reconstruction (PSR) method (Alliez et al., 2021). These algorithms are implemented in the open source Computational Geometry Algorithms Library (CGAL; The CGAL Project, 2021). Consequently, an FE mesh is generated by partitioning the volume enclosed by the surface mesh using tetrahedral cells, which is also known as tetrahedralization. This can be optimally achieved using a Tetrahedral Isotropic Remeshing method (Tournois et al., 2021). The final mesh is converted into a format accepted by MERRILL using the meshio (Schlömer, 2021) library. It should be noted that this reconstruction method can lead to substantial uncertainties in volume and shape of the reconstructed grains, which are due to resolution limitations of the available microCT technology. These problems can be overcome by using high‐resolution tomographic techniques (see Discussion in Section 3.3) and therefore are not essential for the current proof of concept.

One of the grains with a more complex geometry in the synthetic sample of de Groot et al. (2018) is grain 4 shown in area 2 of the scan measurements, thus its shape was chosen for a test model. The original grain has an approximate volume of 3,617.25 μm3 (by summing the cuboid volumes) and is contained in a minimal box of dimensions 13.2 × 16.72 × 13.2 μm3 in the x, y, and z directions, respectively. Because we cannot numerically model a system of these dimensions with the current computational tools and resources, we down‐scaled the grain dimensions to 2% of the original values. The new smaller grain is contained in a box of 0.36 × 0.41 × 0.64 μm3 which allows for the calculations to finish in a reasonable time. Using the down‐scaled coordinates of the cuboid vertices the PSR method was applied to reconstruct a surface mesh for the grain model, which required an estimation of oriented surface normal vectors. In a next step, the volume was tetrahedralized into a mesh with an average tetrahedron edge length of 5.88 nm, which is below the exchange length constant of magnetite. The FE mesh (see Figures 5a and 5b) agrees well with the model using cuboids, where the former has a volume of 0.0265 μm3 (equivalent to a sphere of 0.3699 μm diameter) and the latter a volume of 0.0289 μm3. Details of the mesh construction and statistics of the mesh are discussed in Section S3 in Supporting Information S1.

Figure 5

Micromagnetic modeling results for a down‐scaled magnetite grain 4 of area 2 from de Groot et al. (2018) in a synthetic sample. Panels (a–c) show the finite element mesh obtained from surface reconstruction based on microCT voxel data. The topmost point of the mesh is indicated in (a) as z = −0.11 μm. Each column represents a local energy minima (LEM) magnetic configuration at zero‐field. The corresponding modeled scan signals are depicted in (d–f), for a scan surface at z = 0.4 μm. The grain centroid is indicated by a green dot. Results of multipole inversions up to different orders, and as a function of scan height, are shown in panels (g–i). Here, magnetization is scaled by the room temperature saturation magnetization of magnetite. Dipole moment spherical orientation angles are specified for different multipole expansion orders at a scan height of 0.4 μm, indicated by the vertical dashed line, and the expected magnetization is indicated by a horizontal dashed line.

For this model grain, we obtained a set of six LEM states at zero‐field and room temperature by minimizing the energy of six completely randomized initial configurations. Further information of the minimization process can be found in Section S3 in Supporting Information S1. For these simulations, all the known magnetic interactions in magnetite, which contribute to the total magnetic energy, were considered. These are specified in Table 2 and were obtained from MERRILL (Ó Conbhuí et al., 2018). The three cubic anisotropy axes are defined to be aligned with the x, y, and z‐coordinate axes of the full sample, which correspond to the Cartesian axes labeled in Figure 5b. As a result of the energy minimization all six initial states converged into multi‐vortex configurations where the stray field energy is minimized by the flux closure patterns around the vortices. At this grain size, multi‐vortex domains are likely to be the most characteristic LEM states. Snapshots of three states (see Section S4 in Supporting Information S1 for all six computed states and the corresponding inversion results) are illustrated in Figures 5a–5c, where vortices are displayed as iso‐contours of the helicity scalar field. Darker regions define negative, lighter positive helicity. The normalized magnetization fields are illustrated with streamlines to visualize the magnetic flux around the vortex cores. The top of the modeled grain is at z = −0.11 μm, and scan measurements were calculated at different heights from z = 0 up to z = 1 μm in steps of 0.2 μm. The B signals from the three chosen magnetic states (Figures 5a–5c) at a height of z = 0.4 μm are plotted in Figures 5d–5f. Signals from states b and c exhibit a clear multipole character compared to, for example, Figure 4b. Results of the inversion of these magnetic signals with respect to the centroid of the grain are summarized in Figures 5g–5i, where magnetization values are scaled by the saturation magnetization of magnetite at room temperature M s(20°C). The inversions were computed as a function of the scan grid height and for three multipole expansion orders. In the three cases, the inverted magnetizations are smaller than 9% of M s(20°C). As the scan height increases, the magnetization values converge at different rates toward the values of the octupole expansion which includes all three multipole orders. The effect of higher order magnetic moments on the inversions, for the double vortex state of Figure 5b, is illustrated in Figure 6. Here, the stray field signal matrix, the inverted scan signal matrix and the residual (Zhdanov, 2015) matrix

Table 2

Material Constants of Magnetite at Room Temperature Obtained From MERRILL (Ó Conbhuí et al., 2018)

Magnetic constant
Exchange	A	13.444	pJ m−1
Cubic anisotropy	K 1	−13.039	kJ m−3
	K 2	−3.154	kJ m−3
Saturation magnetization	M s(20°C)	0.482	MA m−1
Isotropic magnetostriction	λ s	0.475	× 10−6
Exchange length	ℓ exch	9.587	nm

Figure 6

Comparison of scan signal inversions between dipole and octupole order expansions at different scan heights. Panels (a) and (b) depict the inverted magnetic field signal and corresponding residual of the double vortex state of Figure 5b, for a scan height of 0.4 μm and using a dipole order inversion. Panels (c) and (d) illustrate similar results but using an octupole order expansion. The pair of panels (f–g) and (h–i) depict equivalent results but from the scan signal at a height of 1.0 μm. Panels (e) and (j) show the scan signals at the two different scan heights.

are computed for two scan heights and for both the dipole and octupole inversions. The residual error can be quantified by the relative error

with denoting the Frobenius norm. The scan signal at a height of 0.4 μm is shown in Figure 6e and the dipole inversion of this signal is shown in Figures 6a and 6b. The dipole inversion exhibits a strong residual (see Figure 6b) of B err = 30.99%. In contrast, for an octupole inversion (see Figures 6c and 6d) the relative error is B err = 3.66% which is an order of magnitude smaller than the dipole case. At a higher scan position of 1.0 μm (see Figures 6f–6i) the same tendency holds: the residual of a dipole inversion produces a relative error of B err = 17.34% (see Figure 6g) while for an octupole inversion this value is substantially reduced to B err = 0.35% (see Figure 6i). This result supports previous claims (Cortés‐Ortuño et al., 2021) that higher order moments are critical to obtain more accurate inverted magnetic moment values for grains carrying a complex magnetic ordering and exhibiting a sufficiently strong stray field signal. Finally, together with the inverted magnetizations, below Figures 5g–5i are shown the spherical angles obtained from the inverted dipole moments for a scan height of 0.4 μm. The most visible change occurs for the state of Figure 5b with both the azimuthal and polar angles varying roughly by 20° for an octupole inversion compared to the dipole inversion.

A summary of the inverted magnetic parameters from the octupole expansion at a scan height of 1.0 μm for the three states of Figure 5 are detailed in Table 3. In addition, expected values for the magnetization and the spherical angles are computed directly from the micromagnetic simulations in order to verify the correct convergence of the inverted values. These results show that each of these magnetic configurations can be uniquely distinguished by their magnetization and orientation angles. Furthermore, the small relative errors B err confirm that the inversion with higher order moments leads to an accurate estimation of these magnetic parameters.

Table 3

Summary of Inverted Parameters for the Three States of Figure 5

		Magnetization			Residual
		Intensity	Direction		Error
State (Figure 5)	Result	|M|/M s(20°C)	θ	φ	B err
Triple vortex (a)	Inversion	0.0391	80.1°	169.7°	0.13%
	M sim	0.0383	79.8°	169.7°
Double vortex (b)	Inversion	0.0234	92.1°	68.7°	0.35%
	M sim	0.0239	91.7°	71.2°
Double vortex (c)	Inversion	0.0524	160.2°	−36.1°	0.06%
	M sim	0.0528	159.6°	−37.0°

Note. Results for the magnetization, orientation angles, and the relative error of the inverted field, are shown for the octupole expansion at a scan height of z = 1.0 μm. Expected values are indicated in the rows.

Natural Magnetite Grains

To further test the micromagnetic modeling method developed in this study, it has been applied to a natural magnetite grain embedded in a pyroxene volume which has been investigated by Nikolaisen et al. (2020). The pyroxene volume has been milled and imaged via focused‐ion‐beam nanotomography (FIB‐nt) from a norite sample obtained from drill cores of the Mount Harcus Intrusion (Nikolaisen et al., 2020). For this study, we analyzed the magnetite particle labeled OPX042 from the public data set of Nikolaisen (2020), which provides FE meshes of all the magnetite inclusions in the pyroxene volume. This grain has a volume of 0.0049 μm3 which is equivalent to a sphere of 0.2107 μm diameter. The tetrahedra for this FE mesh have edge lengths larger than the exchange length, with a median of 11.98 nm (see Section S5 in Supporting Information S1 for details). Hence, at this scale, the focus is on magnetic states dominated by the demagnetizing field. By replicating the workflow used for the other modeled grains, 20 zero‐field magnetic configurations in the OPX042 grain have been obtained using MERRILL by specifying 20 initial states with randomly oriented magnetization vectors at room temperature and then minimizing their total energy. Three of these states were selected and are shown in Figures 7a–7c, which exhibit clear vortex structures (see Section S5 in Supporting Information S1 which shows all the computed states and the corresponding inversion results). In particular, vortices with opposite chirality arise in the states of Figures 7a and 7b. To compute the scan signals and the numerical inversions the coordinate system is defined such that the x, y‐coordinates of the centroid of the grain lies at the origin and the topmost mesh point lies at z = −0.0937 μm. The scan signal of these states at a height of z = 0.6 μm are illustrated in Figures 7d–7f, where the signal of the triple‐vortex configuration (b) is the weakest. Magnetization values obtained from numerical inversions using increasing scan heights are depicted in Figures 7g–7i. The corresponding dipole moment orientations are shown below the magnetization plots for a height of z = 0.6 μm. It can be observed that for a scan height close to the sample, z ≤ 0.2 μm, the recovered magnetizations differ significantly for the three expansion orders, especially for the single and double vortex states. Nevertheless, as the scan height increases the magnetization curves converges to the magnitudes given by the octupole inversion. As for the case of the synthetic sample, a summary of the inverted parameters for an octupole inversion at a height of 1 μm is shown in Table 4. These inverted values are compared with the expected values from the simulations, showing smaller deviations than compared to the ones observed in the synthetic sample grain.

Figure 7

Micromagnetic modeling of a magnetite grain in a natural sample. Panels (a–c) show the mesh model of the OPX042 grain of (Nikolaisen, 2020) obtained from FIB‐nt. The topmost point of the mesh is indicated in (a). Each model represents a LEM magnetic configuration at zero‐field. The scan signal of these states are depicted in the corresponding panels (d–f), using a scan surface located at z = 0.6 μm. The grain centroid is indicated by a green dot. Results of multipole inversions of different orders and as a function of scan height, applied to the scan signals, are shown in panels (g–i) for the three states (a–c), respectively. Here, magnetization is scaled by the room temperature saturation magnetization of magnetite and dipole moment orientation angles are specified for different multipole expansion orders at a scan height of 0.6 μm, as indicated by the vertical dashed lines. The expected magnetizations are indicated by horizontal dashed lines.

Table 4

Summary of Inverted Parameters for the Three States of Figure 7

		Magnetization			Residual
		Intensity	Direction		Error
State (Figure 7)	Result	|M|/ M s(20°C)	θ	φ	B err
Double vortex (a)	Inversion	0.8256	106.1°	70.9°	0.02%
	M sim	0.8258	106.1°	70.9°
Triple vortex (b)	Inversion	0.0770	117.4°	−37.7°	0.05%
	M sim	0.0770	117.5°	−37.6°
Single vortex (c)	Inversion	0.7326	74.2°	−108.0°	0.02%
	M sim	0.7328	74.2°	−108.0°

Note. Results for the magnetization, orientation angles, and the relative error of the inverted field, are shown for the octupole expansion at a scan height of z = 1.0 μm. Expected values are indicated in the rows.

Discussion

This investigation has shown the feasibility of applying micromagnetic modeling to analyze both the magnetic signal of individual grains carrying different magnetic domain states, and recovering their magnetic parameters via multipole inversion. The parameters obtained from grain simulations, such as the dipolar moments, can be directly compared with those obtained from the numerical inversion of the empirical microscopy data, in order to reveal the internal magnetic configurations of the grains.

Single Domain States and Quadrupole Moments

In geometrically simple particles, like the magnetite spheres modeled here, it is possible to uniquely detect the magnetic configuration in the sphere from both the inverted magnetization and the dipole moment orientation. For this particular case, a uniform state and a vortex configuration are the two LEM domain states observed, which could be oriented in different crystallographic directions while keeping the same energy. Accordingly, although there can be infinitely many LEM in the form of vortex or uniform domain states, the cubic anisotropy in magnetite constrains the domain orientations into a finite number of directions. Because of the energy equivalence of either the vortex or uniform domain regardless of their orientation, the inverted magnetization has always the same value and, as expected, the vortex state has associated a smaller magnitude. In consequence the dipole moment orientation is crucial to distinguish the magnetic domain orientation. Nonetheless from the dipole tensor, it is not possible to discriminate between vortices with similar orientation but opposite chirality. As mentioned in Section 2.2, it is not the dipole field but higher order multipole fields, and therefore higher order multipole moments, which contain the information to distinguish this physical property. For instance, state 2 of Figure 3 portrays a vortex with right chirality while the mean dipole moment angles result from averaging over right and left chirality vortices. By analyzing the quadrupole moment tensor components of state 2 (see Table 5) it can be noticed that for vortices of different chirality the strongest moments and exhibit distinct magnitudes with small standard errors. The mean quadrupole tensor components of state 2 are shown in Table 5 and were obtained by, first, normalizing the quadrupole moments of all the simulations associated to this state, which have similar dipole moments. Second, by classifying these vortex states by left and right chirality using the value of the largest normalized quadrupole tensor component. Finally, from these two groups, mean quadrupole tensor values were computed. These results are summarized in Table 5, showing that higher order moments provide additional information about the inverted configurations. In the context of state 2, the difference in the quadrupole moments originate from the corresponding quadrupole field, which can be observed from the residual of a dipole order inversion. A visualization of this residual is shown in Figure S6 in Supporting Information S1. Nevertheless, the physical meaning of the higher order multipole components cannot be interpreted directly as for the case of the dipole moments and their values depend on the center of expansion of the magnetic potential. Furthermore, for these spherical samples, the magnetic signal at the scan grid did not exhibit an evident multipolar character. Hence, the dipolar expansion produces a good estimation of inverted magnetic parameters, which can be seen from Figures 4c and 4d where the magnetizations given by the dipole expansion curve converges toward the higher order expansion curves.

Table 5

Quadrupole Moment Tensor Components of Vortices With Opposite Chirality

State 2 (vortex) in Figure 3
Quadrupole moments (N)	Chirality
	Left	Right
Norm × 10−26	2.1233 ± 0.0050	2.1303 ± 0.0080
Θ1(2)	0.0168 ± 0.0107	−0.0276 ± 0.0150
Θ2(2)	−0.0035 ± 0.0062	−0.0065 ± 0.0075
Θ3(2)	−0.6973 ± 0.0008	0.7425 ± 0.0169
Θ4(2)	0.0041 ± 0.0061	0.0115 ± 0.0099
Θ5(2)	−0.7152 ± 0.0009	0.6667 ± 0.0197

Note. The table shows the mean value of quadrupole moments obtained from the octupole inversion at a scan height of 1.0  μm, for single vortex domains of opposite chirality defined by state 2 of Figure 3. In the snapshot of the figure, it is shown a vortex with right chirality. Quadrupole moments are normalized (N) and then averaged. Left and right chirality states are averaged from 11 and 4 simulation results, respectively.

Multi Domain States and Higher Order Multipole Inversions

The potential of this micromagnetic modeling workflow in more realistic and complex magnetite grain geometries was tested by analyzing particles in synthetic and natural samples. In these cases, the geometric symmetry (e.g., spherical) is broken, thus the number of LEM is further restricted. The size of these samples is equivalent to spheres of diameter above 0.2 μm thus they are computationally more demanding to model. According to this, a few LEM states were selected for analysis. Because of the large dimensions of these structures, stray field effects are stronger and therefore vortex configurations predominate as the minima of energy in order to reduce the demagnetizing field energy. These magnetic configurations produce clear multipolar magnetic signals and it was shown that an octupole expansion efficiently allows to retrieve the dipolar moments and their directions, which is substantiated by the low residuals of the numerical inversion. In this context, dipole order multipole inversions produce noticeable residuals and hence overestimation of the inverted magnetization (see Figures 5g–5i), which is in contrast to the sphere samples where dipole inversions were more accurate. For both the natural and synthetic samples, the inverted total magnetizations and the orientation angles are unique for each magnetic state and therefore, they could be successfully inferred from their magnetic signal.

For the three cases of magnetic grains investigated here the numerical inversions were performed for scan surfaces located at increasing heights up to approximately 1 μm above the grains. In every situation, it was possible to successfully recover the magnetic moments due to a clear modeled stray field signal. In most cases, for sufficiently large scan heights a dipolar inversion gives a reasonable estimate of the grain magnetization because higher order moments decay more rapidly and have a smaller influence on the particle signal. Nonetheless, there are complex magnetic configurations, for example, Figure 5c, where higher order moments are crucial to accurately estimate the grain magnetization for scan heights in the studied range (see Figure 5i) since the magnetization solution from a dipolar inversion converges slowly toward the one produced by a higher order spherical harmonic expansion. This is particularly relevant for shallow grains which usually exhibit more complex magnetic signals (Cortés‐Ortuño et al., 2021) that can be detected via high‐resolution microscopy with a sensor adequately close to the sample surface. Furthermore, in practice, measurement noise can obfuscate the grain signal and it will be necessary to estimate the reliability of the inverted magnetization using a parameter such as the signal‐to‐noise ratio introduced by Cortés‐Ortuño et al. (2021).

Experimental Challenges: Tomography and Magnetometry

The study by Nikolaisen et al. (2020) demonstrated that micromagnetic modeling of individual grains in combination with FIB‐nt is able to find and categorize domain states in magnetite inclusions of rock samples and provides an important methodology for the accurate classification of stable paleomagnetic recorders. However, the mentioned study did not include any magnetic measurement and hence, did not recover the natural remanent magnetization of the magnetic particles. In order to retrieve the in‐situ magnetic configurations of individual grains, it is required to obtain not only position and shape of the particles but also sufficient magnetic information to identify the magnetization state via numerical inversion. Based on our pilot study this may be achieved by a combination of shape microCT data and magnetic surface scanning, if both methods (a) identify the grain shape and position with sufficient accuracy to define the three‐dimensional geometry for the micromagnetic model, and (b) provide sufficiently dense and accurate magnetic field data to identify a possibly unique micromagnetic LEM state that explains this surface measurement.

In the present work, an accurate modeling of the three‐dimensional profiles of the grains is highly determined by the tomographic data. Current technology may partly fall short of reaching sufficient tomographic resolution, in particular for grains in a size range below a few microns, which can be modeled with the current computational capacities. Although commercial microCT instruments do not yet allow for a routine scanning of nm voxel densities, recent advances in tomographic techniques (Lutter et al., 2021; Zabler et al., 2021), such as X‐ray nano tomography (nanoCT), have improved the maximal imaging resolution below the micrometer scale, and hence it might be possible to empirically test this investigation in sufficiently small grains. Furthermore, FIB‐nt already offers the capacity to profile a grain geometry with a few nanometers size resolution (Lutter et al., 2021; Nikolaisen et al., 2020), although it is a destructive technique that requires to perform the magnetometry scanning before applying the tomographic measurement.

In the case of magnetometry scans, recent studies applying the MMT method have used nitrogen‐vacancy magnetometry via a quantum diamond microscope with a spatial resolution of 1.2 μm (de Groot et al., 2021). Although this resolution seems large compared to the current tomography capabilities, it has been sufficient to successfully recover the magnetic moments of individual grains. By improving both the resolution and field sensitivity of the microscope it might be possible to reduce the effects of noise and detect weak signals from complex multi‐domain configurations. In this work, the scan surfaces have been modeled using a 1.5 μm side square region. Recent advances in nitrogen‐vacancy magnetometry indicate that spatial resolutions down to 50 nm have been achieved using scanning QDM (Scholten et al., 2021), which is promising for future experimental testing of a single grain signal. For larger samples with multiple grains, it has been suggested that widefield QDM can be an optimal method to measure the stray field signal. This technology is more simple to implement (no moving parts) and currently achieves spatial resolutions of 400 nm, which is restricted by the optical diffraction limit (Scholten et al., 2021).

Material Parameters and Anisotropy

What remains are substantial challenges to find out further crystallographic or mineralogical information to realistically constrain the micromagnetic modeling. A limitation of the microCT technique is that it cannot be used to directly measure crystallographic orientations beyond estimating it from the particle shape, for example, if it is an octahedral magnetite particle. In this context, Electron backscattered diffraction data may help to directly or indirectly determine the spacial orientation of the magnetocrystalline anisotropy axes (Nikolaisen et al., 2022) which will help to constrain the micromagnetic model.

An alternative method is to invert the crystallographic directions from QDM measurements. This could be done by measuring multiple different magnetization states, for example, after saturating the sample in different directions. From the anomaly as a function of field direction, one may be able to detect the crystallographic orientation. This general approach would allow to write a specific inversion micromagnetic code that also adapts orientation parameters and possibly other material (e.g., magnetostriction) or geometry parameters into the minimization algorithm.

Additional micromagnetic energy contributions, for example, from external stress fields (Béguin & Fabian, 2021) or inhomogeneous mineralogy currently are not taken into account, but may be inferred from the inversion itself. This can provide a very powerful method for obtaining such information on the particle level, especially if magnetic scans from several LEM magnetization states of the same particle—maybe even at different temperatures—are available.

Conclusions

The technique of MMT is a novel tool to study rock magnetic particle ensembles (de Groot et al., 2021). By combining microCT and high‐resolution scanning magnetometry, MMT allows to reconstruct the magnetic potential of individual well‐separated particles. However, for fundamental mathematical reasons, it cannot infer the internal magnetization structure of the measured particles (Fabian & de Groot, 2019). The spatial data of magnetic grains acquired for this technique can be used to effectively construct three‐dimensional numerical models of the individual particles.

Here, we showed that if the material constants of the particles are known, micromagnetic modeling can be used to substantially constrain, or even to uniquely infer, their internal magnetization structure. The method depends on the fact that in sufficiently small particles only a limited, even though in principle infinite, number of LEM states exist, which have distinguishable, and possibly unique, external potential fields. In the case of the magnetite sphere analyzed this study, it was shown that it can host a finite number of single‐domain vortex and uniform configurations. Moreover, their magnetization can be well described using a dipole order inversion, producing accurate and unique values of both the dipole moment intensities and directions. This was also observed for the single vortex configurations of the natural sample grain although at sufficiently high scan surfaces. In case that the easy crystallographic axes of the cubic anisotropy axes are not known, the observation of several single‐domain states might allow to determine their orientation.

More complex LEM magnetization states, such as the multi‐vortex states observed in the synthetic sample grain, exhibit strong quadrupole and octupole field signals and thus, require the reliable reconstruction of higher multipole moments for unique identification. For this work, the higher order multipole inversions showed more accurate results up to distances around a 1 micron order, of the scan surface from the grain top surface.

In this study, the modeled grain signals were introduced without noise and therefore allowed to efficiently recover magnetic moments for varying distances between measurement surfaces surface and particle up to values around 1 μm. Nevertheless, experimentally, the instrumental and geometric signal‐to‐noise ratio can limit an accurate reconstruction of the grain's magnetization, which can be verified in future studies using, for instance, high‐resolution QDM.

The results of this investigation provide the basis for a new method in rock magnetism to infer details of the magnetization structure of a large number of particles as a function of grain size and shape, and to experimentally study transitions between LEMs through repeat measurements during different demagnetization or magnetization procedures. This will allow to perform stability analysis of magnetic rocks at the individual particle level and, hence, accurately select the best paleomagnetic data recorders.

Supporting Information S1

Click here for additional data file.