Edinburgh Research Explorer The magnetic structure and palaeomagnetic recording fidelity of sub-micron greigite (Fe3S4)

We present the results of a ﬁnite-element micromagnetic model of 30 nm to 300 nm greigite (Fe 3 S 4 ) grains with a variety of equant morphologies. This grain size range covers the magnetic single-domain (SD) to pseudo single-domain (PSD) transition, and possibly also the PSD to multi-domain (MD) transi- tion. The SD–PSD threshold d 0 is determined to be 50 nm ≤ d 0 ≤ 56 nm depending on grain shape. The nudged elastic-band method was used to determine the room temperature energy barriers between sta- ble states and thus the blocking volumes. It is found that, in the absence of interparticle magnetostatic interactions, the magnetisation of equant SD greigite is not stable on a geological scale and only PSD grains ≥ 70 nm can be expected to carry a stable magnetisation over billion-year timescales, i.e., all noninteracting SD particles are essentially superparamagnetic. We further identify a mechanism for the PSD to multi-domain (MD) transition, which is of a continuous nature from PSD nucleation up to 300 nm, when structures typical of MD behaviour like closure domains begin to form. © 2017 The Authors.


Introduction
The ferrimagnetic mineral greigite (Fe 3 S 4 ) is the iron sulphide analogue of the iron oxide magnetite (Fe 3 O 4 ). It is commonly formed as a precursor to pyrite in early diagenetic anoxic sulphatereducing sediments (Berner, 1984;Hunger and Benning, 2007) and as a product of biomineralisation by magnetotactic bacteria (Mann et al., 1990). Although thought to be thermodynamically unstable under most sedimentary regimes, it has been found to be stable on geological timescales (Roberts et al., 2011), making greigite a possible natural remanent magnetisation (NRM) carrier in many systems, e.g., lacustrine (Babinszki et al., 2007;Ron et al., 2007) and marine (Roberts and Turner, 1993;Roberts and Weaver, 2005;Rowan and Roberts, 2006;Rowan et al., 2009) sediments, oil fields' shallow overburdens (Abubakar et al., 2015;Donovan et al., 1984;Reynolds et al., 1993) and gas-hydrate-bearing sediments (Larrasoaña et al., 2007). To further our understanding of the potential use of greigite as a proxy for environmental change, hydrocarbon exploration, magnetostratigraphy and in general the contri-* Corresponding author.  bution of this iron sulphide to the magnetic properties of rocks, we have implemented micromagnetic numerical finite-element method (FEM) simulations of greigite. To characterise its basic properties we have modelled its magnetic domain state's shape and size dependence using the DUNLOP package (Nagy, 2016) and its stability on geological timescales at room temperature using the MERRILL package (Nagy et al., 2017).
The model solutions are dependent on a balance between various magnetic forces and thus it is important that the material's magnetic parameters be known as accurately as possible. Past difficulties in producing pure greigite samples on which to determine these parameters has resulted in a lack of accurate models in the literature. However, these difficulties have been overcome, and recent measurements on highly pure, high crystallinity synthetic greigite samples (Chang et al., 2009(Chang et al., , 2008Li et al., 2014;Winklhofer et al., 2014) now allow for numerical micromagnetic models to probe the magnetic microstructure of greigite.
The fundamental magnetic parameters of greigite used in this investigation are the saturation magnetisation M S = 3.51 μ B p.c.u. (Li et al., 2014) or ∼2.7 × 10 5 A/m which is ∼11% higher than the value reported by Chang et al. (2009) of 3.25 μ B p.c.u. (and ∼57% the value of M S for magnetite) and the exchange stiffness constant A = 2 × 10 −12 J/m (Chang et al., 2008)  for magnetite). Winklhofer et al. (2014) estimated a cubic magnetocrystalline anisotropy (MCA) term K 1 = −1.7 × 10 4 J/m 3 (∼42% higher than the value for magnetite) and negligible second MCA constant K 2 to K 1 ratio, i.e., the easy axes are the <111>.
Very few micromagnetic models of greigite have so far been attempted, with only the work of Muxworthy et al. (2013) having been published. They implemented a micromagnetic finitedifference (FD) method that used the earlier value of M S from Chang et al. (2009) to model both individual grains and the effects of magnetostatic interactions between cuboidal particles of greigite arranged in chains. For the single crystals they found good agreement with the analytical calculations of Diaz-Ricci and Kirschvink (1992) (for which crude estimates of the magnetic parameters were used). However, due to their structured spatial discretisation, FD methods are not as well suited for the simulation of the euhedral morphologies seen in natural samples (Snowball, 1997) as FEMs are. Elongation is a common feature of magnetosomal grains but not of particles with non-biogenic origin. Although Muxworthy et al. (2013) have demonstrated that particle elongation and interparticle magnetostatic interactions are important, this study is limited to a variety of equant, isolated particles whose behaviour represents the limit as the particle concentration approaches zero. As such, this study provides a first step towards a more complete picture of the role of non-biogenic greigite in rock magnetics.
The grain size of a magnetic mineral strongly affects its magnetic behaviour and palaeomagnetic recording fidelity. As grain size distributions are usually inferred from bulk magnetic properties, a deep understanding of how a grain's magnetic properties depend on size is needed. We present here the results of a numerical FEM study of the stress-free, zero-field domain states of spherical and equant euhedral single grains of greigite in the single-domain (SD) to pseudo single-domain (PSD) range. To address this question we used the nudged elastic-band (NEB) method  to calculate action-minimising paths (AMPs) between stable magnetic configurations which allow us to determine energy barriers and from these the blocking volumes of naturally occurring equant euhedral grains of greigite.

