Magnetic field-temperature phase diagram of multiferroic (NH4)2FeCl5·H2O

Owing to their overall low energy scales, flexible molecular architectures, and ease of chemical substitution, molecule-based multiferroics are extraordinarily responsive to external stimuli and exhibit remarkably rich phase diagrams. Even so, the stability and microscopic properties of various magnetic states in close proximity to quantum critical points are highly under-explored in these materials. Inspired by these opportunities, we combined pulsed-field magnetization, first-principles calculations, and numerical simulations to reveal the magnetic field–temperature (B–T) phase diagram of multiferroic (NH4)2FeCl5⋅H2O. In this system, a network of intermolecular hydrogen and halogen bonds creates a competing set of exchange interactions that generates additional structure in the phase diagram—both in the vicinity of the spin flop and near the 30 T transition to the fully saturated state. Consequently, the phase diagrams of (NH4)2FeCl5⋅H2O and its deuterated analog are much more complex than those of other molecule-based multiferroics. The entire series of coupled electric and magnetic transitions can be accessed with a powered magnet, opening the door to exploration and control of properties in this and related materials.


INTRODUCTION
With their strong coupling between charge, structure, and magnetism, multiferroics are ideal platforms for exploring quantum phase transitions. 1,2 These transitions are controlled by external stimuli such as voltage, magnetic field, pressure, and light -quite different from traditional phase transitions that rely upon thermal fluctuations. [3][4][5][6][7][8][9] Molecule-based multiferroics offer several advantages over their oxide counterparts for the study of quantum phase transitions. These include overall low-energy scales, flexible molecular architectures, and ease of chemical substitution. [10][11][12][13] Thus, in addition to topological similarities with oxide multiferroics, the molecular analogs tend to have lower critical fields that are more easily reached in the laboratory. 14-17 A significant bottleneck to the investigation of molecule-based multiferroics and their application for ultra low-power multi-state memory, switching devices, and novel computing architectures [18][19][20][21] is the inadequate understanding of the different equilibrium states and their proximity to quantum critical points. Uncovering the full magnetic field-temperature (B-T) phase diagram is a significant step toward addressing these issues. 10,22 At the same time, revealing the nature of the various magnetic quantum phase transitions and the character of the competing phases provides important opportunities for property control and tests our theoretical understanding of quantum magnets. 10,[15][16][17]23 These ideas can be explored in the erythrosiderites, the class of compounds with chemical formula A 2 FeX 5 ⋅H 2 O, where A is ammonium or an alkali metal and X is a halide. We focus on multiferroic (NH 4 ) 2 FeCl 5 ⋅H 2 O and its deuterated analog due to their spontaneous electric polarization and unusual magnetic behavior. 24,25 The orthorhombic unit cell consists of four octahedral [FeCl 5 ⋅H 2 O] 2− groups and eight NH þ 4 ions that are all symmetrically equivalent. 24 Fig. 1c summarizes the temperaturedriven transitions. These include (i) an order-disorder transition involving NH þ 4 [T O/D ≈ 79 K], (ii) the Néel transition below which incommensurate-collinear sinusoidal magnetic order appears [T N ≈ 7.25 K], and (iii) a transition to an incommensurate-cycloidal spin state that supports a spontaneous electric polarization primarily along the a-axis [T FE ≈ 6.9 K]. [26][27][28][29] A fascinating series of spin and polarization flops occur with increasing magnetic field. 24,25 These magnetically-driven transitions will be discussed in greater detail below. There is an abundance of hydrogen and halogen bonding between neighboring [FeCl 5 ⋅H 2 O] 2− octahedra. As shown in Fig. 1a, b, these intermolecular interactions form a set of superexchange pathways linking the Fe 3+ S ¼ 5 2 À Á centers. [28][29][30][31] The full set of magnetic interaction pathways includes Fe-O-H⋯Cl-Fe, Fe-Cl⋯Cl-Fe, and Fe-Cl⋯O-Fe linkages. The magnetic interaction through O-H⋯Cl is the strongest and forms a zig-zag chain along the b-axis. 24 As (NH 4 ) 2 FeCl 5 ⋅H 2 O is a type-II multiferroic, polarization derives from the spin configuration. Therefore, one cannot begin to understand polarization, polarization flops, and magnetoelectric coupling without first examining the underlying spin system and determining the B-T phase diagram.
We combined pulsed-field magnetization with complementary first-principles calculations and numerical simulations to reveal the magnetic properties and complete B-T phase diagram of multiferroic (NH 4 ) 2 FeCl 5 ⋅H 2 O. Significantly, all of the magnetically driven transitions-including that to the fully saturated state-are within the range of experimentally available powered magnets. The overall complexity of the phase diagram arises from the competing exchange pathways (J 1 -J 5 ) produced by an elaborate network of intermolecular hydrogen and halogen bonds. The magnetic quantum phase transition and its satellites are the most striking of the sequence, and the discovery of these structures is one of the major factors that make this system a quantum material. That they occur in a compound where the exchange interactions consist entirely of intermolecular hydrogen and halogen bonds but still sport similarities with other well-known multiferroics such as TbMnO 3 and MnWO 4 is remarkable. [32][33][34][35][36][37] Structure-property relations involving counterion substitution and spin-state noncollinearity are also discussed.

