Bending and Collapse: Magnetic Recording Fidelity of Magnetofossils From Micromagnetic Simulation

Magnetosome chains produced by magnetotactic bacteria are important paleoenvironmental and paleomagnetic recorders. It has been shown that magnetic properties of magnetosome chains are closely related to their morphology and chain structures; however, the in situ structures of magnetosome chains in sediments (magnetofossils) are not known. Magnetosome chains are subject to various deformations after cell dissolution and are therefore unlikely to be fully intact, obscuring their original magnetic signals. Here, we use finite element micromagnetic simulations to quantify changes in magnetic signals in response to chain deformation, in particular, as a function of variable degrees of bending and collapse. Our results indicate that bending/collapse leads to a significant coercivity reduction and domain state transition of the chain. Therefore, hysteresis parameters can be used to assess the degree of chain bending/collapse in magnetofossil‐rich sediments. Calculations of the contributions of chain bending/collapse to the post‐depositional remanent magnetization (pDRM) of magnetofossils indicate that pDRM remains both faithful to the pre‐bending/collapse natural remanent magnetization, and that the remanence of some structurally deformed magnetofossil assemblages remains thermally stable over billion‐year timescales, suggesting that even strongly deformed magnetosome chains in ancient geological materials retain faithful paleomagnetic records and thus have potentials for tracing ancient geomagnetic field variations and microbial activities on early Earth.

2 of 17 et al., 2015,2018). Moreover, folded chains and flux-closure rings have been observed in attempts to create colloidal magnetosome suspensions (Philipse & Maas, 2002). It is likely that magnetofossils in sediments are not preserved in their original intact structures after a range of complex sedimentary processes that can cause magnetosome chains to bend or collapse; for example, during settlement through the water column in profundal facies (Lascu & Plank, 2013) and post-depositional sedimentary diagenesis (McNeill & Kirschvink, 1993). Furthermore, magnetic mineral dissolution caused by diagenesis may affect magnetofossil preservation (Roberts, 2015). Bacterial culture experiments (Kobayashi et al., 2006;Li et al., 2012) and numerical simulations (Amor et al., 2022;Chang et al., 2019;Harrison & Lascu, 2014) revealed that the degree of bending or collapse significantly affects rock magnetic properties of magnetosome chains. However, the in situ magnetofossil chain structures in sediments have not been reported. A careful rock magnetic investigation of magnetofossil-rich sediments suggested an upper limit of the magnetofossil concentration at ∼0.14% by volume and typical distances between chains of >9 times their length (Ludwig et al., 2013). This means that MTB cells were originally separated by sediments matrices, and that chain clumping after cell dissolution did not occur, unlike results from chain disruption experiments using bacterial cultures (e.g., Kopp et al., 2006;Li et al., 2012). Moreover, magnetosome chains may adhere to the surfaces of clay minerals in sediments (Galindo-Gonzalez et al., 2009), thereby stabilizing the chain structure. Since the in situ chain structures are unknown, we carry out systematical forward modeling to explore the effect of two possible modes of chain structure deformation on their magnetic properties in sediments: a regular bending model from straight chains to rings, and a random-walk collapse model.
Early studies have shown that magnetofossils can acquire natural remanent magnetization (NRM) through two possible mechanisms: (a) If magnetosome chains aligning with the geomagnetic field maintain their orientation after MTB death, magnetofossils can acquire a biogeochemical remanent magnetization (BgRM; Abrajevitch & Kodama, 2009;Tarduno et al., 1998); (b) if magnetosome chains reorient due to the geomagnetic field in the water column or through other forces (e.g., bioturbation and compaction) during sedimentary processes after MTB death, magnetofossils can acquire a depositional (or post-depositional) remanent magnetization (DRM or pDRM; Heslop et al., 2013). Chain bending/collapse may occur during sedimentary compaction, possibly leading to magnetofossils acquiring post-depositional remanent magnetization (pDRM) that could obscure an earlier NRM. Furthermore, chain collapse may also reduce the paleomagnetic stability of any remanent magnetization, thereby losing any NRM altogether. However, these effects have not been quantitatively assessed, and hence, the paleomagnetic reliability and stability of bent or collapsed magnetofossil chains over geological times are currently unclear.
Biogenic magnetite has the potential of being used as a paleoenvironmental indicator due to the responsiveness of MTB and their biosynthesized magnetic particles to environmental conditions (e.g., Chang et al., 2018;Roberts et al., 2011;Yamazaki, 2012;Yamazaki & Kawahata, 1998). Hysteresis properties are fast and easy to measure and potentially allow to infer types and amounts of magnetosomes in sediments, which can be used as a magnetic paleoenvironmental proxy (e.g., Berndt et al., 2020;Chang et al., 2018;Egli, 2004). While hysteresis parameters at least partially depend on the size and shape of biogenic particles, that is, factors controlled by the paleoenvironment, micromagnetic simulations in FORCulator indicate that magnetic hysteresis properties of magnetosomes are sensitive to chain structures Harrison & Lascu, 2014). However, the point dipole approximation by FORCulator neglects-and hence cannot predict-particle dimensions. Recently, a systematic extensive three-dimensional finite-element micromagnetic simulation was performed to establish the relationship between hysteresis properties and both particle and chain morphology (i.e., magnetosome size, elongation, intra-chain spacing, and chain length) of straight magnetosome chains (Berndt et al., 2020). Here, we extend this model to bent and collapsed magnetofossils to address two key questions: (a) How does chain bending/collapse affect magnetic hysteresis parameters (which are routinely used in paleoenvironmental studies)? (b) Is any pDRM carried by bending/collapse chain faithful to the pre-bending/collapse NRM (i.e., in the same direction and with a proportional intensity), and thermally stable over geological timescales?
We generated more than 200 magnetosome models with variable degrees of chain bending and collapse, and calculated hysteresis loops, back-field curves, and first-order reversal curves (FORC). We also simulated pDRM after chain collapse/bending and calculated blocking temperatures to assess the paleomagnetic recording reliability of magnetofossil assemblages in sediments and sedimentary rocks.