The micromagnetic algorithm
For a ferromagnetic (in the broad sense, i.e., ferro-or ferri-) material-neglecting external fields, thermal and magnetostrictive effects-the Gibbs free-energy functional can be written as (Brown, 1963) with the ferromagnetic volume. Here, with the reduced (unitary) magnetisation vector m, is the expression for the continuum approximation of the energy density due to the quantum-mechanical exchange forces (Landau and Lifshitz, 1935).
with γ i the direction cosines, is the MCA energy density in the cubic anisotropy system. In the Cartesian system this has the form: where K 2 has been neglected since in the cubic anisotropy system we are assuming, K 1 is the dominant term at room temperature. Finally, is the magnetostatic self-energy density, with H stray the stray field produced by the ferromagnetic body. It is known from thermodynamics that under isothermal and isobaric conditions a system will spontaneously evolve towards an equilibrium state with minimal Gibbs free-energy. It is the aim of a micromagnetic algorithm to find the equilibrium magnetisation m.
The variational principle proposed by Brown (1963) and Therefore, in equilibrium the magnetisation is parallel to an effec- is the local torque produced on the magnetisation by the effective field at each magnetic site, so we have equilibrium via the vanishing of the torque (Brown, 1963). We thus have the two main approaches towards finding a micromagnetic solution, i.e., the minimisation of the Gibbs freeenergy functional, which is faster though prone to false minima, and the vanishing of the torque, which is slower but more robust and physically meaningful (Gilbert, 2004). The former is usually achieved with a variety of gradient descent methods (Fischbacher et al., 2017) while the latter by numerical integration of the Landau-Lifshitz-Gilbert (LLG) equation (Gilbert, 2004) ∂m ∂t where α is a phenomenological damping constant, and γ = γ /(1 + α 2 ) with γ the gyromagnetic ratio.
FEMs allow for an unstructured discretisation of the spatial domain which in our case is decomposed into tetrahedral elements. In the treatment of the micromagnetic theory of Brown (1963) there are some linearisations, which means that there should not be large variations in the direction of m between neighbouring nodes in the FE mesh. To model nonuniform structures it is sufficient that the spatial discretisation in the model be smaller than the exchange length l exch. = 2 A/μ 0 M 2 S (Rave et al., 1998), which for greigite is l exch. ≈ 6.6 nm; a maximum element size of 5 nm has been chosen for all the models. The non-local problem of calculating the stray field is handled via a transformation method (Imhoff et al., 1990) in DUNLOP and by a hybrid finite-element/boundaryelement formulation (Fredkin and Koehler, 1990) in MERRILL.

Choice of morphologies
Based on scanning electron microscopy (SEM) and transmission electron microscopy (TEM) micrographs of synthetic greigite samples (Chang et al., 2008;Li et al., 2014) and of naturally occurring samples (Snowball, 1997;Vasiliev et al., 2008), five octahedral shapes have been chosen with increasing degrees of truncation of their corners: from no truncation at all (octahedron) to a 'completely' truncated shape (cuboctahedron) and three truncated octahedral shapes in-between to which we further refer to as minimally truncated octahedron, truncated octahedron (the regular case), and maximally truncated octahedron (Fig. 1). This series of shapes models the influence of the ratio of {001} to {111} faces, which increases with truncation. Furthermore, spherical shapes which, although unrealistic, serve as 'control subjects' that do not exhibit configurational anisotropy  and cubic shapes which have been modelled before by Muxworthy et al. (2013) and represent the case of only {001} faces. All the volumes are normalised to cubes, i.e., a particle sized 120 nm has a volume of (120 nm) 3 .