RESULTS AND DISCUSSION
Magnetic properties and field-induced transitions in (NH 4 ) 2 FeCl 5 ⋅H 2 O Figure 2 summarizes the pulsed-field magnetization of (NH 4 ) 2 FeCl 5 ⋅H 2 O and its deuterated analog. There are numerous field-induced transitions depending on the isotropic decoration and the direction of the applied field. They naturally separate into two groups: (i) a series of low field spin reorientations below 6 T and (ii) the saturation field and its satellites near 30 T. These satellites are evident in Fig. 2e, f and are associated with competing exchange pathways. They are discussed in detail below. The low field features are small but clearly revealed for B || c, whereas they are more complicated for B || a. The transition to the fully polarized state (B Sat ) is very distinct and occurs near 30 T in each case. This is the energy scale above which all frustration is relieved. Depending on the field direction, we find one or two small satellites in the vicinity of B Sat -evidence for a series of quasi-isoenergetic phases near the quantum phase transition. These values are summarized in Table 1.
The overall shape of the magnetization is consistent with expectations for a three-dimensional (reasonably isotropic with S > 1 2 ) material, with a linear rather than concave rise on approach to B Sat . 17 This is in line with our numerical modeling (discussed below) from which we extract the various exchange interactions.
We can also estimate the value of the primary exchange interaction from the size of the ultimate critical field. Taking B Sat = 30.3 T for the field along the c direction in the hydrogenated sample and using a simple Hamiltonian 38 with one exchange constant H ¼ ÀðJ=2Þ This value is in reasonable agreement with J 1 from our modeling calculations. Here, the primary exchange interaction corresponds to the Fe-O-H⋯Cl-Fe pathway as shown in Fig. 1b. It also compares well with that of the isomorphic K 2 FeCl 5 ⋅H 2 O analog (J 1 = −1.65 K). 31 The saturation fields in Table 1 reveal several interesting structure-property relations. In each case, B Sat ≈ 30 T, with modest variations that depend upon isotopic substitution and field direction. This suggests that the anisotropies and Dyaloshinski-Moriya interactions are not large in (NH 4 ) 2 FeCl 5 ⋅H 2 O. This is consistent with our predictions for small anisotropies, summarized below. Additional evidence for overall threedimensional antiferromagnetism comes from the linear magnetization between 5 and 30 T. We also find that B Sat,D < B Sat,H . This is a very common trend-similar isotope effects were reported in [Ni (HF 2 )(pyz-d 4 ) 2 ]SbF 6 and CuF 2 (pyz)(H 2 O) 2 . 39,40 Deuteration typically reduces the exchange interactions by a few percent, because a heavier atom has a smaller excursion from its equilibrium position in the anharmonic potential. 39,40 Deuterium substitution also tends to smear out the multiple lower field transitions. Moreover, we find that B Sat,a > B Sat,c . This deviation from isotopic behavior is probably related to g-factor anisotropy. 17 There is also a more complex set of magnetically driven transitions when B || a.
Uncovering the spin structure and magnetic exchange interactions To further explore these ideas, we carried out several different numerical simulations, first to understand the zero-field ground state and then to interpret the field-induced transitions. In zero field, the Hamiltonian can be written as . The magnetic coupling is three dimensional: strong quasi-two-dimensional interactions with antiferromagnetically coupled planes. b Calculated J values against distance between Fe 3+ ions. c Schematic depicting the sequence of temperature and magnetic field transitions that occur in (NH 4 ) 2 FeCl 5 ⋅H 2 O. d Ratio of J interactions with respect to J 1 depicting how the collinear and non-collinear states are stabilized with respect to the nature of the counterion. The blue star in the non-collinear state corresponds to (NH 4 ) 2 FeCl 5 ⋅H 2 O. The red star in the collinear state corresponds to K 2 FeCl 5 ⋅H 2 O whose exchange constants are taken from ref. 31 where J ij is the isotropic exchange interaction between S i and S j , K b is the single-ion anisotropy, and S i are classical spin vectors. The actual value of the parameters used in this simulation were obtained from density functional theory (DFT). Details are available in the Methods section and Supplementary Information.
Our numerical simulations produce spin structures that are remarkably consistent with the experimental results (Fig. 3a). Moreover, the experimentally observed wave vector of the spin spiral is Q ≈ 0.228 ± 0.02 reciprocal lattice units 26 , which compares well with the numerically calculated value of Q ≈ 0.205 ± 0.001 (the small difference in Q may be attributed to the effects of single-ion anisotropy along the a-and c-axis). These results demonstrate the validity of the exchange interaction parameters obtained from DFT and our model Hamiltonian. To determine the  origin of non-collinear spin structure, we investigated how exchange interactions change the magnetic ground state. Interestingly, even in the absence of J 1 , J 3 , and J 5 , we can reproduce the non-collinear spin structure ( Supplementary Fig.  S3). However, in the absence of J 2 or J 4 , which form a triangular plaquette (Fig. 1a), a collinear spin structure is stable. Therefore, geometrical frustration with competing antiferromagnetic J 2 and J 4 in a triangular lattice stabilizes the non-collinear spin state propagating along the c direction. The collinear spin structures observed in certain other A 2 FeX 5 ⋅D 2 O erythrosiderites can also be explained within this framework, as J 4 is mediated by the ammonium ions. The latter have very different characteristics compared with alkaline ions. The absence of ammonium ions relieves magnetic frustration by making either J 2 and/or J 4 dominant. Moreover, the interactions between triangular lattice planes are mediated by the strongest antiferromagnetic exchange J 1 . This makes the spin spirals on the different triangular lattice planes almost antiferromagnetic, a prediction that has been experimentally confirmed. 26 In triangular lattice antiferromagnets, exotic magnetic ground states and successive magnetic phase transitions have been observed. 41 Thus, the series of magnetic phase transitions observed in field also arises from the subtly competing exchange interactions in this system. These calculations also suggest that pressure may provide a pathway between the different A 2 FeX 5 ⋅H 2 O materials (Fig. 1d). As the phase diagram indicates, the spin-state transition occurs with increasing J 3 /J 1 . We therefore speculate that uniaxial pressure applied along J 3 could induce a transition from a noncollinear to collinear spin structure. On the other hand, hydrostatic pressure is unlikely to induce a magnetoelastic phase transition, as the J-ratios would be more or less the same.
We also carried out first-principles calculations to obtain the spin densities in the non-collinear antiferromagnetic and ferromagnetic states of (NH 4 ) 2 FeCl 5 ⋅H 2 O. These results are summarized in Fig. 3. In the zero-field (non-collinear antiferromagnetic) state, the spin density resides primarily on the Fe 3+ centers and to a lesser extent on the chlorine and water ligands. The phase of the spin density alternates characteristically 42 and the empty space between ligands is the node between the up-and down-spin states. In the full-field (ferromagnetic) state, the spin-density pattern of Fe 3+ is similar to that of the ligands but with an overall in-phase (rather than out-of-phase) arrangement. Another difference in the spin density can be found on the water ligand extending toward the nearest chlorine through the H-O⋯Cl hydrogen bonding pathway. As discussed previously, the H-O⋯Cl pathway produces the primary exchange interaction J 1 . The additional spin density across this intermolecular hydrogen bond suggests that J 1 increases in the fully saturated state. Hence, the field-induced change in spin density is linked to modifications of the magnetic exchange interactions. This spin density redistribution also contributes to the lack of inversion symmetry by exacerbating the anti-alignment of polarization and magnetization.
Revealing the phase diagram of (NH 4 ) 2 FeCl 5 ⋅H 2 O We can use the pulsed-field magnetization data in Fig. 2 to develop the B-T phase diagram. In order to precisely determine the location of the various phase boundaries and track their dependence on the external stimuli, we calculated (∂M/∂B) T and plotted these curves vs. field (Fig. 2e, f). The position of the spin flop and transition to the fully polarized state is clear. The use of derivative techniques is particularly advantageous when following trends in small features. For instance, careful analysis of the derivative structure in the vicinity of the spin-flop transitions yields three phase boundaries (Fig. 2e). Analysis in the vicinity of the saturation field (indicated with an arrow) reveals a series of transitions as well (Fig. 2f). (∂M/∂B) T also indicates the relative importance of magnetism at each phase boundary. These derivative techniques uncover two important regions: (i) the low-field spin reorientations below 6 T and (ii) the high-field region approaching (and including) the fully polarized state at B Sat ≈ 30 T. With increasing temperature, the position and relative importance of the different peaks in the derivative response evolves. Figure 4 displays the B-T phase diagram of (NH 4 ) 2 FeCl 5 ⋅H 2 O for B || c. We label the different regimes in accord with prior magnetization, polarization, magnetostriction, and neutron scattering. [24][25][26][27]29 Here, AFM = antiferromagnetic, PM = paramagnetic, ICC = incommensurate cycloidal, CS = cycloidal spiral, C1 = distorted cycloid, C2 = quasi-collinear, CLS = collinear sinusodial, FE = ferroelectric, and NE = non-electric. Of course, the transition to the fully polarized state (B Sat ≈ 30 T) along with two additional reorientation transitions immediately preceding it are entirely new -providing the first glimpse of the full complexity of the B-T diagram in (NH 4 ) 2 FeCl 5 ⋅H 2 O. Similar phase diagrams for B || a and for the deuterated crystals are available in the Supplementary Information.
With increasing magnetic field, (NH 4 ) 2 FeCl 5 H 2 O undergoes a sequence of transitions where the quantum fluctuations present in the antiferromagnetic ground state are suppressed as the system is driven into the fully polarized state. The primary highfield transition at 30 T is a quantum phase transition driven by external stimuli rather than thermal fluctuations. [8][9][10]22 The magnitude of B Sat is linked to the largest exchange interaction (J 1 ), with a pathway of Fe-O-H ⋯ Cl-Fe. Moreover, the spin density redistributes to favor this pathway in the ferromagnetic state (Fig.  3d). The other exchange interactions are similar in magnitude but smaller. Thus, although we cannot explicitly link J 2-5 to specific phase boundaries, it stands to reason that the more subtle transitions at 26.9 and 28.6 T can be attributed to the lesser, competing J's. For instance, J 2 and J 4 are the next largest superexchange interactions and may be correlated with the weak 28.6 T phase boundary, whereas J 3 and J 5 may be linked with the 26.9 T boundary. Obviously, all of the expected five spin reorientations should be resolved at the lowest temperatures.
In the lower field regime, (NH 4 ) 2 FeCl 5 ⋅H 2 O undergoes a series of transitions in the vicinity of the 3.7 T spin-flop transition (Fig. 4b).
As before, this cascade of transitions is due to the many competing exchange pathways in this compound. This grouping can be linked to the five exchange pathways, although resolution is lost due to similarities in magnitude and thermal broadening. Triangular frustration created by J 2 and J 4 (Fig. 1a) is alleviated as magnetic field increases. With increasing temperature, the lowfield phase boundaries converge toward a pair of triple points between AFM 1, AFM 2, and AFM 3, as well as AFM 1, AFM 3, and CLS/NE phases. Overall, the low-field portion of the magnetic phase diagram is in excellent agreement with that of Ackermann et al. 24 The main exception is the splitting of the 3.7 T spin flop at base temperature (Fig. 2e).
Although the properties and complex multiferroic phases of (NH 4 ) 2 FeCl 5 ⋅H 2 O have been extensively studied in the low-field regime, [24][25][26][27]29 it is interesting to anticipate the high-field response. In general, for type-II multiferroics, the spins need to break spatial inversion symmetry in order to create an electric polarization. Common spin structures that break spatial inversion symmetry are various kinds of spirals. [43][44][45] If the spins are fully aligned with the magnetic field, then there is no spatial inversion symmetry-breaking spin structure available to create electric polarization. In other words, above the 30 T transition to the fully saturated spin state, the polarization is likely to be quenched. The high-field state is therefore unlikely to be magnetoelectric. Things are different at 5 T. Here, the spin configuration is quasi-collinear AFM 3 and the system is in the FE III phase with polarization along c. 24,27 Whether the spin structure, polarization, and magnetoelectric coupling at 5 T [24][25][26][27]29 are similar to those further away from equilibrium-for instance, at 27 T-awaits further study. That ferroelectricity depends so intimately upon the spin configuration is the overarching motivation for revealing the magnetic phases of (NH 4 ) 2 FeCl 5 ⋅H 2 O and unveiling the overall structure of the phase diagram.
What differentiates the B-T phase diagram of (NH 4 ) 2 FeCl 5 ⋅H 2 O from that of other molecule-based multiferroics is the exceptional level of detail. For instance, in metal-organic framework compounds such as (CH 3 ) 2 NH 2 ]Mn(HCOO) 3 , the phase diagram exhibits a spin flop near 0.3 T, a broad canted phase, and a transition to the fully saturated state at 15.3 T. 22 Likewise, in CaCo 2 As 2 , the critical fields are 3.7 and 7.5 T. 46 The overall simplicity of these B-T phase diagrams arises from the superexchange pathways through formate ligands or Co⋯As interactions, respectively. Magnetic exchange, of course, can also take place through hydrogen and halide bonds. This is the situation in (NH 4 ) 2 FeCl 5 ⋅H 2 O and a number of other molecule-based materials. [47][48][49][50][51] What makes (NH 4 ) 2 FeCl 5 ⋅H 2 O so unusual is that there are five isoenergetic intermolecular hydrogen and halide bonds that function as superexchange pathways. In addition to the extraordinary softness, this gives rise to an unusual degree of frustration that is manifested as a series field-driven transitions and an extremely complex B-T phase diagram. Traditionally, these systems have been under-explored, but our work on (NH 4 ) 2 FeCl 5 ⋅H 2 O is a major step toward changing this situation. Other members of the erythrosiderite family such as K 2 FeCl 5 ⋅H 2 O do not seem to be as soft as the ammonium compound. The K +containing system is also a collinear antiferromagnet rather than a multiferroic. This difference has been a puzzle for some time and appears to be related to the character of the superexchange network as discussed in the Supplementary Information. Unexpectedly, the spin structure of TbMnO 3 is similar to that in (NH 4 ) 2 FeCl 5 ⋅H 2 O, 24,29,[32][33][34] providing yet another reason to more deeply examine the properties in this unusual type-II moleculebased multiferroic.
We combined pulsed-field magnetization measurements with first-principles calculations and numerical simulations of exchange interactions to reveal the magnetic properties of the type-II multiferroic (NH 4 ) 2 FeCl 5 ⋅H 2 O. The B-T phase diagram is surprisingly complex with a spin-flop transition between 3 and 5 T, and a transition to the fully saturated state near 30 T-each of which is preceded by a series of weak reorientation transitions that reflect the many quasi-isoenergetic exchange interactions in this material. The latter arise from an elaborate intermolecular hydrogen and halogen bonding network, and in addition to a full description of the non-collinear ground state, we evaluate field-induced changes to the spin density pattern in terms of these intermolecular exchange interactions across the entire field regime. Structure-property relations in (NH 4 ) 2 FeCl 5 ⋅H 2 O are discussed with an emphasis on how isotopic decoration impacts the intermolecular hydrogen and halogen bond network and saturation fields.

