Linear magnetoelectric effect in göthite, α-FeOOH

By means of symmetry analysis, density functional theory calculations, and Monte Carlo simulations we show that goethite, α-FeOOH, is a linear magnetoelectric below its Néel temperature T N = 400 K. The experimentally observed magnetic field induced spin-flop phase transition results in either change of direction of electric polarization or its suppression. Estimated value of magnetoelectric coefficient is 0.57 μC · m−2 · T−1. The abundance of goethite in nature makes it arguably the most widespread magnetoelectric material.


Results
Goethite, α-FeOOH, crystallizes in the orthorhombic structure with space group symmetry Pbnm (Z = 4) shown in Fig. 1(a) and lattice parameters a = 4.55979 Å, b = 9.951 Å, and c = 3.0178 Å 16 . Upon decreasing temperature it experiences an antiferromagnetic phase transition at temperature T N , which varies in the range from approximately 340 to 400 K depending on the purity of the sample [17][18][19] . Below T N the spins → S i of four iron ions Fe i (i = 1, 2, 3, 4) located at positions (0.0489, 0.8537, 1/4), (0.9511, 0.1463, 3/4), (0.5489, 0.6463, 3/4), and (0.4511, 0.3537, 1/4) 16 , order antiferromagnetically with relative spin arrangement (+ − − +), respectively, as shown in Fig. 1(a) 20,21 . This ordered spin arrangement can be described by the order parameter → A . Other possible spin arrangements with → = k 0 described by the order parameters → F , → G , and → C are summarized in Table 1. The direction of the ordered spins is experimentally found to be along the c axis of the crystal cell. Therefore, the appearing ScIentIfIc REPORTS | 7: 16410 | DOI: 10.1038/s41598-017-16772-w magnetic structure with the wave vector → = k 0 can be described by the order parameter A z . Below we adopt an orthogonal system of axes x, y, and z being parallel to the crystal axes a, b, and c, respectively.
The symmetry of magnetic structure with A z ≠ 0 appearing below T N is Pb′nm 21 and allows linear magnetoelectric effect with magnetoelectric interactions given by z z y where → F and → P are ferromagnetic moment and electric polarization, respectively. Thus, in the antiferromagnetic phase magnetic field applied along the y or z axis induces electric polarization components P z or P y , respectively. It is found, however, that sufficiently strong magnetic field along the z axis results in a spin-flop transition, in which the spins reorient towards either the x or the y axis 18 . This will be discussed in more detail below.
It has to be noted here, that in the case when the initial paraelectric and paramagnetic phase possesses inversion symmetry operation, a magnetic phase transition with → = k 0 occurring according to a single irreducible representation cannot induce electric polarization 22 . However, linear magnetoelectric effect can be possible, as is the case in α-FeOOH: when A z ≠ 0 appears, the inversion symmetry is broken, but spatial inversion together with time reversal operation is a symmetry element, which results in interactions (1) and (2).
Using density functional theory we calculate six magnetic exchange constants, which are summarized in Table 2 and the respective exchange paths are shown in Fig. 1(b,c). It is found that the exchange couplings are mostly antiferromagnetic and the magnetic ground state is described by → ≠ A 0 in accordance with the experiments.    . The origin of frustration is in the presence of triangular arrangements of spins, which interact via the three dominant antiferromagnetic exchange couplings J 1 , J 2 , and J c , as shown in Fig. 1(c).
At H c = 20 T a spin-flop transition occurs in goethite 18 resulting in rotation of the antiferromagnetic vector to either a-or b-axis. The experimental value of the spin-flop magnetic field H c is used to scale our results on magnetic field dependence of magnetization in the antiferromagnetic phase shown in Fig. 2(c), which are in qualitative agreement with the experimental data 18 . In our Monte Carlo simulations we assume D x > D y , which results in appearance of A y at H z ≥ H c and corresponding vanishing of A z . Figure 2(d) shows H-dependence of electric polarization calculated using Eqs (1 and 2) and the ME interaction y z z which is relevant in the spin-flopped phase in the case when D x > D y . (As follows from Eqs (1-3) the components of electric polarization can be calculated multiplying the components of → A and → F obtained by Monte Carlo simulations). It follows that in the antiferromagnetic phase α-FeOOH is a linear magnetoelectric, since external H y and H z induce P z and P y , respectively. Furthermore, at H z = 20 T a flop of polarization from the b-to c-axis may occur. In the case D x < D y the antiferromagnetic vector will flip to A x at H z ≥ 20 T resulting in disappearance of electric polarization.
The microscopic origin of ME effect can be understood rewriting the ME interaction (1) through spins    The interaction w 1 is a single-ion contribution, whereas w 2 , w 3 , and w 4 result from interactions of two spins. Thus, the ME coupling may have both single-ion and two-ion contributions. The single-ion contribution is in accordance with the local non-centrosymmetric crystal environment of Fe atoms, the local crystal symmetry of which is a mirror plane σ z oriented parallel to the xy plane. Thus, it allows local spin-dependent electric dipole moments of electron orbitals d z ~ S y S z 23 .
In order to estimate the value of magnetically induced electric polarization we performed non-collinear DFT calculations. The spins were first relaxed in the stable A z configuration and then constrained to give additional ferromagnetic component F y . Artificially induced ferromagnetic ordering amounted to approximately 0.63 μ B per spin (while the antiferromagnetic component A z to approximately 4.09 μ B per spin), which resulted in rotation of spins away from the z axis by about 8.8°. The resulting electric polarization calculated using the Berry phase approach was found to take the value of 120 μC/m 2 . Taking the experimental magnetic susceptibility of approximately 0.003 μ B /T per Fe 3+ ion 24,25 we can estimate the ME coefficient to be of the order of 0.57 μC · m −2 · T −1 , which is comparable to that of LiNiPO 4 23,26 .
When D x > D y and the external magnetic field is higher than the spin-flop field H z > H c , the antiferromagnetic vector changes to A y ≠ 0 and P z appears. DFT calculations in this phase give value of the ME coefficient ∂P z /∂H z very close to the value of ∂P z /∂H y in the phase A z ≠ 0 calculated above. These ME coefficients together with the spin-flop field value allow to set scaling of P z and magnetic field in Fig. 2(d). However, the ME coefficient ∂P y /∂H z in the low field phase A z ≠ 0 depends on the unknown magnetic susceptibility χ = ∂ ∂ M H / z z , which should be much lower, though, than χ ⊥ = ∂M y /∂H y , precluding from setting reliable scale for P y in Fig. 2(d).
Relative values of different contributions to ME effect can be estimated from DFT calculations. For this purpose one can use the ME interactions Performing calculations using the magnetic configurations C z C y , A y F z , and C z C y similar to above and evaluating P z using the Berry phase approach we find that the biggest contribution to ME effect is w 3 and the other contributions relative to w 3 are w 1 /w 3 ≈ −0.034, w 2 /w 3 ≈ −0.21, and w 4 /w 3 = 0. Therefore, it follows that w 1 and w 2 act in the direction opposite to w 3 .