Crystal growth model
The magnetic domain structure of a ferromagnetic nanoparticle obtained by a micromagnetic algorithm is not only a function of its mineralogy, size and shape, but also of its history. In particular, it is known that a SD particle can grow and remain SD up to a threshold size d max after which it will become PSD (Enkin and Williams, 1994). (In this study we are concerned with the zero-field microstructure and properties; in this context, the onset of PSD behaviour is marked by the formation of a single-vortex structure.) If the particle then decreases its volume it will remain PSD down to a threshold d min (< d max ) below which it will become SD again. This defines a size range in which a ferromagnetic grain can be either SD or PSD dependent on its history (Muxworthy and . This phenomenon has been modelled for the seven morphologies ( Fig. 1) in the 30 nm to 300 nm size range.
Starting from a 30 nm particle, we obtain the micromagnetic solution and extrapolate it to a larger grain. This becomes the new initial condition for which we solve and repeat, growing the particle in steps of 2 nm. Since a very fine incremental size step is used, much smaller than the exchange length, we can be certain that we are not missing any features from one step to the next. This process accurately models grain growth.
Once the particles have grown to 120 nm, the procedure is then reversed in decreasing steps of 2 nm until the initial size of 30 nm is reached. Since chemical alteration usually proceeds by alteration of the surfaces, the volume decreasing process can be thought of as a model for chemical alteration to a non-magnetic phase. This growth from 30 nm to 120 nm followed by the reverse process is referred to as the main loop (ML).
Solutions on the size-descending curve not found on the sizeincreasing curve were also subject to growth; the secondary loop (SL). The ML and SL allow us to investigate the different domain states with size, shape and history. These micromagnetic solutions were performed by numerical integration of the LLG equation (Eqn. (9)). Stable solutions at 120 nm, whether found on the ML or SL, were further grown up to 300 nm. Energy minimisation was used for calculations larger than 120 nm, as integration of the LLG equation is prohibitively slow at these sizes. The parallel DUNLOP package (Nagy, 2016) was used to model these grain size dependencies.
These models overlook the effect of thermal fluctuations. However, this limitation is addressed by calculation of the relaxation times. These allow the determination of the particles domain state at a given size beyond the capabilities of standard micromagnetics.

Relaxation times
Over a sufficiently long observation time, termed the relaxation time, a ferromagnetic particle will switch between different stable states due to thermal activation. The relaxation time is given by (Néel, 1955) where K B is Boltzmann's constant, T the temperature at which the transition occurs and τ 0 the attempt time, commonly with a value of ∼10 −9 s (McNab et al., 1968). Any transition between stable states must occur along an AMP ; E B is the energy barrier along such a path. When E B is of the order of the thermal energy available on a timescale of interest, a particle is in a superparamagnetic (SP) state, spontaneously switching back and forth between different SD orientations. SD grains with an energy barrier larger than ∼60 K B T have a relaxation time in the order of billions of years and are thus considered stable SD (SSD) reliable palaeomagnetic recorders (Dunlop and Özdemir, 1997). The SP-SSD threshold is therefore one of the key parameters to assess the palaeomagnetic significance of a NRM. A NEB method  has been implemented in the micromagnetic package MERRILL (Nagy et al., 2017) (such methods are currently unavailable in DUNLOP) to calculate the energy barriers between the stable configurations for the (truncated) octahedral particles, from 30 nm increasing in steps of 2 nm until the blocking volume is reached, which is here taken to be the volume for which the relaxation time surpasses four billion years.
Unlike early NEB methods (as applied to magnetic systems) which approximated each particle as a single dipole (Berkov, 1998), here we are concerned with full micromagnetic solutions and thus, our results extend beyond SD configurations and coherent rotations. A similar method was implemented by Dittrich et al. (2002) to study the energy barriers in patterned granular media.

Magnetic structure in the SD-PSD regime
Starting at 30 nm from a randomised initial condition, all shapes relax to a SD state along an easy direction. Extrapolating onto larger grains, this state remains stable up to a threshold d max (Figs. 2a, c). When the SD state solutions are grown beyond d max , numerically stable up to d max = 92 nm (c). Growing this solution to a 94 nm grain it is found to relax to a [111] EAV (d), stable up to 120 nm. The EAV is then interpolated into smaller grains, stable down to d EH = 54 nm (a). At 52 nm the EAV goes to a [001] HAV (e), stable down to d min = 40 nm. At 38 nm, the solution relaxes to the original SD state (a). This sequence is referred to as a Type 1 main loop (ML). Growth of the HAV from 52 nm forms the Type A secondary loop (SL): the HAV is found to be stable up to d HE = 68 nm (a) and to realign with the easy direction at 70 nm. Vortex states can not only be easy or hard-aligned, but also <011> intermediate-aligned vortices (IAVs) (f) and distorted IAV configurations (dIAV) (g). MLs in which IAVs are found are referred to as Type 2. Colour represents the MCA energy normalised by |K 1 |. The vortex cores are highlighted by obtaining a helicity (K = m · ∇ × m) isosurface and reducing the opacity of the rest of the arrows. they are found to relax to an easy aligned vortex (EAV) (Figs. 2a,d). For all shapes, this vortex configuration is stable up to 120 nm. On reversal, the EAV is stable down to a threshold d EH (Fig. 2a). Below d EH , the EAV goes to a hard-aligned vortex (HAV) preserving its chirality (Figs. 2a,e). Further decreasing the volume, the HAV is stable down to d min and below that, the HAV relaxes back to SD (Fig. 2a). This general behaviour on the ML is referred to as Type 1 (Figs. 2a, 3a, e, g, k), for which we introduce the notation: with 30 nm < d min < d EH < d max < 120 nm.
Another sequence was observed. On the size-descending curve, the EAV is stable down to a threshold d EI below which the EAV Reduced magnetisation (a, c, e, g, i, k) and energy density (b, d, f, h, j, l) against size. All shapes relax to an easy aligned SD state at 30 nm from a randomised initial condition. The SD state is numerically stable up to d max (black lines). Growth beyond d max results in the magnetisation relaxing to an EAV, stable up to 120 nm (red lines). The solution is then interpolated into smaller grains. This forms the main loop (ML) (opaque lines). The EAV is stable down to a threshold beyond which a HAV is nucleated on the size-descending curve. The HAV is stable down to d min , below which it relaxes to a SD state. This is a Type 1 ML (a, e, g, k). When the EAV goes through an IAV (green lines) before going to the HAV configuration the ML is Type 2 (c, i). Growth of the HAVs and IAVs found on the size-descending curve forms the secondary loop (SL) (translucent lines). When these realign with the EAV or with the vortex state they nucleated from the SL is Type A (c, e, g, i (translucent blue line)). When they are stable up to 120 nm the SL is Type B (a, i (translucent green line), k). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) stable down to d IH . Below that, the IAV goes to a HAV which remains down to d min . Finally, below d min the HAV relaxes back to SD. This ML behaviour is referred to as Type 2 (Figs. 3c, i), with the formula: Growth of the vortex configurations found on the sizedescending curve (HAVs, IAVs) forms the SL. When the ML is Type 1, the HAV is grown up to a threshold d HE beyond which it realigns with an easy direction (Fig. 2a). When the ML is Type 2, the HAV goes to either an EAV (Fig. 3c) or to the IAV it nucleated from (Fig. 3i). When the IAV is grown beyond a threshold d IE it realigns with an easy direction (Fig. 3c). This behaviour on the SL is referred to as Type A (Fig. 2a, 3c, e, g), with the general formula (in the case of growth of the HAV and realignment to an EAV): There is also the possibility that the HAV/IAV is stable up to 120 nm. This SL is a Type B (Fig. 3a, i, k), with the general formula (in the case of growth of the HAV): Fig. 4. Micromagnetic structures of octahedra, minimally truncated octahedra and regular truncated octahedra. Left column (a, d, g) shows the largest SD solutions at d max , obtained by interpolating from solutions for smaller grains starting at 30 nm. On interpolating to a grain 2 nm larger, the structure relaxes to an EAV (b, e, h), stable up to 120 nm. From 120 nm interpolation is carried out into smaller grains. Eventually the vortex aligns with a hard direction (c, f, i) which is stable down to d min after which the solution becomes SD again down to 30 nm. Top row (a, b, c) shows the structures for the octahedra; centre row (d, e, f) for the minimally truncated octahedra and bottom row (g, h, i) for the regular truncated octahedra. Colour represents the MCA energy normalised by |K 1 |.
3.1.1. Spheres Fig. 2a shows the reduced magnetisation and Fig. 2b the energy density (against size) for the spherical shapes. Spheres showed a Type 1 ML (Eq. (11)), with a specific formula: The SD state at d max = 92 nm and the EAV at 94 nm are shown in Figs. 2c-d. The SL is formed by growing the HAV found at 52 nm (Fig. 2e). This is a Type A SL (Eq. (13)), with formula: Fig. 3 shows the reduced magnetisation and energy density plots for the rest of the shapes. The octahedra (Figs. 3a-b)