Crystal growth
Hydrogentated and deuterated (NH 4 ) 2 FeCl 5 ⋅H 2 O single crystals were grown by solution method using HCl/DCl, FeCl 3 , and NH 4 Cl/ND 4 Cl, respectively. 24,29 A sealed saturated solution was kept in a sample environment chamber at 38°C and allowed to slowly evaporate. Large crystals were obtained. The crystals were characterized by magnetic susceptibility and specific heat measurements, and no significant deuteration-induced effects were recorded between the hydrogenated and deuterated samples. Samples were oriented using morphological faces as references. 24 Magnetization measurements High-field magnetization measurements were performed using a 65 T short-pulse magnet at the National High Magnetic Field Laboratory in Los Alamos using a 1.5 mm bore, 1.5 mm long, 1500-turn compensated-coil susceptometer, constructed from 50-gauge high-purity copper wire as described in ref. 22 . In addition to full field pulses, we carried out 10 and 35 T shots to resolve the spin flop and saturation. The critical fields (B SF , B Sat , and all satellite transitions) were determined using first derivative techniques. The spin transitions were brought together to create the B-T phase diagram.

First-principles density functional theory
First-principles calculations were performed using DFT within the local density approximation LDA + U method as implemented in the Vienna ab initio simulation package (VASP 5.4.1). [52][53][54] We use the Dudarev 55 implementation with on-site Coulomb interaction U = 2.0 eV to treat the localized 3d electron states in Fe. Our value for U is chosen to bring Q closest to the experimental value. The projector augmented wave (PAW) potentials 56,57 explicitly include eight valence electrons for Fe (3d 7 4s 1 ), five for N (2s 3 2p 3 ), one for H(1s 1 ), seven for Cl (3s 2 3p 5 ), and six for O (2s 2 2p 4 ). Before calculating the exchange interactions and spin density, the crystal structure was taken from ref. 58 and the Crystallography open database (COD ID: 9012597). The atomic structure was optimized for both lattice parameters and atomic positions with collinear-antiferromagnetic spin ordering. A plane-wave basis set with a cutoff energy of 500 eV was used. The k-point sampling used the Monkhorst-Pack scheme 59 and employed 2 × 3 × 4 and 2 × 3 × 2 meshes for the unit cell and the supercell of (NH 4 ) 2 FeCl 5 ⋅H 2 O, respectively. The atomic positions were optimized until the interatomic forces were smaller than 1 meV/Å.
To calculate exchange interactions, we used an energy-mapping analysis for localized spins without spin-orbit coupling. 60,61 A 1 × 1 × 2 supercell was chosen to evaluate the exchange interaction J 4 , which cannot be included in a smaller unit cell. A detailed explanation for the computation of J ij is provided in section 2 of the Supplementary Information. In order to simulate the spin state with Q ≈ 0.23 along the c direction, 26,27 a very large supercell would be required. To make this problem more tractable, we assumed that Q = 0.25 and used a 1 × 1 × 4 supercell to reproduce the non-collinear spin ordering. On the other hand, the ferromagnetic spin density was obtained from a collinear spin-polarized scf calculation. The kpoint sampling uses the Monkhorst-Pack scheme 59 and employs a 2 × 3 × 1 and a 2 × 3 × 4 mesh for non-collinear and collinear spin density calculations, respectively.