Micromagnetic Models
Two types of geometrical models were built: (a) a bending model where the chain forms a regular arc (Figures 1a and 1b), and (b) a collapse model (Figures 1c and 1d), where particles follow a three-dimensional random-walk (Harrison & Lascu, 2014). All models contain 10 cubic or cuboctahedral particles spaced at an interparticle spacing d and forming a bending angle b or collapse angle c. Based on a grain size analysis of magnetosomes from transmission electron microscopic images (Berndt et al., 2020), we chose only the most common crystal sizes for modeling, including equant and elongated (with an aspect ratio of 1.2) 40-nm grains with tight (i.e., d = 3 nm) and loose (i.e., d = 8 nm) interparticle spacing, as well as tightly spaced (i.e., d = 3 nm), equant 100-nm particles (Table S1 in Supporting Information S1). The center angle formed by 10 particles arranged in a complete ring is 324°, so regular bending angles ranging from 0° to 324° with a step size of 30° were modeled. We modeled collapse angles in the range of 0°-70° with a step size of 10° (since higher bending degrees cause overlapping of particles). Ten random-walk collapse models were generated for each collapse angle.
The models were meshed with the finite-element software Trelis 16.3 (Trelis, 2021) with a mesh size of 9 nm. Micromagnetic simulations were carried out in a finite element micromagnetic modeling package, MERRILL version 1.3.5 (Ó Conbhuí et al., 2018). The cubic anisotropy axes of all particles in a chain were rotated with the corresponding bending/collapse angles to align with the chain axis. Simulation results were visualized using ParaView 5.5.2 (Ahrens et al., 2005).

Hysteresis Simulations
As hysteresis loops of magnetosome chains depend strongly on applied field direction (Li et al., 2013), we use 100 random field directions (in 3-D space) for 40-nm bending models, 25 for 100-nm bending models, 10 for 40-nm collapse models, and 5 for 100-nm collapse models due to the limitation of computing time (Table S1 in Supporting Information S1). The applied field ranges from −300 to +300 mT with a step size of 1 mT at low fields (0-5 mT), 5 mT at medium fields (5-100 mT), and 20 mT at high fields (>100 mT). The same random magnetic field directions were used for the back-field curve calculations, as follows: (a) calculating the saturation remanence with a positive 300-mT saturation field (B s ), (b) applying fields from 5 to 300 mT successively in (2) rotating the new particle by the collapse angle c, around the axis perpendicular to the plane formed by the existing and new particle; (3) rotating the new particle by a random angle around the axis parallel to the axis formed by the existing and new particle when it was first created (Harrison & Lascu, 2014). the negative field direction and obtaining remanence at each field step. The step size is −5 mT from −5 mT to −100 mT and −20 mT from −100 mT to −300 mT. The same particle size (i.e., 40 nm) and crystalline forms (i.e., cuboctahedra) and 100 random field directions were used for FORC simulation. The applied field ranges from −200 to +200 mT with a step size of 2 mT. We used 100 random field directions to simulate FORCs for straight and bending models. 20 random-walk collapse models were established for each collapse angle and FORCs in 20 random external magnetic field directions were simulated.