Octahedra and truncated octahedra
The SD state at d max = 68 nm (Fig. 4a) shows significant flowering-deflection of the magnetisation onto edges and vertices. The EAV nucleated at 70 nm is shown in Fig. 4b. Growth of the HAV nucleated at 50 nm (Fig. 4c) showed a Type B SL (Eq. (14)), i.e., the HAV was found to be stable up to 120 nm, The SD state shows greater stability and less flowering at d max = 74 nm (Fig. 4d) than for the octahedra, before relaxing to an EAV at 76 nm (Fig. 4e). Growth of both the IAV from 48 nm (Fig. 2f) and HAV from 46 nm (Fig. 4f) forms the composite SL. Both showed a Type A SL (Eq. (13)), specifically: The next two degrees of truncation, the regular truncated octahedra (Figs. 3e-f) and the maximally truncated octahedra (Figs. 3g-h), both showed Type 1 MLs (Eq. (11) The SD states show increasing stability and less flowering (Figs. 4g, 5a). The EAVs at 80 nm and 82 nm are shown in Figs. 4h, 5b. Growth of the HAVs nucleated at 48 nm ( Fig. 4i) and at 46 nm (Fig. 5c) The SD state is more stable, with d max increasing to d max = 114 nm ( Fig. 5d) before relaxing to an EAV at 116 nm (Fig. 5e). Growth of the HAV from 60 nm (Fig. 5f) shows a Type A SL (Eq. (13)), with the HAV going to the dIAV it nucleated from. The dIAV (Fig. 2f

Cubes
The cubes (Figs. 3k-l) showed a Type 1 ML (Eq. (11)). The EAV nucleated from the SD state is a distorted EAV (dEAV). The dEAV mostly aligns with an easy direction, but its ends deflect from the vertices to a hard direction, attaching to opposite faces. The ML is then: The SD state at d max = 80 nm is largely flowered (Fig. 5g). The dEAV initially nucleated at 82 nm (Fig. 5h) could also be interpreted as a distorted HAV. However, the alignment of such state with an easy direction becomes more obvious for larger cubes. Growth of the HAV from 62 nm (Fig. 5i) results in a Type B SL:

Discussion
The values for d max of the spheres are significantly larger than those found for the euhedral shapes (except for the cuboctahedra which have an anomalously large value), perhaps due to the absence of corners to act as nucleation points. Truncation of the octahedral particles increases the (numerical) stability of the SD solutions which is expressed as the increase in d max from 68 nm for the octahedra (Fig. 3a) to 114 nm for the cuboctahedra (Fig. 3i) ( Table 1). A less pronounced effect is the decrease of d min with truncation from 46 nm for the octahedra to 38 nm for the cuboctahedra (Table 1). HAVs are more (numerically) stable at smaller sizes the more truncated the particle is. This is due to the large stray field energy created by a vortex pointing towards a vertex. It is energetically favourable for a vortex to attach its ends to flat surfaces large enough to accommodate its core. This is seen in the distortion of the vortex structures shown in Figs. 2g, 5h. These avoid the production of large stray fields by attaching their ends to grain faces, at the cost of anisotropy and exchange energies needed to distort the otherwise straight structures of the vortices.
The energy plots (Figs. 2b,3b,d,f,h,j,l) show the SD energy density is fairly constant with size for all shapes. For the octahedra, the intersection of the EAV and HAV energy curves occurs above the SD curve. This means that it is then the EAV energy curve which first intersects the SD curve at d 0 and thus this PSD state becomes the GEM thereon (Fig. 3b). With truncation, this intersection moves closer to the SD curve as can be seen in Fig. 3d in which all the different PSD states (EAV, IAV and HAV) and the SD energy curves meet at roughly the same point. Further truncation causes this intersection to eventually occur below the SD energy curve. This creates a narrow range of sizes for which the HAV is the lowest energy state (Figs. 3f, h). Completely truncated, the cuboctahedra show a split of this intersection into distinct crossings of the HAV/IAV and IAV/EAV energy curves (Fig. 3i, compare with Fig. 3d), creating a broad range of sizes, from ∼50 nm to ∼66 nm, for which the HAV has the lowest energy. A range for which the HAV has the lowest energy was also found for spheres and cubes. The overall effect of truncation on d 0 is to decrease this threshold (Table 1).
For the cubes we found d max = 80 nm, much smaller than the value of 107 nm obtained by Muxworthy et al. (2013). We found d min = 38 nm, in agreement with the value by Muxworthy et al.
(2013). We found the intersection of the SD and HAV energy curves d 0 ≈ 52 nm (Fig. 3l), lower than the value of 58 nm by Muxworthy et al. (2013). Modelling the ML for cubes with the M S value used by Muxworthy et al. (2013) we obtained d max = 92 nm, d min = 42 nm smaller and larger, respectively, than the values by Muxworthy et al. (2013). The differences in d max , d min can be due to Muxworthy et al. (2013) using a FD method as opposed to a FEM used here. However, excellent agreement of d 0 ≈ 60 nm was found with the value of 58 nm by Muxworthy et al. (2013). The difference between d 0 obtained with the different M S is significant as it is larger than the exchange length and thus, unlikely to be an effect of discretisation.

Identifying the PSD-MD transition
Unlike fine SD grains, bulk ferromagnetic materials can possess a (roughly) null net magnetisation. This is because bulk ferromagnets, in their lowest energy state, have a multi-domain (MD) magnetic structure (Dunlop and Özdemir, 1997). MD structure is characterised by the coexistence of multiple magnetic domains: small regions saturated in different easy directions and separated by narrow planar regions called domain walls where the magnetisation vector transforms continuously from the direction of one domain to that of its neighbour. Near the surface of a bulk sample the magnetic domains lie in directions parallel to the surfaces, which do not necessarily coincide with the magnetocrystalline easy axes; these are called closure domains as they enclose the magnetic flux. Then, MD grains minimise the stray field energy at the expense of exchange and MCA energies associated with the domain walls and closure domains. This balance becomes more energetically favourable as a particle grows. To identify a plausible PSD-MD transition, EAV states found to be stable at 120 nm were further grown up to 300 nm in steps of 2 nm.

Growth of EAVs
The octahedron EAV (Fig. 4b) was found to remain stable up to 300 nm-as was the case for the EAVs for all shapes. However, although the overall structure is preserved as the grain grows, there is a gradual, continuous transformation towards MD structure. As the crystal becomes larger, the vortex core (the cylindrical region encompassing the highest helicity K = m · ∇ × m) width remains the same. The regions radially far from the core become screened from its influence and more influenced by the effects of MCA, giving way to six ever larger expanses of the grain that become magnetised along easy axes. The sections with higher MCA energy become increasingly narrower, more planar and confined between the easy-aligned regions. With further growth, the regions close to the edges become magnetised along the edges thus completing a picture of MD structure: magnetic domains magnetised along easy directions (Fig. 6a); flat, narrow magnetic domain walls dividing these and closure domains formed along edges of the body (Fig. 6b). The spheres (not shown) show this same general behaviour of formation of magnetic domains and domain walls, but without the formation of closure domains as there are no edges to nucleate these type of domains.
At 300 nm, the regular truncated octahedron shows the same basic structure as the 300 nm octahedron except for the widening of the domain walls as they approach the square {001} faces . Completely truncated, the cuboctahedra have a somewhat different structure. In the 300 nm cuboctahedron solution the domain walls have become so wide as they draw closer to the square faces that they engulf three of the magnetic domains (Fig. 6e), and hard-aligned Néel walls appear which are visible as blue streaks on the square faces (Fig. 6f).

Discussion
A tentative mechanism for a PSD-MD transition has been identified. This proceeds by a gradual formation of magnetic domains, surfaces) is to widen the domain walls as they approach the surfaces as this reduces the stray field energy (d). Closure domains along the edges are still pronounced (d). Fully truncated, the effect on the structure of the cuboctahedron EAV grown to 300 nm is to reduce the proportion of the volume occupied by magnetic domains (e). The domain walls are so wide close to the surfaces that they engulf three of the domains, while also hard-aligned Néel walls are formed as seen from the blue streaks on the surface (f). Colour represents the MCA energy normalised by |K 1 |. (For interpretation of the references to colour in this figure, the reader is referred to the web version of this article.) domain walls and closure domains 'seeded' by the vortex structures nucleated in the SD-PSD transition regime. However, this transition is more obvious for the octahedra than the truncated octahedra and cuboctahedra.
Since our models do not account for the effects of thermal fluctuations, the stability of the solutions is only numerical, e.g., the energy landscape of a micromagnetic solution is somewhat flat in the vicinity of the solution, thus small numerical perturbations can be insufficient for a micromagnetic algorithm to drive the solution to a new local energy/torque minimum depending on the sensitivity or control parameters of the algorithm. We find for the octahedra and the cuboctahedra two solutions for each morphology up to 300 nm (Figs. 3a, e). Although the EAVs for both have the lowest energies, this leaves open the question of whether we can expect to find metastable grains with the higher energy structures. Likewise, in Section 3.1 we find that grains remain SD beyond d 0 . In theory, a particle can remain in a metastable SD state beyond this threshold if the energy barrier separating the SD from the PSD state is higher than the thermal energy available. Knowledge of the energy barriers near the SD-PSD transition is needed to answer these questions.

Energy barriers and blocking volumes
A nudged elastic-band (NEB) method  was implemented for the calculation of the energy barriers at room temperature as a function of size and shape for the (truncated) octahedral morphologies. To calculate the energy barrier for a given shape and size we obtain many solutions from randomised initial conditions. From this set of solutions the state with the lowest energy is identified. For the smaller particles we expect the GEM to be SD and for larger grains PSD. Once the GEM has been identified from the set of initial solutions, two appropriate solutions must be chosen to calculate the action-minimising path (AMP) between them and thus the energy barrier.
When the GEM is a <111>-aligned SD state, the pair of appropriate solutions are magnetised normal to contiguous {111} faces.
Above d 0 the GEM is usually a HAV, then, the paths to test are Fig. 7. Typical action-minimising paths. The regular truncated octahedra illustrate the typical AMPs found for all the (truncated) octahedral shapes. Well below the SD-PSD threshold the AMP between SD configurations is a coherent rotation (left column). The energy along the AMP is a single bump (a) and the energy barrier is an intermediatealigned SD state (d). At sizes just below the SD-PSD threshold the AMP is still between SD states, but via a curling mode (vortex nucleation) (centre column). Starting from SD (e), a vortex is nucleated (h) then unwinds to an intermediate-aligned SD state which has the highest energy along the AMP (not shown). Then, the final SD state (k) is reached via coherent rotation. This makes the energy along the AMP more complex and asymmetric (b). Well above the SD-PSD threshold, once the EAV has the lowest energy, the AMP is via a structured rotation of the vortex core (right column). The transition is between two EAVs with the same chirality (f, l) as a change of chirality produces prohibitively large energy barriers. In the AMP there is an IAV sitting in its own shallow LEM (i). The energy barrier is the one that separates the EAVs from the IAV. However, as the grain grows, the IAV LEM becomes more shallow and the energy along the AMP flattens (c). Colour represents the MCA energy normalised by |K 1 |. The energies along the AMPs are plotted in units of K B T , with T = 300 K. The left column of Video 1 (supplementary material) shows these transitions. between vortices at 90 • from each other, e.g., between [100] and [001] vortices. For some shapes the GEM is a distorted vortex, for these, it is important to calculate several paths between different configurations to find the transition with the lowest energy barrier. At larger sizes, for all shapes the GEM is an EAV, this means that the paths to calculate are between vortices pointing towards contiguous {111} faces, much like for SD particles. For all PSD to PSD transitions the vortices must have the same chirality as a change produces prohibitively large energy barriers one to two orders of magnitude larger than the AMP barrier and thus we can neglect the possibility of such transitions. It is not necessary to calculate the energy barriers for transitions other than the ones with the lowest energy as these dominate the behaviour and higher energy transitions usually proceed via lower energy ones (Nagy et al., 2017). For perfectly regular, equant grains as those modelled in this study, the lowest energy transition is degenerate and so, the relaxation time has to be divided by three (for the three distinct degenerate paths). Fig. 7 shows three examples of transitions for the regular truncated octahedra which were found to be typical for all the (trun-cated) octahedra, but for the cuboctahedra. At the smaller sizes (Fig. 7, left column) the transition is from SD to SD via coherent rotation. The energy along such a path is plotted in Fig. 7a; the energy barrier is the difference between the easy (Figs. 7d, j) and intermediate-aligned SD (Fig. 7g) energies. The energy barrier for such transitions, for all shapes, is quite low, with relaxation times from a few microseconds to a few days for the largest of these. A few nanometres before the SD-PSD transition the AMP is not a coherent rotation, but a transition via a curling mode (Fig. 7, centre column). The energy barrier is also the difference between the easy and intermediate-aligned SD energies.

Octahedra and truncated octahedra
Once the GEM is an EAV, the transitions are between isochiral EAVs directed towards contiguous {111} faces (Fig. 7, right column). The transition is a structured rotation of the vortex, through an IAV (Fig. 7i), which maintains its overall shape along the path.
The energy along such paths is a double-bump where the IAV sits in a shallow local energy minimum (LEM). This means that the energy barrier is not the difference between the EAV and IAV energies. Rather, it is given by the barrier separating the EAV from the IAV. However, with increasing size the IAV LEM becomes more shallow and thus the energy barrier is better approximated by the difference between the EAV and IAV energies-at 64 nm the IAV LEM is so shallow that the double-bump feature becomes flattened (Fig. 7c). In general, the PSD to PSD transitions are very similar to coherent rotations between SD states in that they are structured rotations of the vortex core. The left column of Video 1 (supplementary material) shows these typical transitions.

Cuboctahedra
The cuboctahedra show a very different behaviour. Their relaxation times increase exponentially for the SD to SD transitions, but do not drop for the first PSD to PSD transitions, at 48 nm, which are hard-aligned to hard-aligned (Fig. 8, left column) and pass through a dIAV (Fig. 8g). The energy barrier is an IAV (Fig. 8j). The relaxation times for these transitions then decrease as the HAV and dIAV energies get closer until, at 66 nm, the GEM is a dIAV. Then, since the straight IAV has a very high energy, transitions in which the dIAV ends reattach to contiguous faces are unfavourable. The transition with the lowest energy barrier is between two distorted vortices like the ones shown in Figs. 8e, k, which are (though distorted) [011]-and [110]-aligned. The transition preserves the shape of the distorted vortex as it rotates in a structured manner, keeping its ends attached to the square faces.
The relaxation times for this type of transitions plateau at a few microseconds up to 74 nm. In Fig. 3j the energies of the dIAV and the EAV intersect at roughly 90 nm and are very close up to 120 nm. Since the transition between two EAVs is likely to go through a dIAV, it is not expected that the relaxation time for EAV to EAV transitions will rapidly increase until sizes much larger than 90 nm. Indeed, we find that only for the particles larger than 110 nm the relaxation times start to grow exponentially-though at a rate slower than for (truncated) octahedra transitions-surpassing 4 billion years at roughly 134 nm (Fig. 9c). A typical transition of this type at 132 nm is shown on the right column of Fig. 8 between [111]-and [111]-aligned vortices via vortex distortion (not shown). The highest energy state along the path is a straight IAV (Fig. 8i). The right column of Video 1 (supplementary material) shows these transitions. Fig. 9 summarizes the results of the relaxation times obtained for the five (truncated) octahedral shapes in the 30 nm to 74 nm range and the relaxation times for the cuboctahedra in the 120 nm to 134 nm range (Fig. 9, inset). For all shapes, except the cuboctahedral, the general behaviour is an exponential increase of the relaxation times for the SD to SD transitions from a few tens of microseconds to a few seconds for the regular truncated octahedra and up to 12 days for the octahedra, followed by a sharp drop at 46 nm (for the regular truncated octahedra) to 50 nm (for the octahedra). This drop in the relaxation times marks the SD-PSD threshold. Then, the PSD to PSD transitions have increasingly lower relaxation times until reaching a global minimum at ∼50 nm to ∼54 nm. Once the GEM is the EAV the relaxation time shoots up exponentially with crystal growth, surpassing 4 billion years at ∼68 nm to ∼72 nm.

Discussion
These calculations show that equant, non-interacting, roomtemperature SD grains of greigite are on a geological timescale, SP given that the largest relaxation time found for SD particles (octahedra) of 12 days is only stable on laboratory timescales. Therefore, the expected SP-SSD threshold does not exist but rather a SP-PSD threshold since it is only the PSD grains which can hold a remanence for extended periods of time. Muxworthy et al. (2013) estimated the energy barriers for non-interacting cubes of greigite from the anisotropy energy surface and found a SP-SSD threshold of 56 nm, for a relaxation time of four billion years. In this study the blocking volumes for cubes and spheres were not calculated as these are not naturally occurring morphologies for greigite. However, given the results for the five (truncated) octahedral shapes, it is unlikely that a SP-SSD threshold should exist for either cubes or spheres.
Since the relaxation times for the larger SD particles are so small and since the transitions from SD to SD just before the GEM becomes PSD are through a curling mode (vortex nucleation), we conclude that the SD-PSD threshold is defined only by the first intersection of a PSD energy curve with the SD energy curve, d 0 (Figs. 2b,3b,d,f,h,j,l). That SD-SD transitions just below the SD-PSD threshold occur via a curling mode was also observed for magnetite by Enkin and Williams (1994) and shows that greigite particles are likely to show PSD aspects even below the SD-PSD threshold. It is highly unlikely that metastable SD grains exist beyond d 0 in the absence of grain-grain interactions. The relaxation times between PSD states just above d 0 are the shortest of all and only start to grow exponentially for the EAV to EAV transitions. We thus conclude that the PSD-MD transition through a crystal growth mechanism, must occur through an EAV path.
The splitting of the energy curves' intersections exhibited by the cuboctahedra not only has repercussions on the SD-PSD threshold, but also on its stability against thermal fluctuations. All the (truncated) octahedral particles have blocking volumes of ∼72 nm while for the cuboctahedra this volume is ∼134 nm. This could be explained by the fact that the first PSD to PSD transitions, between hard-aligned PSD states, traverse through IAVs (Fig. 8, left column): as the cuboctahedral particle grows, the energy difference between the IAV and HAV becomes smaller (Fig. 3j) and this is reflected in a diminishing of the relaxation times for this type of transition (Fig. 9). Since for the rest of the shapes the PSD and SD energy curves all meet in a relatively short span, this effect of diminished stability with increasing size only lasts for a much smaller range of sizes, and once the EAV is the GEM, the relaxation times for transitions between this type of domain state grow very rapidly. However, for the cuboctahedra the EAV becomes the GEM at a larger size than for the rest of the shapes and the energy difference with the dIAV, through which the AMP goes through, is very small.

Conclusion
We have presented calculations of the SD-PSD threshold for realistic shapes of equant greigite grains. We have found that octahedral shapes have the largest SD-PSD threshold d 0 ≈ 56 nm, a value lower than the threshold obtained for magnetite octahedra by Witt et al. (2005) of ∼73 nm. This value decreases by the Fig. 8. Action-minimising paths for the cuboctahedra. The cuboctahedra showed different AMPs above the SD-PSD threshold from the rest of the shapes. Above the SD-PSD threshold, and up to 66 nm, the lowest energy state for the cuboctahedron is a HAV. Transitions between HAVs (left column) occur via a distortion (g) and structured rotation (j) of the vortex core. The energy barrier is an IAV (j). The dIAV along the AMP (g) sits in its own LEM, forming a three-bump energy barrier (a). The dIAV becomes the lowest energy from 66 nm. Transitions between dIAVs (centre column) form a single bump energy barrier (b). The transition is a structured rotation of the distorted vortex, keeping attached to the same surface (e, h, k). Once the lowest energy is an EAV, transitions between these (right column) occur via a distortion and structured rotation of the vortex core. The energy barrier is an IAV (i). Colour represents the MCA energy normalised by |K 1 |. The energies along the AMPs are plotted in units of K B T , with T = 300 K. The right column of Video 1 (supplementary material) shows these transitions. effects of truncation to d 0 ≈ 53 nm for the regular truncated octahedra and further down to d 0 ≈ 50 nm for the cuboctahedra. The SD-PSD threshold for greigite cubes was lower (d 0 ≈ 52 nm) than that of Muxworthy et al. (2013) (= 58 nm); however, we used a different value for M S . When we modelled using the same parameters as Muxworthy et al. (2013) we found excellent agreement, obtaining a d 0 ≈ 60 nm.
NEB method calculations of the room-temperature blocking volumes of greigite show the importance of thermal fluctuations in determining the magnetic structure of a ferromagnetic grain beyond simpler micromagnetic models. We found that equant SD greigite grains cannot be expected to be reliable palaeomagnetic recorders: only PSD grains larger than ∼72 nm are able to hold a remanence in geological timescales, though there is a strong grain shape dependence.
We found that cuboctahedral particles have a blocking volume of ∼134 nm, a volume almost seven times larger than for the (truncated) octahedral particles. This highlights the importance of addressing not only the size distribution in sedimentary rock samples, but also the shapes of the greigite grains present.
We found that the transitions between PSD states occur via well-defined states by a structured rotation of the vortex cores. This could be useful for future analytical models beyond the SD coherent rotation theory. Fig. 9. Relaxation times for the euhedral shapes. The relaxation times are calculated with τ 0 = 1 ns, T = 300 K. τ increases exponentially for the SD-SD transitions, but reaches only laboratory scale stabilities before the SD-PSD threshold. The (truncated) octahedral shapes show a dip in the stability in the lower end of the PSD range for HAV-HAV and IAV-IAV transitions, reaching a minimum at ∼52 nm. Once the EAVs have the lowest energies, the relaxation times grow exponentially with size, reaching four billion years from ∼70 nm. The cuboctahedra, however, show a different behaviour: they do not show a dip in the stability for the first HAV-HAV transitions, but with increasing size they produce lower relaxation times as the energy of the dIAV gets closer. The relaxation times for the transitions between dIAVs plateau up to 74 nm. Exponential growth of the relaxation times occurs from ∼110 nm (inset), surpassing four billion years from ∼134 nm.
Although in this study we have disregarded the effects of particle elongations, which are sure to effect the blocking volumes (Muxworthy et al., 2013), our results are representative of widespread authigenic greigite grains of abiotic origin. The effects of grain-grain magnetostatic interactions have also been excluded from this study. Isolated greigite grains are not very common, rather, natural samples show that greigite is more commonly aggregated in tight clusters. However, this study provides a steppingstone towards understanding the more complex greigite clusters.