Numerical simulations
To find the spin ground state at zero magnetic field from Eq. (1), the exchange interactions were taken from the DFT calculations. In zero magnetic field and with no anisotropy in the cycloidal plane, the spin cycloid can be approximated by a single sinusoidal function. 62,63 As singleion anisotropy along the b-axis preserves the harmonic spin cycloid in the ac plane, no S b i R ð Þ component was considered. We can safely assume that the spins are classical vectors, as the Fe 3+ ions have a large spin value 5/2. Thus, the trial functions for S i (R) are with S y i ðRÞ ¼ 0. The index i labels the four distinct spins in the magnetic unit cell (see ref. 26 for spin number assignment in the magnetic unit cell) and ϕ i are the phases. For our model, ϕ 1 = ϕ 2 and ϕ 3 = ϕ 4 .
To evaluate Q and ϕ i , we minimized the energy E ¼ H h i using periodic boundary conditions along all directions while fixing ϕ 1 = 0. For the A.J. Clune et al.
incommensurate spin state, we minimized the energy within a magnetic unit cell with 2000 sites along c. After minimization, we checked that the classical forces on each spin vanish and we tried multiple Q and ϕ i values as initial values in the optimization to avoid metastable states. The detailed procedure to determine the field-dependent magnetization is described in section 3 of the Supplementary Information.

DATA AVAILABILITY
The data that support the findings of this study are available from corresponding authors upon reasonable request.