pDRM Simulations
To investigate the paleomagnetic reliability and stability of bent or collapsed magnetofossils over geological times, an initially straight chain with a saturation remanent magnetization state (SIRM; i.e., the pre-bending/ collapse NRM) was progressively bent/collapsed in steps of 30°/10°, respectively. Then, the magnetic moment vectors at all vertices of the mesh were rotated and loaded into the model with the next higher bending angle. Next, the energy of this magnetic state was reminimized in zero field to obtain the magnetic state. The magnetization calculated at each step is used to simulate pDRM with increasing bending or collapse angle. Moreover, the effect of changes in pDRM intensity relative to straight chains due to the geometric rotation of magnetic moments in each particle was calculated analytically. If the pDRM was perfectly faithful to the pre-bending/collapse NRM, the new state should be identical to a purely geometric rotation of all magnetic moments; if on the other hand, a change in the domain state of the chain or individual grain occurs, the initial NRM signal can become corrupted. Similarly, we simulated pDRM in the same random directions as hysteresis simulations.
Energy barriers between two local energy minima (LEM) states can be calculated in MERRILL using the Nudged-elastic-band (NEB) method, which can determine the minimum energy transition between two given LEM states (Nagy et al., 2017;Ó Conbhuí et al., 2018). Considering that the particles in our chain model have strong interaction, we regard the chain as a single system for the NEB calculations. When calculating the minimum energy path, we set two antiparallel stable remanent states along the chains, generated in the above pDRM simulations as two endpoints. Energy barriers between two antiparallel endpoints were calculated to obtain the thermal stability of the pDRM. Moreover, we calculated the energy of 100 random minimized magnetization states to check whether there are any other possible stable states and whether the chain-aligned remanence state has the lowest energy.
The NEB method for calculating energy barriers using the constraint of minimum action (Fabian & Shcherbakov, 2018) requires that there should be no intermediate minima between the LEM endpoints. For the energy paths with obvious intermediate minima, we extracted these intermediate minimum states as the new endpoint to split the path, and then recalculated the minimum energy paths until there are no intermediate minima in the new paths. From these splitting paths, we chose those that provided the lowest energy barrier to calculate the relaxation time τ (Nagy et al., 2017;Néel, 1949). The blocking temperature (T B ) was then calculated by repeating the calculation at multiple temperatures up to just below the Curie temperature and selecting the temperature where the relaxation time is of the order of minutes (i.e., τ = 100 s; Nagy et al., 2017).
We calculated minimum energy paths of some typical magnetosome chain models at different temperatures using the NEB method ( Figure S1 in Supporting Information S1). These minimum energy paths contain fewer intermediate minima at high temperatures (500°C) compared to lower temperatures (<400°C; Figure S1 in Supporting Information S1), reflecting the increased path stability as the magnetostatic interactions, which is proportional to the saturation magnetization decrease with temperature. Minimum energy paths of the straight chain model (b = 0°) do not contain intermediate minima at 500°C (Figure S1d in Supporting Information S1). Minimum energy paths of a ring-shaped model (b = 324°) and a medium collapse model (c = 30°) have only one obvious intermediate minimum at 500°C (Figures S1l and S1p in Supporting Information S1). Minimum energy paths of a strong collapse model (c = 70°) have two obvious minima at 400° and 500°C (Figures S1s and S1t in Supporting Information S1). Only the minimum energy paths of a bending model (b = 180°) still contain much noise at high temperatures (Figures S1g and S1h in Supporting Information S1). Hence, we chose straight chain model, ring-shaped model, 30° collapse model, and 70° collapse model to calculate energy barriers at high temperature ranges using the path splitting method. Then, we used polynomial fitting to estimate their blocking temperature (Nagy et al., 2017). 10.1029/2021JB023447 5 of 17