Conclusions
Based on the symmetry analysis of the available crystal and magnetic structures of goethite, α-FeOOH, we suggest that it is linear magnetoelectric below its Néel temperature. Using density functional calculations and Monte Carlo simulations we find main exchange constants in goethite and calculate its magnetic and magnetoelectric behavior.
Goethite belongs to the α-AlOOH diaspore structural type, which is also shared by, for example, α-MnOOH, Fe(OH)F, and Co(OH)F. The latter compound is also antiferromagnetic below ~40 K with the spin arrangement similar to α-FeOOH 27 and should, thus, display similar linear ME properties below its T N .
Nature creates beautiful polycrystalline goethite samples, which are encountered in significant amounts in various deposits. However, synthesis of single crystals in laboratory or preparation of good ceramic samples can be a challenge, as α-FeOOH starts to decompose at temperatures higher than 200 °C to form hematite, α-Fe 2 O 3 In this respect it may be easier to show the magnetoelectric behavior experimentally in the aforementioned isostructural compounds with similar magnetic structure, e.g., in Co(OH)F.

DFT calculations. Density functional theory calculations were performed using the Vienna Ab-initio
Simulation Package (VASP) 28 and the projected augmented wave method 29 . We used the GGA exchange correlation approximation corrected by means of the GGA + U formalism for the Fe atoms with U eff = U − J = 3 eV within the Dudarev approach 30 . This value of U eff was shown earlier to properly account for the structural and magnetic properties of α-FeOOH 31,32 . The energy cutoff was 850 eV, whereas the Brillouin zone integration was done using the 8 × 4 × 12 set of k-points determined by the Monkhorst-Pack scheme, which provided both the energy and k-points convergence.
Spin polarized collinear calculations were used to determine the exchange constants. For the determination of 6 exchange couplings the Hamiltonian was fitted to relative total energies of 7 different collinear magnetic structures. Magnetic cell sizes included a × b × c, 2a × b × c, a × b × 2c, and 2a × 2b × c, with the corresponding changes in the k-points grid. The crystal structure was relaxed in the most stable → A magnetic structure using the stopping criterion for absolute values of forces on atoms of 10 −3 eV/Å, while fixing the atomic positions for total energy calculations of other magnetic structures.
Electric polarization was calculated using the Berry phase approach as implemented in VASP while including spin-orbit coupling and performing fully non-collinear magnetic calculations. Since electric polarization appears in non-collinear magnetic structures the directions of local magnetic moments were constrained to form slightly non-collinear magnetic structure, while allowing for atomic and structural relaxation.
The calculated lattice parameters a = 4.638 Å, b = 10.037 Å, and c = 3.038 Å are within 1% of the experimentally determined values 16 where → S are classical vectors of unit length, D α (α = x, y, z) are anisotropy constants, and → H is magnetic field. In the Hamiltonian (12) we account only for the anisotropy terms D α ≠ 0 pursuing a minimal model, which reproduces the ground state A z ≠ 0. However, it has to be noted that there exist other anisotropy terms, which can generate some G-type contribution to the magnetic structure in the spin-flopped phase (i.e., when A x ≠ 0 or A y ≠ 0) when magnetic field is applied along the z direction, but not in the ground state, since A z is the only order parameter transforming according to the irreducible representation Γ 2− . We assume that such interactions are small and will only slightly modify the results quantitatively. In our simulations we tentatively use D x = −D z = 1.5 meV and D y = 0, which reflects the easy axis direction parallel to the c-axis.
The calculations are performed using the Metropolis scheme and a simulation box with dimensions 12 × 12 × 12 unit cells. After every change in temperature or external magnetic field the system is allowed to relax for 5·10 3 Monte Carlo steps per spin (MCS), whereas the statistical information is subsequently gathered over the next 15·10 3 MCS. Simulations of larger systems with dimensions 18 × 18 × 18 and 24 × 24 × 24 resulted in only slight increase of T N by approximately 2% and 3%, respectively.