Hysteresis Loops and Back-Field Curves
Our simulations of hysteresis loops ( Figure S2 in Supporting Information S1) and back-field curves ( Figure  S3 in Supporting Information S1) show that both the coercivity (B c ) and the coercivity of remanence (B cr ) decrease gradually by about 50% with progressive bending/collapse angle for completely bent/collapsed chains (Figures 2a-2d). The ratio of saturation remanence to saturation magnetization (M rs /M s ) increases slightly from 0.5 (the theoretical value for uniaxial particles) to ∼0.6 for 40-nm bending models, indicating a deviation from perfect uniaxiality, but decreases significantly for the 40-nm collapse model (Figures 2e and 2f). However, both bent and collapsed 100-nm models have strongly decreasing values (Figures 2e and 2f), because the domain states of the individual particles switch from the single domain (SD) to single vortex (SV) with increasing bending/ collapse angle.
Hysteresis loops of bent and collapsed magnetofossils differ from those of straight chains, which only have one major switch (Berndt et al., 2020): Chains with large bending angles can switch twice (Figure 3), while strongly collapsed chains can switch multiple times (Figure 4). For strongly bent chains, which are almost circular, the first switching is characterized by domain alignment approximately along the chain direction (i.e., from ① to ② in Figure 3); the second switching corresponds to domain alignment along the external field direction (i.e., from ③ to ④ in Figure 3). These two switches manifest as two distinct jumps in a hysteresis loop (Figure 3). Collapsed models may have multiple switches, leading to segmented jumps in the hysteresis loop (e.g., Figure 4).
Saturation remanent states of grains in all 40-nm bending and collapse magnetosome models exhibit SD behavior but show multiple orientations of magnetic moments: (a) Similar to straight chains (Berndt et al., 2020), moment directions in slightly bent models are parallel to the chain axis in all directions of external fields (Figures 5a and 5c); (b) some moments in strongly bent chains are perpendicular to the chain axis when the external field direction is perpendicular to the chain ( Figure 5b); (c) collapsed chains have multiple magnetic moment directions depending on chain structures (Figure 5e). Particles at both ends of 100-nm chains with slight bending are in SV states (Figure 5c), but all particles in a chain are in SV states for 100-nm models with strong bending (Figures 5d and 5f).

FORC Diagrams
Simulations reveal a significant effect of bending and collapse on FORC diagrams. Straight chains show the typical central ridge in FORC diagrams, along with a negative distribution in the lower-left region due to reversible magnetization of randomly oriented chains (Figure 6a; Chang et al., 2019;Harrison & Lascu, 2014;Newell, 2005). FORC diagrams of ring-shaped models (b = 300°) contain three distinct peaks (Figure 6b), which is due to transient behavior of the nucleation and annihilation of the super vortex that spans the entire ring (Carvallo et al., 2003;Zhao et al., 2017). Similarly, collapse leads to a shift of the central ridge toward lower coercivities along with a vertical broadening of the signal (c = 30°, Figure 6c) and eventually to the formation of "wings" similar to those caused by an increase of the 3-D interaction fields (c = 60°, Figure 6d).

Contributions of Bending/Collapsed Chains to pDRM and Their Thermal Stability
When magnetofossil chains are bent/collapsed during deposition, the magnetic remanence is reduced due to (a) the geometric rotation of particles away from the original deposition axis, and (b) changes in domain states of Figure 3. Simulated hysteresis loop and the switching mode of domain states for a strongly bending magnetosome chain (40 nm, cuboctahedra, loose spacing, and b = 300°). Blue arrows highlight the domain direction. Black arrow represents the direction of the external magnetic field. Small arrows within particles show magnetic moments and were colored by helicity. either the individual particle or the whole chain. Our models show that the reduction of pDRM intensity in a single chain due to bending/collapse can almost entirely be explained by the geometrical rotation of particles, which leads to 100% and 60% reduction for ring-shaped and fully collapsed chains respectively (Figures 7a,  7b, and S4a in Supporting Information S1). For the 40-nm models, their magnetic states are always SD and are consistent with the geometric bending/collapse direction (Figures 8a, 8c, and 8e). Only the 100-nm models deviate slightly from this, which is because the domain states change from SD to SV for the 100-nm particles with increasing bending or collapse (Figures 8b, 8d, and 8f). Simulated change in the normalized pDRM intensity (pDRM/SIRM) for deformed chains relative to straight chains gradually decreases to zero with increasing bending degree (Figures 7c and S4c in Supporting Information S1) but first increase slightly and then decrease continuously for collapse angle > ∼20° (Figure 7d), which is caused by the relatively faster SIRM decrease with collapse than the pDRM for the collapse model (Figures 2f and 7b).
Simulation results show that the magnetization states with magnetic moment aligned along the chain (i.e., endpoint states in our study) are the lowest energy states ( Figure S5 in Supporting Information S1). There are some other possible stable magnetization states, but their energy is higher than that of chain-aligned states (Figure  (Figures 9b and 10). The T B of the ring-shaped model, 30° collapse model, and 70° collapse model are 551°, 519°, and 354°C, respectively ( Figure 11). Therefore, chain collapse does lead to a reduction in T B , but it still remains high enough to be stable over geological time scales. Hence, these four models with different chain structures have high T B , indicating the high thermal stability of magnetosome chain remanence.

Assessing Collapse Degree From Rock Magnetic Properties
Our simulations indicate that bending and collapse significantly affects the magnetic properties of magnetosome chains, especially the reduction of coercivity (Figures 2a-2d), which is consistent with previous micromagnetic simulations Harrison & Lascu, 2014) and experimental data (Kobayashi et al., 2006;Li et al., 2012). A change in the structure of the interaction fields, which become less aligned in the chain direction, leads to the reduction of coercivity through a number of effects. First, the external magnetic field is at an angle to the magnetization so that the magnetostatic interaction energy is reduced. Furthermore, particle gaps in a chain increase with increasing bending and collapse, and therefore significantly reduce coercivity (Berndt et al., 2020). Moreover, a gradual increase in bending or collapse degree disrupts the ordered chain arrangement in the chain direction, which makes it easier for individual grains to switch their magnetization, resulting in multiple switching in some chains (Figures 3 and 4), hence, reducing the coercivity.
Currently, it is difficult to directly quantify magnetosome chain structures within sediment samples. Extracted magnetofossil assemblages imaged by electron microscope observations do not represent their in situ microstructures in sediments due to deformation of chain structures during sample preparation. Some studies proposed to evaluate the bending or collapse degree of magnetosome chains using comprehensive rock magnetic parameters. It has been shown that collapse of magnetosome chains lowers the ratio of anhysteretic remanent magnetization susceptibility to saturation isothermal remanent magnetization (χ ARM /SIRM; Lascu & Plank, 2013;Li et al., 2012), and produces a lower delta ratio between losses of field cooled (FC) and zero-field cooled remanent magnetization across the Verwey transition in low-temperature experiments (δ FC /δ ZFC ; Carter-Stiglitz et al., 2002;Li et al., 2012).
Here, we propose a simple strategy to estimate chain collapse from hysteresis experiments of magnetofossil-dominated samples. Our simulations show that data of collapse models move significantly to the lower-left region with progressive collapse angle in the Néel diagram (Néel, 1955, Figure 12). Previous experimental data show a consistent trend ( Figure 12): Intact magnetosome chains in MTB have a ratio of ∼0.5 and higher coercivity (25-40 mT;Jovane et al., 2012;Li et al., 2012;Lin & Pan, 2009;Moskowitz et al., 1988;Pan et al., 2005), but collapsed and clumped chains from culture experiments have a lower M rs /M s ratio (the lowest ratio is close to 0.2) and coercivity (even less than 10 mT; Li et al., 2012). Some magnetofossil-rich sediments have a similar M rs /M s ratio and coercivity to the medium and strong collapse model (Figure 12; Chang et al., 2018;Ludwig et al., 2013), indicating that magnetosome chains in natural samples may have medium to high collapse angles. In addition, data of bending models move slightly to the upper-left region with progressive bending angle in the Néel diagram ( Figure S6 in Supporting Information S1). In particular, the M rs /M s ratio of ring-shaped chain reaches ∼0.6, which is different from data trend of bacterial culture experiments and natural sediment samples. Hence, we speculate that while slight bending occurs in magnetosome chains in intact MTB cells, magnetofossils in compacted sediments are unlikely to experience such a systematic bending. Rather, the chain collapse model is more likely to be representative of deformed magnetofossils in sediments than the regular bending model. Furthermore, our simulations reveal that FORC distributions are also sensitive to chain bending and collapse ( Figure 6). However, these results are obtained by forward simulation of pure magnetosome chains. It is necessary to consider the influence of magnetic minerals of other origins in sediments on the bulk magnetic properties. Therefore, we can first separate the magnetofossil component using magnetic unmixing techniques such as isothermal remanent magnetism curves (Egli, 2004), central ridge analysis of FORC diagrams (Egli et al., 2010)   diagram and FORC diagrams to qualitatively evaluate chain collapse in natural magnetofossil-dominated samples based on our forward simulation.

Paleomagnetic Recording Fidelity of Collapsed Magnetofossils
Magnetofossils can effectively acquire NRM in suitable depositional environments (Mao et al., 2014;Paterson et al., 2013) and can dominate the NRM in some sediments . However, magnetofossils may acquire pDRM through chain bending/collapse during sedimentary processes. Our simulation results indicate that the pDRM intensity contributed by magnetofossils is sensitive to the degree of chain bending and collapse (Figures 7a-7d). Importantly, our micromagnetic and analytical calculations indicate that the geometry of particle rotation-rather than changes in domain state-is the major cause of pDRM reduction in collapsed magnetosome chains (Figures 7a and 7b), indicating that pDRMs are faithful to the pre-bending/ collapse NRM, albeit with a weaker intensity. Therefore, we can first obtain the magnetofossil concentrations in sediments using rock magnetic unmixing techniques, and then changes in relative paleointensity due to collapse can be corrected if the chain collapse degree in magnetofossil-dominated sediments can be evaluated.
Biomineralization of magnetosomes in bacteria is thought to have existed for a long time in geological history (Lin et al., 2017(Lin et al., , 2018(Lin et al., , 2020 with the earliest geological evidence of magnetofossils having an age of ∼2 Ga (Chang & Kirschvink, 1989). Despite the billion-year timescales of remanences carried by collapsed magnetofossil chains, our simulations suggest that even some strongly bent or collapsed magnetosomes have high blocking temperatures and are capable of retaining stable remanence over geological timescales. Hence, although the pDRM is weakened by post-depositional processes, magnetofossils in ancient rocks can still effectively preserve the paleomagnetic fields of the early Earth.

Conclusions
Our quantitative micromagnetic simulations indicate that bending and collapsing of magnetosome chains have a strong effect on rock magnetic properties and paleomagnetic records of magnetofossil-bearing samples, with apparent coercivity reductions, transition in the domain state of the whole chain, and pDRM intensity reductions with increasing bending or collapse. More importantly, our simulations of blocking temperature indicate that bent or collapsed magnetosome chains can retain a stable pDRM over billion-year timescales. We also show that pDRMs are not only stable but also faithful to the pre-bending/collapse NRM, both in direction and in intensity through a simple geometric reduction in magnitude. The collapse degree of chains can be estimated qualitatively from Néel diagram and FORC diagrams, which may allow correction of the reduction in paleointensity of magnetofossil-dominated samples due to chain collapse. Our simulation provides new quantitative constraints on understanding the magnetic signatures of magnetofossils, which are important for reconstructing paleoenvironmental conditions and geomagnetic field variations.

Data Availability Statement
The complete source code used to produce the data is available on https://doi.org/10.5281/zenodo.5881792. The calculated hysteresis data and energy barriers data are available on the Peking University Open Research Data Platform (https://doi.org/10.18170/DVN/P7KDKR).  Néel, 1955) of modeled hysteresis parameters in this study and experimental data of magnetotactic bacteria (Jovane et al., 2012;Li et al., 2012;Lin & Pan, 2009;Moskowitz et al., 1988;Pan et al., 2005) and magnetofossil-rich sediments (Chang et al., 2018;Ludwig et al., 2013). Gray arrow indicates the direction of the increasing collapse angle. Numbers along the data lines represent collapse angle.