K Beauvois
email: beauvois@ill.fr

V Simonet

S Petit

J Robert

F Bourdarot

M Gospodinov

A A Mukhin

R Ballou

V Skumryev

E Ressouche

Dimer Physics in the Frustrated Cairo Pentagonal Antiferromagnet Bi 2 Fe 4 O 9

The research field of magnetic frustration is dominated by triangle-based lattices but exotic phenomena can also be observed in pentagonal networks. A peculiar noncollinear magnetic order is indeed known to be stabilized in Bi2Fe4O9 materializing a Cairo pentagonal lattice. We present the spin wave excitations in the magnetically ordered state, obtained by inelastic neutron scattering. They reveal an unconventional excited state related to local precession of pairs of spins. The magnetic excitations are then modeled to determine the superexchange interactions for which the frustration is indeed at the origin of the spin arrangement. This analysis unveils a hierarchy in the interactions, leading to a paramagnetic state (close to the Néel temperature) constituted of strongly coupled dimers separated by much less correlated spins. This produces two types of response to an applied magnetic field associated with the two nonequivalent Fe sites, as observed in the magnetization distributions obtained using polarized neutrons.

Magnetic frustration, expected to occur when all spin pair interactions cannot be simultaneously satisfied, is one of the major ingredients at the origin of the flurry of discoveries in magnetic studies for the last 20 years. One of the basic experimental signatures of frustration is the difficulty of a system to order magnetically in spite of significant magnetic interactions, the extreme case being the absence of magnetic order at zero temperature as a consequence of the macroscopic degeneracy of the ground state [1][2][3][START_REF] Moessner | Introduction to Frustrated Magnetism[END_REF]. In materials where the magnetic moments order eventually at finite temperature, this leads in the region above the ordering temperature but well below the temperature characterizing the strength of the interactions, to a classical spin liquid state also called cooperative paramagnetic, where the magnetic moments are highly correlated although fluctuating. This disordered state can sustain not only well-defined excitations [START_REF] Robert | [END_REF] but also zero energy modes, which are the signature of local motions connecting the ground-state spin configurations. These modes acquire a gap on entering the ordered phase [6] and are alternatively described in a molecular approach [7]. This rich physics has been well established for triangle-based lattices, for instance the Heisenberg kagome antiferromagnet with nearest-neighbor interactions [START_REF] Robert | [END_REF]8]. Many other exotic manifestations of magnetic frustration have been revealed in quantum systems or with additional ingredients such as strongly anisotropic Hamiltonians [9][10][11].

Another direction has been opened with the identification of an equivalent for the pentagonal Cairo lattice in the real material Bi 2 Fe 4 O 9 [12]. The Cairo lattice is not based on triangles but on edge-sharing pentagons. It is then still prone to magnetic frustration due to the odd number of bonds in the elementary pentagonal units. This pentagonal lattice has a complex connectivity with three-and fourfold connected sites at variance with triangle-based lattices, which fosters alternative ways to accommodate frustration. This leads in Bi 2 Fe 4 O 9 to an unconventional ground state consisting of an orthogonal arrangement of the magnetic moments. This classical ground state was also obtained theoretically, as well as other interesting phases, including, in the presence of quantum fluctuations, a resonating valence bond liquid or an orthogonal dimer ground state (valence bond crystal) [13,14]. The latter recalls the exactly solvable dimer ground state of the Shastry-Sutherland lattice [15] largely investigated for its exotic physics [16,17]. These findings stimulated further theoretical [18][19][20][21][22][23][24] and experimental [25][26][27][28] studies on pentagon-based physics, even spreading beyond the field of magnetism. In spite of this interest, an experimental determination of the Hamiltonian of the prototypical material Bi 2 Fe 4 O 9 has never been reported yet that would solely ascertain the crucial influence of frustration on its exotic properties.

In this Letter, we present the experimental determination of the magnetic interactions in Bi 2 Fe 4 O 9 materializing a Cairo lattice with spins 5/2 using inelastic neutron scattering. A minimum of five exchange interactions allows us to account for the magnetic order and for the as-sociated excitations including a peculiar quasiflat mode. We also show the magnetic density maps measured using polarized neutron scattering under a magnetic field. They reveal a correlated state above the ordering temperature, resulting from the hierarchy of interactions and characteristics of the underlying dimer physics. The unit cell of the orthorhombic oxide Bi 2 Fe 4 O 9 contains eight magnetic Fe 3+ ions equally distributed on two different Wyckoff sites of the P bam space group: 4h for Fe 1 and 4f for Fe 2 . These sites have a different connectivity and different oxygen coordination, tetrahedral for Fe 1 and octahedral for Fe 2 . They form a lattice closely related to the Cairo pentagonal one with noticeable differences (see Fig. 1): the site with fourfold connectivity in the perfect lattice is actually constituted by a pair of Fe 2 atoms located below and above the pentagonal plane. The projected lattice in the ab plane is also slightly distorted compared to the perfect Cairo lattice as far as the bond lengths and bond angles are concerned. The Heisenberg Hamiltonian with isotropic exchange interactions is a good starting point to describe the magnetic properties of this frustrated lattice of Fe 3+ ions (J = S = 5/2, L = 0)

H = i,j J ij S i • S j (1)
where S i is a spin operator and J ij = J 1 -J 5 are the five superexchange interactions between pairs of spins, inferred from the structure [12]. There are three exchange couplings in the ab plane (J 3 , J 4 , and J 5 ) and two additional ones out of the plane of the pentagons (J 1 and J 2 ), which connect only the stacked Fe 2 [see Fig. 1 (b)]. Bi 2 Fe 4 O 9 exhibits below T N ∼ 240 K a long-range antiferromagnetic order characterized by a propagation vector k = (1/2, 1/2, 1/2). The resulting spin configuration was elucidated by neutron diffraction (see Fig. 1) and is made of two sets of orthogonal pairs of antiferromagnetic spins in the ab plane corresponding to each Fe site, with a global rotation by an angle α = 155 • relative to each other [12]. Despite the deviation of the experimental system from the perfect Cairo lattice, this orthogonal magnetic structure matches the one identified in theoretical studies on the perfect lattice and is shown to be rather robust while varying the ratio of the exchange interactions and increasing the quantum fluctuations [14]. This peculiar magnetic structure actually results from both frustration and complex connectivity.

All our experiments were performed on a single crystal of Bi 2 Fe 4 O 9 of dimensions ∼ 2.5 × 2 × 1.5 mm 3 , grown by the high flux temperature solution method using a flux of Bi 2 O 3 . The magnetization distributions were obtained from two neutron scattering experiments performed on the CRG D23 two-axis diffractometer at Institut Laue Langevin (ILL) in its polarized neutron mode, with an incoming neutron wavelength λ = 2.37 Å from a graphite-Heusler double monochromator configuration. Measurements were performed with a magnetic field of 6 T applied along the sample c axis at T = 250 K (paramagnetic state) and along the a-b crystallographic direction at T = 250 K and 15 K (ordered state) (see Supplemental Material [START_REF]See supplemental material for details on the magnetization density maps, spin waves analysis, comparison with ab initio calculations, quasiflat excitation, and[END_REF]). Inelastic neutron scattering (INS) experiments were performed on the CRG IN22 triple-axis spectrometer at ILL in an orange cryostat at a constant final wave vector k f = 2.662 Å-1 , with an energy resolution equal to about 1 meV. The sample was oriented in order to access the (h, h, ) scattering plane.

In order to measure the magnetic excitations, several constant-Q energy scans have been performed at 1.5 K along the reciprocal space directions (1/2, 1/2, ), (h, h, 0) and (2, 2, ) sketched in Fig. 2(a). The excitations along (2, 2, ) were fitted using the takin software [35,36], which includes the spin wave model as input and the instrument resolution [see Fig. 2(b)] (see Supplemental Material [START_REF]See supplemental material for details on the magnetization density maps, spin waves analysis, comparison with ab initio calculations, quasiflat excitation, and[END_REF]). This treatment was necessary to disentangle two overlapping modes. The excitations in the other two directions were simply fitted by the sum of two Lorentz functions,

S(Q, ω) = bg+A(ω, T ) 1 (ω -) 2 + σ 2 + 1 (ω + ) 2 + σ 2 (2) where bg is the background, A(ω, T ) = [1 + n(ω, T )]Zωσ/2, n(ω, T ) = 1/(e ω/k B T -1)
is the Bose factor, Z is the weight of the excitations, is their energy and σ is the half width at half maximum. The resulting energy position of the excitations is reported as black dots in Figs. 2(c)-2(h). All Q scans have been combined into the experimental maps presented in Figs. 2(c)-2(e). The magnetic nature of the excitations has been checked through the temperature dependence. Well-defined spin waves are observed, as expected for this ordered compound. An acousticlike mode emerging from the antiferromagnetic Bragg position (1/2, 1/2, 1/2) is clearly visible, as well as a high-energy branch along (h, h, 0). Additionally, an almost nondispersive mode located at the energy of about 19 meV is observed along (2, 2, ). Higher energy modes are inferred from Raman spectroscopy with two magnetic excitations identified at the energies 32.2 and 58.5 meV, out of the energy window of our neutron experiments [37].

Our INS measurements were then compared with spin wave calculations performed using the SpinWave software [38,39] based on the linear spin wave theory using the Holstein-Primakoff formalism [40]. The starting point was the model hamiltonian of Eq. 1 involving five isotropic superexchange interactions (see Fig. 1). Two additional constrains were used to limit the number of refined parameters. First, it was shown in 12 that the rotation angle α =155 °between both iron sublattices is obtained for J 3 /J 5 = 2.15. A second relation between the exchange interactions was inferred from the Curie-Weiss temperature θ CW ≈ -1670 K estimated from magnetic susceptibility measurements [12]. We used the local Weiss molecular field model on the stabilized magnetic structure and the equipartition theorem [START_REF] White | Quantum Theory of Magnetism[END_REF] 2

× 3 2 k B θ CW = S 1 j J 1j < S j > +S 2 j J 2j < S j >, (3)
where S 1 (S 2 ) is the spin on site Fe 1 (Fe 2 ). This allowed us to further reduce the model to three independent parameters, which were systematically varied in the calculations. We checked the capability of each of the sets of parameters to reproduce the measured spin waves, as well as the magnetic structure of Bi 2 Fe 4 O 9 through a real-space mean-field energy minimization of the spin configuration. Finally, a model Hamiltonian compatible with the experiments was obtained with the values of the exchange constants given in Table I. The calculated spin waves are displayed in the lower panels of Figs. 2(f)-2(h) and show a very good agreement with the experimental data. The calculations indicate that the spin waves dispersion extends at higher energies up to 80 meV (see Supplemental Material [START_REF]See supplemental material for details on the magnetization density maps, spin waves analysis, comparison with ab initio calculations, quasiflat excitation, and[END_REF]), with, in particular, two zone center modes at the energy positions of the Raman excitations, which further validates our model.

Our analysis establishes that all five interactions are antiferromagnetic. Because of the geometry of the pentagonal lattice, this implies competition between the J 3 , J 4 and J 5 interactions within the ab planes confirming the role of magnetic frustration in the stabilization of the 90 • magnetic order previously reported [12].

The antiferromagnetic interaction J 1 between the Fe Along (1/2, 1/2, ) and (h, h, 0), this mode was not visible in our experimental data. However, we investigated in purpose the (2, 2, ) direction where this mode was strongest in the calculations. We indeed observed it experimentally along this reciprocal space direction. This excited state corresponds to the out-of-phase precession of both Fe 2 spins around their local field, which does not influence the neighboring Fe 1 spins and thus remains localized. A weak dispersion is actually calculated due to the J 2 interplane couplings in the c direction, which is beyond the instrumental resolution of the experiment (see Supplemental Material [START_REF]See supplemental material for details on the magnetization density maps, spin waves analysis, comparison with ab initio calculations, quasiflat excitation, and[END_REF]). Also noticeable in Table I is the fact that the J 4 interaction is significantly stronger than the other ones, which is compatible with the 180 • superexchange path through the central oxygen ion according to the Goodenough-Kanamori rules [START_REF] Goodenough | Magnetism and the Chemical Bond[END_REF]. This hierarchy of interactions results in a lattice with dominant pairs of antiferromagnetically coupled spins Fe 1 on almost orthogonal bonds [see Fig. 1(a)], a picture that is expected to survive above the Néel temperature. Note that the exchange interactions deduced from our spin wave analysis are slightly different from those obtained from ab initio calculations [START_REF]See supplemental material for details on the magnetization density maps, spin waves analysis, comparison with ab initio calculations, quasiflat excitation, and[END_REF][START_REF] Pchelkina | [END_REF] in the LSDA+U approximation. However, the positive sign of all exchange couplings (antiferromagnetic) and the overall hierarchy of the interactions are identical in both cases. In order to investigate the fingerprint of these Fe 1 dimers in the paramagnetic state, we measured the magnetization distributions under a magnetic field of 6 T with two orientations with respect to the crystallographic axes, as shown in Fig. 3. The magnetization per site has been extracted from these maps in the dipolar approximation implying a spherical electronic distribution around the atoms. Several fitting processes were performed taking into account the presence of magnetic moments on the iron sites only, or on all the iron, oxygen and bismuth sites. A magnetic contribution can indeed be present on the oxygen atoms and, due to their 6s lone electron pairs, on the Bi atoms. The final averaged magnetizations per Fe sites are given in Table II.

Interestingly, in the paramagnetic state [Figs. 3(a) and 3(b)], the two iron sites have radically different behaviors: whereas the Fe 2 ions carry an induced magnetic moment of 0.045(5) µ B aligned along the field, the induced magnetization on site Fe 1 is vanishingly small (see Table II). This is in contrast with the ordered magnetic moments on both sites refined from previous neutron diffraction experiments below T N which are rather similar, equal to 3.52 and 3.73 µ B respectively [12]. Our measurements performed with two different directions of the magnetic field yield the same result, which show that the anisotropy is not responsible for this behavior as expected for Fe 3+ ions with zero orbital angular momentum. The most noticeable difference between the two maps in Figs. 3(a) and 3(b) is the presence of magnetic density on the Bi sites for H c and not for H ab. This could be real or an artifact due to an imperfect reconstruction for H c since a smaller number of observations has been used compared to H a-b. The magnetization distribution has also been measured at low temperature (see Fig. 3(c)). A larger field-induced polarization on the Fe 2 compared to the Fe 1 one is actually preserved in the ordered state. As detailed in Table II, the magnetization on the Fe 1 is actually slightly larger at 15 K than at 250 K and also more delocalized (see Supplemental Material [START_REF]See supplemental material for details on the magnetization density maps, spin waves analysis, comparison with ab initio calculations, quasiflat excitation, and[END_REF]).

These results point to an original paramagnetic state. A calculation assuming free spins 5/2 in a 6 T field at 250 K yields a field-polarized magnetization of 0.19 µ B , which is 4 times larger and almost 200 times larger than the one measured on the Fe 2 and Fe 1 respectively. This suggests, that slightly above T N , the Fe 1 form an assembly of strongly correlated antiferromagnetic dimers in agreement with the dominant J 4 , while the Fe 2 spins are much less correlated. Since the Fe 1 ions form a spin arrangement with great similarities to a Shastry-Sutherland lattice, we suggest that the temperature regime slightly above T N could be reminiscent of this original physics [15][16][17]. At higher temperature, the correlations among dimers should vanish, while at lower temperature, correlations involving Fe 2 spins grow, driving the system to the physics of the Cairo pentagonal lattice. Below T N , the long-range magnetic order is finally triggered by the weakest J 2 interaction, connecting the pentagonal planes. Note that Bi 2 Fe 4 O 9 is not the unique materialization of the Cairo lattice as it is related to a wide family of compounds including the multiferroic RMn 2 O 5 (R a rare-earth/Y, Mn occupying the pentagonal lattice) whose complex magnetodielectric phase diagrams could be investigated in the renewed perspective of pentagonal physics [27,44,45].

Our neutron scattering investigation of Bi 2 Fe 4 O 9 allowed us to achieve a complete determination of its complex magnetic interactions and to unveil various facets of unconventional magnetism including frustration and dimer physics, with distinct behaviors associated with the two inequivalent Fe sites of the pentagonal lattice. The Fe 1 ions produce strongly coupled antiferromagnetic pairs of spins dominating the correlated paramagnetic state, whereas in the ordered state, the pairs of Fe 2 spins produce original spin dynamics, associated with protected local motions, coexisting with dispersive spin waves. Beyond the canonical examples of frustrated systems like kagome or pyrochlore lattices with firstneighbor interactions, our Letter discloses novel behaviors that should be more generally observed in materials where the frustration is interlocked with complex connectivity and hierarchal interactions.

To obtain the magnetic density maps, the flipping ratio method was used. Bragg reflexions were collected up to sin θ/λ = 0.36 Å-1 and the intensities were reduced and corrected for extinction. The polarization for the incident neutron beam was 0.90 (1). Accurate magnetic structure factors were extracted from the flipping ratios using the Cambridge Crystallography Subroutine Library [S1]. The magnetization distribution in real space was reconstructed using a uniform prior density and the Maximum of Entropy technique [S2] in order to minimize the noise and truncation effects. 28 and 49 independent magnetic structure factors F M (Q) were used in the reconstruction for H c and H ab respectively, as well as a value of the macroscopic magnetization measured with an extraction magnetometer at Institut Néel. In order to fit the magnetic excitations along (2, 2, ), we used the open-source software, Takin, for triple-axis experiments analysis [S3, S4]. We simulated the experimental spectra by convoluting a theoretical dynamical structure factor S(Q, ω) model with the instrumental resolution of the IN22 triple axis spectrometer. The resolution function has been calculated using the Popovici algorithm. We used as theoretical model two simple dispersions, one sinusoidal and one flat, for excitations with a lorentzian shape of full width at half maximum equal to 4 meV. The input of the software were the spin wave dispersions and intensities calculated with the SpinWave program allowing a first fit for each constant-Q ω-scan of the global intensity scale factor and of the offset in ω only. In a second step, all the parameters were let free to vary and a fit of the dispersions, intensities and width of the excitations was performed. This procedure allowed us to distinguish the quasiflat mode from the lower dispersive mode and to determine precisely its energy position at about 20 meV. The excitations in the other two directions (1/2, 1/2, ), (h, h, 0) were simply fitted by the sum of two Lorentz functions as explained in the main paper.

III. ADDITIONAL INFORMATION ABOUT THE QUASIFLAT EXCITATIONS

As explained in the main text, the spin excitation spectrum includes a peculiar nearly flat optical mode. It involves deviations of the Fe 2 spins from their equilibrium position, that remain decoupled from the fluctuations of the four Fe 1 to which they are in principle connected in triangular pathways (see Figure S5). For the sake of simplicity, we consider in this reasoning a single pair (S 1i , S 2i ) in each unit cell i. We assume this pair is located on the y axis, hence the index i refers to the x axis. The case of the second pair (along x) can be treated in the same way, the index i being associated to the y axis. Furthermore, we set J 2 = 0.

For the purpose of this derivation, we define the equilibrium and fluctuations as: S 1i = S o 1i + s 1i . The local field H i experienced by the spins of the same pair is also written: H i = H o i + h i where H o i is the equilibrium mean value and h i the corresponding fluctuations.

Owing to the magnetic structure, we have: S o 1i = S o 2i = S o (the two Fe 2 are parallel) and

H o i = H o i-1 =
H o (according to the propagation vector, the spins of the Fe 1 pairs belonging to neighboring cells and connected to a given Fe 2 spin are parallel). Up to first order, the equations of motion write:

ds 1i dt = (H o i × s 1i + H o i-1 × s 1i + h i × S o 1i + h i-1 × S o 1i ) + J 1 S 2i × S 1i ds 2i dt = (H o i × s 2i + H o i-1 × s 2i + h i × S o 2i + h i-1 × S o 2i ) + J 1 S 1i × S 2i
We then define auxiliary variables:

D i = S 1i -S 2i = s 1i -s 2i S i = S 1i + S 2i = 2S o + s 1i + s 2i
hence:

dD i dt = d(s 1i -s 2i ) dt = (H o i × D i + H o i-1 × D i + (h i + h i-1 ) × (S o 1i -S o 2i )) + 2J 1 S 2i × S 1i = 2H o × D i + 0 + 2J 1 S 2i × S 1i = 2H o × D i + J 1 S i × D i = (2H o + 2J 1 S o ) × D i
This equation shows that D i is decoupled from the H and S fluctuations and only feels the effective magnetic field (2H o + 2J 1 S o ). As a result, the spectrum encompasses an infinite number of resonances, (with no dispersion), whose energy essentially depends on the molecular field (2H o + 2J 1 S o ). Similarly, we can derive the equation of motion for S i , leading to the acoustic counterpart:

dS i dt = d(s 1i + s 2i) dt = (H o i × s i + H o i-1 × s i + (h i × +h i-1 ) × (S o 1i + S o 2i )) + 2J 1 S 2i × S 1i = 2H o × s i + 2(h i + h i-1 ) × S o
For J 2 = 0, the calculations performed with the SpinWave software [S5, S6] show that the dispersion of this mode only depends on J 2 , connecting the pentagonal planes. For the sake of illustration, Figure S6-b displays the dispersion calculated along (2, 2, ) for several values of J 2 ranging from J 2 = 0.1 up to J 2 = 3.1 meV. Importantly, with increasing J 2 , the bandwidth of the mode increases and finally connects to the acoustic modes stemming from zone centres. Indeed, the localized character of this mode is ensured by small values of J 2 . Figure S6-c shows constant energy cuts (at 5, 10, 15 and 19 meV) across the dispersions. Classical cones are observed stemming from the magnetic Bragg peaks, corresponding to the acoustic modes. Those modes become essentially 1-dimensional, forming lines along . The flat mode has a peculiar dynamical structure factor, which can be deduced from the 19 meV cut. As expected, it essentially depends on , with stronger intensities at = 0, 1, 2, ...

V. CLASSICAL SPIN DYNAMICS CALCULATIONS OF THE MAGNETIC STRUCTURE FACTOR

We have computed the dynamical structure factor through classical spin dynamics simulations. More specifically, S(Q,ω) was solved numerically by combining a Monte Carlo method to an integration of the equations of motion, therefore taking into account non-linear effects associated with thermal fluctuations (see Ref. [S10-S12] for technical details about the method). The calculations have been performed using the Hamiltonian determined by our spin wave analysis. Intensity maps of S(Q,ω) are shown in Figure S8 for different relative temperatures. The evolution of the structure factor highlights strong spin pair correlations persisting even at temperature much higher than the Néel temperature and several regimes in the paramagnetic state. Just above the Néel temperature, the Fe 1 and Fe 2 are correlated within the pentagonal planes but the correlations along the c axis through J 2 is lost. At higher temperatures, typically 7 × T N , only the Fe 1 dimers remain correlated forming a classical dimer state recalling the Shastry-Sutherland physics, but where the Fe 1 quasi-orthogonal dimers are separated by paramagnetic free F e 2 spins. Finally, the system is expected to become completely uncorrelated at very high temperature. The dimer regime is identified from a comparison with a calculation taking into account only the Fe 1 spins, which are responsible for the high-energy features, as shown in Figure S9.

FIG. 1 .

 1 FIG. 1. (a) ab plane projection of the Fe atoms in Bi2Fe4O9 forming an equivalent of the Cairo pentagonal lattice. The Fe1 (in blue) and Fe2 (in orange) have a different connectivity. Blue ellipses underline strongly coupled antiferromagnetic Fe1 spins. (b) Magnetic arrangement stabilized below TN : the orange and blue rectangles materialize the two sets of orthogonal antiferromagnetic pairs in the ab plane [12]. The Fe2 form ferromagnetic pairs sandwiching the Fe1 planes. The five exchange interactions within and between the pentagonal planes are labeled.

2

 2

FIG. 2 .

 2 FIG. 2. (a) Sketch of the scattering plane (h, h, ) investigated by INS in Bi2Fe4O9, with the nuclear (blue spheres) and magnetic (red spheres). Bragg peak positions and the measurement cuts along the (1/2, 1/2, ), (h, h, 0) and (2, 2, ) reciprocal space directions. (b) Along (2, 2, ), series of measured constant-Q energy scans and fit of the excitations (black lines). Measured (c)-(e) and calculated (f)-(h) dynamic structure factor S(Q, ω) using the exchange constants of Table I and a small single-ion anisotropy term constraining the spins in the ab plane. The color scale of the calculations was truncated for the spectra along (1/2, 1/2, ), (h, h, 0) in order to emphasize the weaker flat mode along (2, 2, ). The width of the calculated excitations was taken as the energy resolution of 1 meV. The empty (filled) black points on top of the measured (calculated) S(Q, ω) give the fitted energy positions of the experimental spin wave dispersion.

  Values of the magnetization measured in Bi2Fe4O9 at T = 250 and 15 K and under a magnetic field µ0H = 6 T: macroscopic magnetization m per unit cell (second column) obtained with an extraction magnetometer with the field direction along a-b at 15 K and along a-b or c at 250 K; magnetization per site mFe 1 and mFe 2 (third and fourth columns) deduced from the flipping ratio neutron experiments. At T = 15 K the magnetic density was obtained for H a-b. At T = 250 K, the magnetic density is the average of the values obtained for H c and H a-b since the same macroscopic magnetization is measured.

Fe2FIG. 3 .

 3 FIG. 3. Spatial distributions of the magnetization of Bi2Fe4O9 in µB projected along the c axis measured with polarized neutrons under a magnetic field µ0H = 6 T applied (a) along c at T = 250 K and (b) along a-b at T = 250 K and (c) along a-b at T = 15 K (c).

FIG. S2 .FIG. S3 .

 S2S3 FIG.S1. Same spatial distribution density maps of Bi2Fe4O9 as those presented in the main paper (Figure3) but drawn with a different intensity scale for the magnetization in order to emphasize the weak features. This highlights the delocalization of the magnetization in the ordered state compared to the paramagnetic state.

B

  FIG. S4. Calculated dynamic structure factor of Bi2Fe4O9 up to 75 meV along the reciprocal space directions (1/2, 1/2, ), (h, h, 0) and (2, 2, ) showing the complete magnetic spectrum. This was obtained using the SpinWave software [S5, S6] with the five exchange interactions reported in Table I of the main article and a single-ion anisotropy within the ab-plane of 0.03 meV. The red crosses indicate the two magnetic excitations identified in Raman scattering [S7].

FIG. S5 .

 S5 FIG. S5. Schematic description of the frustrated mechanism responsible for the quasiflat mode observed at 19 meV in the inelastic spectrum. It involves the out-of phase precession of both Fe2 spins of one ferromagnetic pair connected via a triangular exchange path to four Fe1 neighbors.

  FIG. S8. Calculated S(Q,ω) by classical spin dynamics for different relative temperatures T /TN through (a) Q-cuts and (b) constant ω-cuts, including zero energy ω0 and the energy corresponding to the nearly flat mode ω1 (indicated by the arrows). The intensity scale is the same for all the (a) maps and varies for the (b) maps in order to emphasize the main features.

2 FIG. S9 .

 2S9 FIG. S9.Comparison of the calculated S(Q,ω) for Bi2Fe4O9 using the set of exchange coupling obtained from our spin waves analysis with the one taking into account only the F e1 spins coupled via J4, at T /TN 7 through (a) Q-cuts and (b) constant ω-cuts over an energy range δω2 (between the dotted lines). The intensity scale is the same for all the (a) maps and varies for the (b) maps in order to emphasize the main features.

ACKNOWLEDGMENTS

We are grateful to J. Debray for the orientation of the sample, the Institut Laue Langevin for providing us with the neutrons, and Navid Qureshi for his help in the data analysis. J. R. thanks B. Canals for the joint development of the software used for the calculations of the spin dynamics presented in the Supplemental Material. Part of this project was supported by Bulgarian National Science Fund BNSF DN -08/9.

Supplemental Material: Dimer Physics in the Frustrated Cairo Pentagonal Antiferromagnet Bi 2 Fe 4 O 9 K. Beauvois, 1,2, * V. Simonet, 2 , S. Petit, 3 J. Robert, 2 F. Bourdarot, 1 M. Gospodinov, 4 A. A. Mukhin, 5 R. Ballou, 2 V. Skumryev, 6,7 and E. Ressouche We discuss the difference mentioned in the main text between the exchange interactions deduced from our spin wave analysis of the magnetic spectra measured by inelastic neutron scattering and the ones obtained from ab-initio calculations [S8]. In Table I, one can see that the values for the interactions are in the same energy range. There are however some discrepancies and, in particular, our J 4 value is twice larger. Moreover, the ab-initio J 3 /J 5 ratio does not yield the correct angle between the Fe 1 and Fe 2 magnetic sublattices. Figure S7 shows the calculated spin waves from the ab-initio interaction values. Contrary to the spin waves calculated with the set of values deduced from our spin wave analysis, the energy and shape of the dispersing branches do not match well with the experimental results. Due to the lower J 1 value, the energy position of the quasiflat mode is too high, at about 30 meV instead of 19 meV. These discrepancies might be accounted for by some limiting assumptions in the ab-initio calculations that are recalled in the reference [S9].


Xavier Fabrèges

Sylvain Petit

J.-B Brubach

Pascal Roy

M Deutsch

A Ivanov

L Pinsard-Gaudart

V Simonet

R Ballou

Sophie De Brion
email: sophie.debrion@neel.cnrs.fr

A Pimenov

A A Mukhin

V Yu Ivanov

V D Travkin

A M Balbashov

Interplay between spin dynamics and crystal field in the multiferroic compound HoMnO 3

Keywords: THz spectroscopy, inelastic neutron scattering, dynamical magnetoelectric effect

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

INTRODUCTION

Multiferroics continue to attract much attention in condensed matter physics, because of the promising applications in spintronics they offer, but also because these materials raise fundamental questions. One of these concerns the understanding of the microscopic origin of the coupling between magnetic and electric degrees of freedom, and how new excitations resulting from this coupling, for instance electromagnons, can emerge. Electromagnons were first reported in 2006 as magnetic excitations excitable by the electric field component of the electromagnetic wave [1]. Since then, electromagnons have been observed in several compounds with non colinear magnetic order such as perovskites [1, 2], hexaferrites [3] or even the more simple cupric oxide [4]. Various mechanisms at the origin of these excitations were invoked and identified, among which the inverse Dzyaloshinskii-Moriya interaction [5], symmetric exchange-striction mechanism [6] or higher harmonics of the magnetic cycloids [7]. Other types of excitations involving electric and magnetic degrees of freedom have also been reported such as atomic rotation mode excited by the magnetic field component of the terahertz (THz) light [8] or hybrid excitations arising from the interplay between rare-earth crystal field (CF) excitations and 3d transition metal magnons [9,10], in particular in the hexagonal manganite ErMnO 3 [11]. All these hybrid excitations often lie in the THz range and are also investigated using inelastic neutron [12] or X ray scattering [13] and Raman spectroscopy [14].

The hexagonal manganites h-RMnO 3 where R= Y, Er, Ho, etc, crystallizing in the P6 3 cm space group, have been widely studied [15][16][17][18] since they combine remarkable properties associated to multiferroism and magnetic frustration. They belong to type I multiferroics, where ferroelectricity and magnetic orders occur at different temperatures: onset of an electric polarization along the c-axis at much higher temperature (T C 800K) than the ordering of the Mn magnetic moments (T N 80K). Despite these different energy scales, non negligible coupling between the electric and magnetic degrees of freedom are present. For instance, the Mn magnetic ordering is accompanied by ionic motions inside the unit cell [19,20] giving rise to a giant magneto-elastic coupling, connected with an increase of the ferroelectric polarization. The Mn sublattice consists of triangular planes (Fig. 1) with a peculiar stacking along the c-axis, each Mn 3+ ion sitting close to the middle of the triangle of Mn 3+ belonging to adjacent layers. The actual materials have then in common a 120 • spin structure for the Mn magnetic moments in the triangular planes, but the stacking and relative phase difference of these patterns depend on the R atom ((magnetic R=Ho, Yb, or not R= Y, Sc), and correspond to the four one-dimensional irreducible representations shown in Fig. 1. A key parameter tuning the sign of the exchange interaction and the Mn 3+ single-ion anisotropy seems to be the shift away from 1/3 of the x position of the Mn 3+ ions within the triangular plane [19,21]. The ions R 3+ are situated on two Wyckoff sites 2a and 4b. For magnetic rare earths, they also order at lower temperature, either polarized by the Mn molecular field (site 4b) or from their mutual interaction in a temperature range of few K (site 2a) [22,23].

In this study, we focus on the Ho compound, which has been intensively studied for its rich H -T phase diagram and its strong magnetoelectric effects [17,24]. It becomes ferroelectric at 875 K and the Mn magnetic moments order at T N =75 K in the Γ 4 representation (see Fig. 1). Remarkably, a spin reorientation occurs at T SR =37 K which is concomitant with the ordering of the Ho(4b) magnetic moments, while the Ho(2a) moments order around T Ho =5 K, both along the c-axis [19,[25][26][27][28][29]. This Mn reorientation transition, unique in the hexagonal manganites, shows up also in the electric properties with a drop in the electric polarisation below T SR which is further recovered at T Ho [30]. Signatures of the rareearth Mn coupling in the spin dynamics in an applied magnetic field have also been reported [31,32], although no rare-earth CF-Mn magnon hybridization has yet been observed.

We report here a detailed study on the dynamical properties of this compound in the THz range using both electromagnetic waves and neutrons as complementary probes. We show that the THz dynamics is particularly complex in HoMnO 3 due to the presence of several Ho CF excitations in addition to the Mn spin waves that surprisingly become non-dispersive along the c-axis at low temperature. To interpret the observed spectra and their remarkable temperature dependence through the spin reorientation transition and below, we have performed calculations both for Mn and Ho dynamics. Our work brings new insight on the complex spin dynamics of this compound including possible signatures of the couplings between the different magnetic ions, but also with their lattice surroundings.

EXPERIMENTAL DETAILS

Single crystals have been grown using the floating zone method, yielding rods of about 7 mm in diameter and 4 cm long. From these crystals, two small plaquettes were cut for the THz measurements (Typically 2 mm in diameter and 300 µm in thickness): one with the c-axis perpendicular to the plaquette, the other one with the c-axis within the plaquette. The polarized THz emission of the synchrotron at SOLEIL was used to probe the THz properties similarly to the study on ErMnO 3 [11].

Inelastic neutron scattering (INS) measurements were performed at the thermal spectrometer IN8 (ILL) as well as at the 4F1 and 4F2 cold triple-axes spectrometers installed at LLB. The final wave vector used was k f = 2.662, 1.55 or 1.3 Å-1 depending on the desired energy resolution (0.8, 0.25 and 0.12 meV respectively). Higher order reflections were removed using a PG (IN8) or nitrogen cooled Beryllium filter (4F). The sample was mounted to have access to Q = (h, 0, ) scattering wavevectors and attached to the cold finger of an orange cryostat. Several Energy-scans taken along Q = (h, 0, ) were collected to construct the maps shown in Fig. 2 and3.

INELASTIC NEUTRON SCATTERING

The magnetic excitations have been first investigated by INS. Figures 2a and3 give an overview of those results, showing the spectra along a * , along c * , and the temperature dependence measured at the zone centers (1, 0, 0) and (1, 0, 1). Dispersive spin waves are clearly visible along a * with a maximum of the branches around 17 meV (140 cm -1 ) while the amplitude of dispersion along c * is much weaker, reflecting the hierarchy of the magnetic interactions much stronger within the ab-planes than in between. In addition to the spin waves, nondispersive excitations are visible corresponding to transition between CF levels of the two Ho 3+ ions on 4b and 2a sites. Remarkably, a significant temperature dependence of the magnetic excitations is visible at T SR as shown at the zone centers (1, 0, 0) and (1, 0, 1) in Fig. 3 where changes of spectral weight are observed for all the excitations below 30 cm -1 . More surprisingly, along c * , the spin waves dispersion disappears gradually below T SR and finally the corresponding branches remain flat below ≈ 10 K.

The coexistence of both spin waves and CF transitions, along with a subtle temperature evolution, make however the interpretation of those spectra quite complex. Fortunately, previous studies [12,19,22,[33][34][35][36][37] performed for several members of the RMnO 3 provide valuable assistance. They especially concluded to the relevance of the I. This is responsible for the rise of the small gaps at zone centers. All the spin wave calculations have been performed using the SpinWave software [38].

following Hamiltonian to describe the Mn spin waves:

H = 1 2 p, i,j J i,j S i,p S j,p + 1 2 i,j K i,j S i,1 S j,2 + i,p D z (S z i,p ) 2 -D x (S x i,p ) 2 -D y (S y i,p ) 2 (1)
In this expression, S i,p is the Mn spin at site i, p = 1, 2 is an index that labels the two Mn planes in the unit cell, D z is an easy plane anisotropy (D z ≥ 0), and D x,y ≥ 0 are effective easy axis anisotropies, forcing the spins to point along the x or y directions (see Fig. 1). J i,j is the antiferromagnetic coupling between neighboring spins in the same plane, and K i,j the antiferromagnetic coupling between neighboring spins in adjacent planes. As shown in Fig. 1a, the x M n position allows to distinguish two different paths for K, a long and a short one, labeled hereafter as K 1 and K 2 . To proceed, we first assume D x = 0 and D y = 0. Straightforward classical minimization of H shows that the ground state structures are planar 120 • spin configurations. Using the convention of Figure 1, the moments are given by :

S i,p = S   cos(φ p + (i -1) 2π 3 ) sin(φ p + (i -1) 2π 3 ) 0   (2)
and the classical energy writes:

E = -3JS 2 + 2(K 1 -K 2 )S 2 cos ∆φ, with ∆φ = φ 2 -φ 1 (
3) The value of φ 1 along with the relative phase difference ∆φ = φ 2φ 1 between two adjacent planes fully determine the final structure. These structures can be classified using the Γ i=1,2,3,4 irreducible representations listed in Table I. Note that two 2-dimensional irreducible representations are also allowed in the P6 3 cm space group but were not considered in this work. The sign of K 1 -K 2 determines the stability of each phase: Γ 1,2 are stabilized for K 1 -K 2 ≥ 0, while Γ 3,4 are stabilized for K 1 -K 2 ≤ 0.

In principle, six spin wave branches are expected in this model. However, since K 1 ≈ K 2 , these branches should be considered as three almost doubly degenerate

Γ i φ ∆φ Conditions 1 π/2 π D y ≥ D x , K 1 -K 2 ≥ 0 2 0 π D y ≤ D x , K 1 -K 2 ≥ 0 3 0 0 D y ≤ D x , K 1 -K 2 ≤ 0 4 π/2 0 D y ≥ D x , K 1 -K 2 ≤ 0 TABLE I.
Condition for the Mn 3+ single-ion anisotropy and the interplane couplings stabilizing one 3 dimensional magnetic order of the Mn out of the four irreducible representations Γi.

modes. Moreover, due to the hexagonal symmetry, these three branches can be described by a single one, stemming from (0, 0, ) and going soft at (3, 0, ). The two other branches are identical but stem from (1, 0, ) and (2, 0, ) and go soft respectively at (4, 0, ) and (5, 0, ).

As shown in Fig. 2(c,d), spin wave calculations indicate that the former corresponds mostly to correlations between spin components along the c-axis (Fig. 2c), while the two latter correspond mostly to in-plane correlations (Fig. 2d). As a result, two acoustic-like and one opticallike spin wave emerge from each zone center. The gap ∆ of the latter occurs at around 40 cm -1 (5 meV) and is related to the easy-plane anisotropy by ∆ ≈ 3S √ D z J. As explained above, the degeneracy is lifted as soon as K 1 = K 2 . This becomes especially clear when looking at the dispersion along (1, 0, ). This point is illustrated in Fig. 2e-f, showing the calculated spin wave dispersions for the different Γ i magnetic configurations. We notice that the spectral weight of the lowest energy mode is strong at (1, 0, 1) for Γ 1,3 , and at (1, 0, 0) in the case of Γ 2,4 . This change has been put forward in Ref. 19 to explain that the spin reorientation process in HoMnO 3 involves a symmetry change from Γ 2,4 below T N to Γ 1,3 below T SR with a 180 • rotation of the spins in one of the triangular layer. This change was attributed to a change of sign of K 1 -K 2 induced by the change in position of the Mn ions from x M n = 0.325 above T SR to x M n ≈ 0.334 below. In Ref. 19, the authors assumed Γ 1 to be the ground state. Actually, non linear optical measurements as well as other probes suggest rather a Γ 3 magnetic order below T SR [26,27,29]. Since the Γ 1,3 and Γ 2,4 pairs are homometric, i.e. they cannot be distinguished by unpolarized neutron scattering if x M n is strictly equal to 1/3, we will assume that Γ 3 is the actual ground state in the following.

With these results in hand, we come back to the analysis of the experimental results. Fig. 2b displays the calculated inelastic neutron intensity for Q wave vectors along (h, 0, 1). Using J =2.5 meV and D z =0.65 meV, the high energy behavior of the spectrum is well reproduced. For instance, the spin wave labeled M originates from the optical-like mode at ∆. It is clearly visible on Figure 2. In Fig. 3 (panel A), it looks like a flat band along c * due to the weakness of the interplane couplings K 1 and K 2 . Finally, as shown in panel B, ∆ hardens towards 45 cm -1 when decreasing temperature. Note that some renormalisation of the Mn excitations has been reported in HoMnO 3 as due to magnon-phonon coupling [37].

The second dispersive spin wave labeled M' is the acoustic-like spin wave mode stemming from the zone center and consisting of two interlaced branches. It can be also identified in Fig. 3, thanks to the strong intensity expected for Γ 2,4 structures at (1, 0, 0). Below T SR , this spectral weight jumps to (1, 0, 1), as a result of the transition to Γ 1,3 . From panel B in Figure 3, the change of spectral weight from (1,0,0) to (1,0,1) at T SR can also be noticed. The bandwidth of M' is well captured for K 1 -K 2 = 0.003 meV. The gap of this mode originates from small single-ion anisotropies of the Mn 3+ (terms D y and D x ). An additional origin of anisotropy is expected to come from the Ho molecular field that should lead to a scaling of the gap with the Ho(4b) ordered moment. Indeed, at the mean field level, the coupling between the 3d and 4f magnetic species writes:

V = m,i,p, S m,i g m,i,p, J p,
where g is the coupling tensor, m, p and i, are indices that label the unit cell and the atom for the 3d and rare earth sublattices. J p, is the average magnetic moment at the rare earth site. Combined with the exchange Hamiltonian described in Eq. 1, V polarizes the Mn spins, hence plays the role of anisotropy. In principle, this should result in a gap in the spin wave spectrum as observed in YbMnO 3 [39]. However, in the present case, the expected scaling seems not to fully operate probably due to the coupling of the Mn spin waves with Ho CF excitations.

Indeed, in addition to the spin waves, non-dispersive excitations are also visible in Fig. 2 and3. These features correspond to CF transitions of the Ho 3+ ions on 4b and 2a sites and are labeled 0-6. As we shall see in the next part, those labels are also used to describe the THz data (see Figure 4). Note that the mode 5 is not visible in our neutron scattering experiments probably because it overlaps with the magnon M or because its neutron cross section is intrinsically weak. Panel B of Figure 3 shows the temperature dependence of CF transition modes 1 to 4. The intensity of those modes gradually increases with decreasing temperature while their positions slightly change, but the most remarkable feature is the formation of an additional line labeled 1 just below T SR . Note that the line labeled 0 appears at very low temperature only, hence probably due to the Ho ordering on site 2a. In addition, quite strikingly, the dispersion of M' eventually merges gradually below T SR into those CF lines (at least into those below 30 cm -1 ). Finally, the spectrum encompasses only flat bands below ≈ 10 K.

We point out that this very peculiar temperature evolution of M' shows that this branch crosses several Ho crystal field excitations, which makes very likely hybridization mechanisms (via a resonant process) between the Mn spin wave dispersion and the CF lines at low temperature.

SYNCHROTRON BASED THZ SPECTROSCOPY

To get a better insight into the spin dynamics, we have used THz spectroscopy as a complementary probe of both Mn spin waves and Ho Crystal field excitations. In this electromagnetic absorption technique, due to the high velocity of photons in the material, the measurements are considered in the long wave length limit and comparison with neutron data is therefore relevant for the (0, 0, 0) center of the Brillouin zone. There, all Ho crystal field excitations should be present. As for the the Mn spin waves, calculations similar to those presented in the previous section show that the four different spin configurations Γ i have the same signature at (0,0,0). It consists of a high energy branch (M) expected around 45 cm -1 which involves correlations between spin components perpendicular to the c-axis, and a single low energy branch M' corresponding to spin correlations along the c-axis that is pushed to finite energy when some spin anisotropy is present (see above and Fig. 2).

The THz spectrum of HoMnO 3 was first studied in [31] where the Mn (M) magnon mode as well as several Ho crystal field excitations were observed as a function of external magnetic field. The Mn magnon was more recently investigated under light irradiation [40] while its dependence on the Ho-Mn interaction is reported in [32]. Here, thanks to synchrotron based spectroscopy, we are ) Energy (cm -1 ) Energy (cm -1 ) able to measure with great precision the whole spectra from 0.3 to 1.8 THz (10 to 60 cm -1 ).

All the spectra recorded at different temperatures from 130 K or 100 K down to 6 K are shown in Fig. 4 for the three possible orientations (as regards the c-axis) of the THz electric e and magnetic h fields. Several absorption bands are observed from 11 cm -1 up to 55 cm -1 , quite differently from other hexagonal manganites such as YMnO 3 or ErMnO 3 , suggesting the influence of Ho ions in these THz spectra. Further inspection shows that, for the two different THz polarizations with the same magnetic field direction h⊥c, the spectra are very similar and differ substantially from the case where h c. This is a strong indication that all the observed excitations are of magnetic origin. The Mn magnon (M) expected around 45 cm -1 is clearly observed for h⊥c and seems to behave as expected for a conventional magnon by emerging below T N . It shifts to higher energies when the temperature is lowered (see Fig. 5). This magnon involves spin components perpendicular to this axis, in agreement with the THz magnetic selection rule h⊥c. The other Mn magnon (M') involving the spin component along the caxis is expected at lower energies for h c but is not easily distinguishable. All the other observed bands are specific to the Ho compound and should be related to Ho Closed and open symbols refer to h⊥c and h c respectively. Note the complexity of the spectra with several components that are all affected by TSR. crystal field excitations. Quite surprisingly, all spectra are strongly affected by the spin reorientation process at T SR .

A detailed description of the temperature dependence of the spectra are given in Fig. 5 Left: absolute and relative absorption measured at different temperatures. Right: results of gaussian fits for the absorption peaks with the position, area and width as a function of temperature. Two peaks are present. The lower energy peak appears smoothly when lowering the temperature and can be fitted below 50 K. The higher energy one, the stronger, is affected at TSR with a small shift in energy and increased width and area.

sorption peaks have been fitted with gaussian profile and their position, spectral weight (area) and lifetime (width) are reported. They all show, at the reorientation transition T SR , a slight change of their position in energy and more or less pronounced change of spectral weight.

In particular band 1 is observable only below T SR and for h⊥c (see Fig. 6). Band 2, always present for h⊥c, is observable only below T SR when h c (see Fig. 6). Band 3, always present for h⊥c, disappears below T SR when h c, while band 4 has a strong increase below T SR (see Fig. 7). Note that band 3 below T SR might have a shoulder and that band 4 encompasses certainly several modes. Finally, band 6 and maybe also band 5 have a side band that increases below T SR (see Figs. 8 and9). The Mn magnon (M) is also affected by the spin reorientation at T SR : an abrupt change, although weak, occurs in its position and line width while its area, which was increasing below T N , starts to decrease and becomes temperature independent below 10 K (see Fig. 5).

A comparison with the neutron data can be made keeping in mind that neutron scattering measures the imaginary part of the dynamical susceptibility that is equal to the THz absorption divided by the energy. Both techniques show remarkable agreement: At 20 K, below T SR , the observed five flat bands along (1, 0, )) in the neutron data (see Fig. 3) correspond nicely to band 1 -4 and to the magnon M observed in THz spectroscopy. The band 6 is also visible as a non dispersive branch at low temperature in the neutron data of Fig. 2. At 45 K, above T SR , five branches are visible in the neutron spectra, three of Left: absolute and relative absorption measured at different temperatures. Right: results of gaussian fits for the absorption peaks with the position, area and width as a function of temperature. Two peaks are present that are both affected at TSR. them correspond to bands 2-3-4 in THz, while the dispersive magnon M' is not distinguishable in THz, probably because it is below the energy window of the THz spectroscopy or mixed inside bands 3-4 (see Fig. 3). To go further into the analysis of the different magnon/crystal field contributions and their possible hybridization, we propose additional modeling of the complex Mn and Ho dynamics. We first focus on the crystal field excitations associated to the Ho ions.

CRYSTAL FIELD CALCULATIONS

Crystal field effects were calculated in a point charge model taking into account the whole oxygen coordination shell of the Ho 3+ ions (4f 10 ) for site 4b in C 3 symmetry and site 2a in C 3v symmetry [11,41,[START_REF] Hutchings | Solid State Physics: Advances in Research and Applications[END_REF]. The atomic coordinates were obtained from a refinement of neutron diffraction data at 10 K [19]. The associated Hamiltonian for each Ho 3+ ion in the ground multiplet J=8 can be expressed as:

H CF 4b = A 0 2 O 0 2 + A 0 4 O 0 4 + A 3 4 O 3 4 -(A 3 4 ) * O -3 4 + A 0 6 O 0 6 +A 3 6 O 3 6 -(A 3 6 ) * O -3 6 + A 6 6 O 6 6 + (A 6 6 ) * O -6 6 ( 4
)
H CF 2a = A 0 2 O 0 2 + A 0 4 O 0 4 + A 3 4 (O 3 4 -O -3 4 ) + A 0 6 O 0 6 +A 3 6 (O 3 6 -O -3 6 ) + A 6 6 (O 6 6 + O -6 6 ) (5)
for sites 4b and 2a, in terms of O m n Racah operators. A global screening factor of 20% was estimated in order to obtain the transition between the CF levels in the experimental energy range. The CF parameters resulting from this calculation, labeled set 1, are given in Table II. The deduced energy scheme and selection rules are given in Fig. 10. For Ho 4b, the energy scheme contains singlet (A) and degenerate doublets (E 1 E 2 ). Four different energy levels are calculated below 100 cm -1 . For Ho 2a, two singlets (A 1 , A 2 ) and one doublet (E) are present with four different energy levels in the same energy range. They are all allowed to occur for any orientation of the THz magnetic field. The calculated absorption is shown in Fig. 13(a) and (b) for both orientations of the THz magnetic field h at 40 K and 6 K respectively. The energy and relative intensities of the highest energy bands 5 and 6 are correctly reproduced. However the agreement between the calculations and the observed complex set of mixed excitations below 30 cm -1 is less accurate.

A 1 ↔#A 2# E↔#E# A 1,2 ↔#E# E↔#E# Ho#Site#4b#(C 3 #symmetry)# Ho#Site#2a#(C 3v #symmetry)# E 1 ⊕E 2 ↔!E 1 ⊕E 2 ! A↔#A# E 1 ⊕E 2 ↔#E 1 ⊕E 2 # A↔E 1 ⊕E 2 # 2# 2# 4# 4# A# A# E 1 ⊕E 2 # E 1 ⊕E 2 # 99#cm ;1# 0;0.6#cm ;1# E 1 ⊕E 2 # 3# 5# 6# 3# 5# 6# A 1# A 2# E# E# E# 5# 0;1.6#cm ;1# 15#cm ;1# 39#cm ;1# 49#cm ;1# 6# 2# 4# 15#cm ;1# 30#cm ;1# 3# h k c h ? c h k c h ? c
FIG. 10. Crystal field scheme and selection rules of the parameters set 1 for the two different Holmium sites with a tentative matching of the calculated CF transitions with the observed THz bands: Two transitions are expected for site 4b corresponding to band 2 and 3 and three for site 2a corresponding to band 4, 5 and 6. The selection rules and the symmetry of the levels are indicated below.

To go further, we took into account the coupling of the Ho 4b magnetic moments with the Mn ones expected to occur below T SR (discarding a possible magnetic coupling of the Ho 2a with Mn) [29]. This was achieved by including in the calculations the influence of the local magnetic field produced by the Mn moments on the Ho 4b site, similarly to the case of ErMnO 3 [11]. This field was chosen along the c-axis according to the observed cdirection of the Ho ordered magnetic moments. Its value was adjusted at 2 T so that the ordered Ho 4b magnetic moment extracted from the neutron diffraction data is reproduced at low temperature (see Appendix I). The influence of this magnetic field is to change slightly the energy of some peaks and to produce additional splitting of some CF levels (see Fig. 11b andc for the parameter set 1). We can notice in particular that the lowest energy CF peak corresponding to Ho 4b (in blue) in zero field splits into two peaks under magnetic field. However, experimentally, there are abrupt changes at T SR for all the THz bands, both on their position and their intensity, which cannot be accounted for by the sole effects of temperature and molecular field. This is a strong indication that electric effects with modification of electric charge screening are at play at T SR . This is further corroborated by the observation simultaneously of a sharp anomaly in the dielectric constant as well as a jump in the electric polarisation along the c-axis [30].

With little additional changes of some CF parameters (set 2 in Table II), we can reproduce qualitatively most of the observed features in the THz spectra (see Fig. 13d for T=6 K and a field of 2 T). In particular, bands 2, 3, 4, 5 and 6 are predicted at energies closer to the measured ones, although there are still some discrepancies with the calculated intensities. The presence of peak 3 above T SR for h c and its disappearance below T SR while a new peak labeled 2 appears at lower energy could be due to the field splitting ascribed to the Ho-Mn coupling. Band 4 for h c below T SR could then result from a mixture of several peaks including the higher energy mode from the split peak (see horizontal bracket in Fig. 13d).

Strikingly, for h ⊥ c, band 1 appearing at T SR in the measurements is missing in the calculations and therefore cannot be accounted for by the effect of a magnetic field. This discrepancy between calculations and experiment is another indication that hybridization is at play between Ho CF excitations and the lowest energy Mn magnon (M').

DISCUSSION

Our analysis and calculations for the Mn magnetic order and associated spin waves, using the spin Hamiltonian of eq.1, can account for most of the experimental data. Mainly, two spin branches M and M' are observed with INS that are correctly calculated (Fig. 2 and3). In particular, the spin reorientation at T SR is clearly seen as a change in the spin wave dispersion along c * for the Energy"(cm -1 )"

Energy"(cm -1 )"

Absorp5on"(arb."units)" ) µ 0 H=0" T=40"K" Site"4b" Site"2a" Parameters"1" µ 0 H=0" T=6"K" Parameters"1" µ 0 H(4b)=2"T"" T=6"K" Parameters"1"

µ 0 H(4b)=2"T" T=6"K" Parameters"2" { (a)" (b)" (c)" (d)"
FIG. 11. Calculated CF THz absorption spectra for both orientations of the THz magnetic field h and for Ho sites 4b (blue) and 2a (yellow), using parameters set 1 at 40 K and zero magnetic field (a), at 6 K and zero magnetic field (b) and at 6 K with a 2 T magnetic field parallel to the c-axis. The dashed line and the two grey arrows for h c indicate the splitting of the lowest energy Ho 4b CF excitation under magnetic field. Panel (d) corresponds to the calculations performed using the parameter set 2 at 6 K and under a magnetic field of 2 T along the c-axis with a tentative matching of the calculations with the measured bands of excitations.

branch M' and is correctly captured when the change of Mn spin configuration from Γ 4 to Γ 3 is introduced. The second branch M is also observed with THz spectroscopy (Fig. 5), in agreement with the model. As for the non dispersive Ho CF excitations observed with INS and THz spectroscopy, our calculations using the CF hamiltonian described by eq.4 and 5 is able to reproduce the experimental data qualitatively when two ingredients are included: first, a change in charge screening at T SR that explains the observed shifts and changes in intensity of the whole CF spectrum, and second, a molecular field along the c-axis produced by the Mn order. Similarly, a molecular field produced by the Ho 4b order may be introduced to account for the increased anisotropy gap in the Mn spin wave branch M'. However our calculations do not agree quantitatively with the whole experimental data: they cannot describe properly how the Mn spin branch M' become non dispersive along (1, 0, ) and seems to merge into Ho CF excitations below 20 K. Moreover, they cannot explain the THz band 1 that appears at T SR and is locked at 3 cm -1 below band 2 (Fig. 6). Ad-ditional ingredients should be at play there, that couple the 3d and 4f magnetic species dynamically.

It is possible for Mn Magnons and Ho crystal field excitations to hybridize, if they are close in energy, which is the case for the Ho manganite where 6 different Ho crystal field excitations lie in the energy range of the Mn spin waves (see Fig. ??). They also should not be forbidden by symmetry. This can be checked using a group representation analysis for both the Mn spin waves and the Ho crystal field excitations : if they contain at least one common irreducible representation, the coupling is allowed. The details are given in appendix II. According to [43], for the Mn magnetic orders with zero propagation vector corresponding to the four different one dimensionnal Γ i , the spin wave representation may be decomposed into τ = A 2 B 1 E 1 E 2 where we use the irreducible representations of the C 6v (6mm) point group. For Ho 4b and 2a crystal field excitations, their decompositions are already given in Fig. 10, in the local point groups associated to the two different sites: C 3 (3) for Ho 4b, and C 3v (3m) for Ho 2a. Their energy scheme contains singlets (A for Ho 4b, A 1 , A 2 for Ho 2a) and doublets (E 1 , E 2 for Ho 4b, E for Ho 2a). To infer compatibility with those of the Mn magnons, we must consider their decomposition into the irreducible representations of the C 6v point group. For Ho 4b, singlets A are decomposed into {A 1 , A 2 , B 1 , B 2 } and degenerate doublets (E 1 E 2 are decomposed into {E 1 , E 2 }. As for the Ho 2a, singlets A 1 and A 2 are decomposed into {A 1 , B 2 } and {A 2 , B 1 } respectively, and doublets E into {E 1 , E 2 } as for Ho 4b. It can be easily seen that they all have in common at least one irreducible representation contained in the Mn spin wave representation except for Ho 2a A 1 singlet that is calculated around 40-45 cm -1 and corresponds to band 5.

We conclude then that, in the energy range relevant for the experimental data reported in this study, all the CF excitations, except one (band 5 associated to Ho 2a), are allowed to couple to Mn spin waves. Looking back at neutron experiments and the spin wave calculations reproduced in Fig. 12, it can be seen that, when the temperature is lowered close to T SR , the magnon M' with small gaps at the zone centers (1,0,0) and (1,0,1) become close in energy to Ho (4b) CF excitations (bands 2 and 3-4), with different spectral weight for the magnons depending on the magnetic order Γ 1,3 or Γ 2,4 . A coupling mechanism may be at play here that contributes to the spin reorientation process and the lost of dispersion for M' along c * .

It is worth noting that Ho and Mn ordered magnetic moments below T SR are orthogonal. This suggests that, in the Mn -Ho coupling Hamiltonian of eq., the coupling tensor g shows non diagonal terms, as it is the case for instance for the pseudo dipolar, the Dzyaloshinkii-Moriya interaction [39] or the highest order trigonal anisotropy [44]. Furthermore, the 2 T effective magnetic field expe- rienced by the 4b Ho suggests a coupling strength g ∼ 0.1 K. In this context and given the neutron scattering data and the CF analysis detailed above, we may to consider a possible dynamical coupling between i) in plane Mn spin wave modes with Ho CF transitions polarized along the c-axis and ii) out of plane Mn spin wave modes with Ho CF transitions polarized in the (a,b) plane. In this view, it is tempting to assign band 1 to the Mn spin wave hybridized with band 2 (and also likely with band 3). In this scenario, the in plane character of the CF transitions would be transferred to the hybrid mode which would now appear in the h⊥c channel as observed experimentally in the THz measurements. This hybridisation occurs only below T SR when the spectral weight of these excitations is enhanced and their energy matches. The energy difference between band 1 and 2, around 3 cm -1 experimentally, is in agreement with the coupling strength g when all the Ho-Mn links are considered (12 of them).

CONCLUSION

To summarize, we have reexamined the THz dynamics present in the multiferroic compound HoMnO 3 . The peculiarity of this Ho compound in the hexagonal manganites is the occurrence of a spin reorientation below T SR =37 K. Remarkably, several Ho CF excitations are present in the same energy range as the Mn spin waves. Thanks to the complementary neutron and electromagnetic wave probes, we have identified five Ho 4b and Ho 2a CF excitations and shown how Mn spin waves as well as Ho CF are affected by the reorientation process at T SR . Using spin wave and CF calculations, we are able to reproduce qualitatively most of the observed features. At T SR , the Mn spin configuration switches from the high temperature Γ 4 configuration to Γ 3 . A change in the spectral weight is observed along the (1, 0, ) direction that is correctly calculated in a simple hamiltonian for the Mn win waves. At the same time, changes occur in the CF excitations that can be accounted for when coupling to the lattice and surrounding oxygen charges are introduced. Then, Ho CF excitations start to match in energy the Mn branch at different Brillouin zone centers. Further inspection of the neutron and THz data points towards additional Ho-Mn coupling mechanisms that show up below T SR in two different ways : first as an additional excitation observed in THz spectroscopy that is locked 3 cm -1 below a CF excitation and secondly the progressive disappearance of the Mn dispersion along c * . Our detailed study on HoMnO 3 gives finally new insight on the very rich rare earth -transition metal dynamical couplings that show up in the THz range.

This work was financially supported by the ANR-13-BS04-0013-01.

APPENDIX I

In the Ho CF calculations, we introduce a constant magnetic field to simulate the effect of the Ho-Mn coupling. We used the measured Ho 4b and Mn ordered moment extracted from neutron diffraction data [39]: Ho 4b moments are observed along the c-axis; a molecular field of 2 T along this direction is able to polarize Ho spins with the correct magnitude and temperature dependence at least below 10 K (see Fig. 13). The Mn moment varies only slightly below T SR thus validating the assumption of a constant molecular field.

APPENDIX II

For Mn Magnons and Ho crystal field excitations to hybridize, there may be symmetry requirements that can be derived through a group representation analysis for the Mn spin waves and the Ho crystal field excitations. The Mn order with zero propagation vector is described by the isomorphic space group C 3 6v . The Mn spin wave representations can be reduced into representations of the associated point group C 6v (6mm). This group contains the lower symmetric C 3 (3) and C 3v (3m)that are best describing the Ho CF shemes. In the present analysis, to look at possible coupling mechanisms with the Mn spin waves, we will rewrite the CF representations in the more symmetric C 6v . Its character table is given in table III. The symmetry of the different Mn spin waves have been reported in [43] with Kovalev notation. We will rather use Mulliken notation. To simplify calculations, the spin wave analysis has been performed in a more symmetric space group than C 3 6v , namely D 4 6h , which means neglecting some oxygen displacements and therefore some weak magnetic anisotropies. From [43], it follows that the spin wave representation can be decomposed into A 2g B 1g E 1g E 2g . Going back to the lower symmetric C 3 6v space group, gerade und ungerade modes are no longer relevant so that the spin wave decomposition becomes A 2 B 1 E 1 E 2 . When spin anisotropies are introduced, the degeneracy of the low energy branches is lifted as can be seen in our spin wave calculations (Fig. 2), dictating the use of magnetic Shubnikov groups.

For Ho CF schemes, The character tables for the associated point groups C 3 (3) and C 3v (3m) are reported in table IV and V. Their correlations to the more symmetric C 6v group that we will use are given in table VI. For Ho 4b, in its local C 3 symmetry, two different kind of symmetry are present for the crystal field states: singlets are described by the irreducible representation A and degenerate doublets described by E 1 E 2 . These irreducible representations induce, in the higher symmetric C 6v group, states that are decomposed into linear combinations of {A 1 , A 2 , B 1 , B 2 } for the singlets and {E 1 , E 2 } for the degenerate doublets. For the Ho 2a, in the local C 3v symmetry, three different symmetries are present: singlets A 1 and A 2 , and doublets E. They induce, in C 6v symmetry, states that are decomposed into {A 1 , B 2 } and {A 2 , B 1 } respectively for the singlets and {E 1 , E 2 } for the doublets.

From these considerations, we can now look at possible coupling between the Mn spin waves and the Ho 4b and 2a crystal field excitations that are allowed by symmetry, that is if their representations contain at least one common irreducible representation. This is the case fot all CF excitations except for the Ho 2a A 1 singlets.

C 6V
6mm 1 2 3 6 m d m v functions Multiplicity -1 1 2 2 3 3

.

A 1 Γ 1 1 1 1 1 1 1 z, x 2 + y 2 , z 2 A 2 Γ 2 1 1 1 1 -1 -1 J z B 1 Γ 3 1 -1 1 -1 1 -1 . B 2 Γ 4 1 -1 1 -1 -1 1 . E 2
Γ 5 2 2 -1 -1 0 0 z, x 2y 2 , xy E 1 Γ 6 2 -2 -1 1 0 0 x, y, xz, yz J x , J y All the tables reproduced here where established thanks to [45].

C 6v C 3 C 3v 6mm 3 3m A 1 A A 1 A 2 A A 2 B 1 A A 2 B 2 A A 1 E 2 E 1 , E 2 E E 1 E 1 , E 2 E

FIG. 1 .

 1 FIG. 1. Sketch of the HoMnO3 crystallographic and magnetic structures. (a) Projection of the different atomic positions onto the ab-plane, omitting the oxygens for clarity. The Mn ions on the 6c Wyckoff site have a x coordinate within the ab-plane close to 1/3 and a z coordinate along the c-axis equal to 0 (red spheres) or 1/2 (black spheres). The rareearths occupy two different Wyckoff sites (2a and 4b). The main intra-plane (J) and inter-layer (K1,2) exchange interactions between the Mn magnetic moments are indicated. (b) Magnetic arrangement of the Mn ions with zero propagation vector corresponding to the four one-dimensional irreducible representations Γi of the P63cm space groupe. The homometric pairs, indistinguishable through neutron scattering if xMn=1/3, are indicated by the boxes.

FIG. 2 .

 2 FIG. 2. Spin dynamics in HoMnO3 (a) Mn 3+ spin waves and Ho 3+ CF excitations measured by INS along the (h, 0, 1) reciprocal space direction. These data have been taken on IN8 at 2K with a medium energy resolution of about 0.8 meV. M and M' denote the Mn dispersive spin waves that can be seen at the zone center. Labels 1-6 refer to non dispersive Ho CF lines. (b) Comparison with spin wave calculations using the model described in the text for the Γ1,3 configurations of the Mn magnetic structure. (c,d) Spin wave calculations along (h, 0, 0) in the Γ1,3 configurations for spin fluctuations respectively along the c-axis and within the basal hexagonal plane. (e-h) Spin wave calculations along (1, 0, ) for the 4 different possible Γi. In the calculations of panels (b-f), modest but finite values of Dx and Dy were taken into account to ensure the stability of the various Γi, see TableI. This is responsible for the rise of the small gaps at zone centers. All the spin wave calculations have been performed using the SpinWave software[38].

FIG. 3 .

 3 FIG. 3. Panel A: Inelastic neutron scattering color maps in HoMnO3 for the (1, 0, ) reciprocal space direction at different temperatures, ranging from 65 K down to 1.5 K. The spin reorientation occurs at TSR =37 K. Dispersing spin waves branches (labeled M and M')a ssociated to Mn order can be observed at 65, 45 and 27 K. The spectral weight is typical of the Γ2,4 configuration at 65 K and 45 K, but it becomes typical of Γ1,3 at 27 K. In addition, flat modes (0-4), likely Ho crystal field excitations, are also observed. For the temperatures below 27 K, the spin waves gradually fade away into the CF excitations and only non-dispersive excitations remain at 1.5 K (below THo). The dots added on the color maps were obtained by gaussian fits of the individual constant-Q E-scans. Panel B : inelastic neutron scattering color maps in HoMnO3 for the (100) and (101) reciprocal space positions as a function of temperature. The dots added on the color maps were obtained by gaussian fits of the individual constant-Q E-scans. The dashed vertical lines indicate the transition temperatures.

5 FIG. 4 .FIG. 5 .

 545 FIG.4. THz absorption spectra of HoMnO3 for the 3 different orientations of the electromagnetic wave e, h fields with respect to the crystal c-axis, using two different samples. Upper panels: absolute absorptions. Lower panels: relative absorptions with a reference above the magnetic ordering temperature at 100 K or 130 K.

FIG. 6 . 4 FIG. 7 .

 647 FIG.6. Details of the THz spectra around band 1 and 2. Left: relative absorption measured for two different orientations of e and h on the same sample. Right: results of the gaussian fits for the absorption peaks with the position, area and width as a function of temperature. Closed and open symbols refer to h⊥c and h c respectively. Below TSR, two new excitations appear abruptly: band 1 for h⊥c and band 2 for h c.

5 FIG. 8 .

 58 FIG.8. Details of the THz spectra around band 5 when h c. Left: absolute and relative absorption measured at different temperatures. Right: results of gaussian fits for the absorption peaks with the position, area and width as a function of temperature. Two peaks are present. The lower energy peak appears smoothly when lowering the temperature and can be fitted below 50 K. The higher energy one, the stronger, is affected at TSR with a small shift in energy and increased width and area.

6 FIG. 9 .

 69 FIG.9. Details of the THz spectra around band 6 when h⊥c. Left: absolute and relative absorption measured at different temperatures. Right: results of gaussian fits for the absorption peaks with the position, area and width as a function of temperature. Two peaks are present that are both affected at TSR.

FIG. 12

 12 FIG.12. (a) Mn spin wave calculations along (h, 0, 0) for spin fluctuations within the (a, b) plane (upper panel) and along the c-axis (lower panel). The calculated CF levels are added for the Ho 4b site (red and blue lines) as well as the Ho 2a site (yellow and green lines). The coupling of these CF excitations with the Mn magnons is always allowed by symmetry for all but for the dashed one.

FIG. 13 .

 13 FIG.13. Temperature dependence of the ordered magnetic moments Ho 4b extracted from the refinement of the neutron diffraction data (black dots) compared to the calculated magnetic moment (blue line) induced by a contant magnetic field of 2 T due to the coupling with the Mn magnetism. The red stars show the Mn ordered moment also extracted from the neutron diffraction data.

TABLE III .

 III C6v (6mm) character table

TABLE IV .

 IV C3 (3) character table + y 2 , J z E 2 Γ 2 1 ω 2 ω x, y, xz, yz, x 2y 2 , xy,J x ,J y E 1 Γ 3 1 ω ω 2

	C 3	3 1 3 + 3 -	functions
	A Γ 1 1	1	1	z, x 2

TABLE V .

 V C3v (3) character table , y, xz, yz, x 2y 2 , xy,J x ,J y

	C 3v	3m 1 3 m	functions
	Multiplicity -1 2 3
	A 1	Γ 1 1 1 1	z, x 2 + y 2 , z 2
	A 2	Γ 2 1 1 -1	J z
	E	Γ 3 2 -1 0 x

TABLE VI .

 VI Correlations between the global frame point group C6v (6mm) and the two local frame point groups C3 (3) and C3v (3m)


Anthony Carpentier

Nicolas Galopin

Olivier Chadebec

Gérard Meunier

Modeling of Magneto-Mechanical Coupling Using Magnetic Volume Integral and Mechanical Finite-Element Methods

published or not. The documents may come   L'archive ouverte pluridisciplinaire

HE modeling of electromagnetic devices using deformable materials involves the resolution of the mechanical and electromagnetic equations including coupling terms in each physic. A classical method to solve numerically such problems is the finite-element method [START_REF] Galopin | Finite element modeling of magnetoelectric sensors[END_REF]. However, for devices with a huge volume of free space compared with the active structure or high size ratio between the geometric objects like micro-electro-mechanical systems (MEMS) [START_REF] Cugat | Magnetic micro-actuators and systems (MAGMAS)[END_REF], problems of accuracy could be present [START_REF] Rakotoarison | Formal sensitivity computation of magnetic moment method[END_REF] for the resolution of the magnetostatic problem. To avoid this limitation, this paper investigates the magneto-mechanical coupling of devices composed of ferromagnetic and deformable materials in which the magnetic problem is solved with an integral formulation while the mechanical problem is solved by the finite-element method.

Considering only static case and no magnetostrictive effects, the equations of the magneto-mechanical problem can be solved separately. The magnetostatic problem is an open boundary problem for the magnetic field whereas the mechanical part is restricted to the solid domain. In this case, the integral formulations of the magnetostatic field problems are particularly advantageous for the numerical solution of open boundary problems which include magnetic materials. Indeed only the active regions containing these materials need to be discretized [START_REF] Trowbridge | Integral equations in electromagnetics[END_REF]. Therefore, a volume integral method is used to solve the magnetostatic problem.

With the assumptions of no magnetostrictive phenomena, only the magnetic forces are involved in the magnetomechanical coupling. A local application of the virtual work principle has been proposed in [START_REF] Carpentier | Application of the virtual work principle to compute magnetic forces with a volume integral method[END_REF], to compute the nodal integral method. This method is used here in order to compute a source term of a finite-element description of the mechanical problem. Advantage of this approach is to keep the same mesh for both resolutions and to consider only the active regions.

This paper presents in the first part a description of the problem with the associated equations. The second part is dedicated to the magneto-mechanical formulation in which a volume integral method is used to solve the magnetostatic problem. A local application of the virtual work principle, adapted to the framework of the integral methods, is proposed to compute the local magnetic forces. The finite-element method to solve the mechanical problem is also introduced. Finally, an application is given in the last part to illustrate the potential of the developed method.

II. DESCRIPTION OF THE PROBLEM

Let us consider the magneto-mechanical problem defined by the 2-D regions, f and m , filled, respectively, with ferromagnetic and nonmagnetic materials (Fig. 1). Each of them has an isotropic linear elastic mechanical behavior. Both regions do not overlap and their union is noted = f ∪ m . Let f , m , and be, respectively, the boundaries of the regions f , m , and . The regions f and m are in contact over the boundary c = f ∩ m . The region 0 is associated to the free space.

Modeling of Magneto-Mechanical Coupling Using Magnetic Volume Integral and Mechanical Finite-Element Methods

Anthony Carpentier 1 , Nicolas Galopin 1,2 , Olivier Chadebec 1,3,4 , and Gérard Meunier 1,3 The magneto-mechanical coupling with deformable bodies is studied using different numerical methods for both physics. The magnetic behavior is considered linear and without magnetostrictive effects; then equations of the magneto-mechanical coupling can be solved separately. A volume integral method using the magnetic scalar potential formulation is used. This method is particularly well suited to compute fields without meshing the air domain. A local application of the virtual work principle, adapted to the integral formulation, is used to compute the nodal magnetic force distribution at the origin of the magneto-mechanical coupling with deformable bodies. An advantage of this approach is that the same mesh can be used by both methods during the resolution steps.

Index Terms-Finite-element method, integral method, magneto-mechanical coupling, nodal magnetic forces.

I. INTRODUCTION

T magnetic force distribution in the framework of the volume

Considering only the static case, the magnetic behavior is governed by the Maxwell's equations

∇×H = J (1) ∇ • B = 0.
(

The magnetic soft behavior is considered linear. Therefore, in an isotropic homogeneous ferromagnetic material without magnetostrictive phenomena, the magnetic behavior law is defined by

B = μ 0 (H + M) (3)
where μ 0 is the magnetic permeability of free space and M the magnetization. The fields M and H are linked by

M = χH ( 4
)
where χ is the magnetic susceptibility.

The mechanical behavior is governed by the application of local equilibrium equations for a static case

∇ • [σ ] + f = 0 ( 5 )
where [σ ] is the stress tensor and f the volume force density.

Only the magnetic interaction is considered (all other potential energy fields are neglected as well as the frictional forces), thus the volume force density is associated to the magnetic forces acting on the ferromagnetic material f . The mechanical behavior law is defined by

[σ ] = C [s] ( 6 )
where [s] is the strain tensor and C the stiffness tensor. The materials are considered homogeneous and isotropic with a elastic mechanical behavior. In this case, the stiffness tensor C can be expressed in function of Young's modulus E and Poisson's ratio ν

C i j kl = E 1 + ν ν 1 -2ν δ kl δ i j + 1 2 δ ik δ jl + δ il δ j k ( 7
)
where δ is the Kronecker's symbol. With the assumption of small deformations, the strain tensor is linked to the displacement field u (specifying the position of a point from a reference frame) by the linearized equation

[s (u)] = 1 2 ∇u + ∇u T . ( 8
)

III. MAGNETO-MECHANICAL FORMULATION

A. Magnetostatic Formulation

The magnetic field is written as the sum of the reduced magnetic field created by the magnetic material H red and the magnetic source field created by currents flows H 0

H = H red + H 0 . ( 9
)
The Maxwell's equations ( 1) and ( 2) are satisfied when the reduced magnetic field H red is expressed at any point of coordinates r by [START_REF] Durand | Magnétostatique[END_REF] H

red (M, r) = - 1 4π ∇ a M r • r -r |r -r | 3 d . ( 10
)
The magnetic field H [START_REF] Henrotte | The eggsheel approach for the computation of the electromagnetic forces in 2D and 3D[END_REF], solution of the magnetostatic problem, verifies the (10) and ( 4) at any point of the ferromagnetic region f H (r) = H red (χH, r) + H 0 (r) . [START_REF] Henrotte | The eggsheel approach for the computation of the electromagnetic forces in 2D and 3D[END_REF] The simply connected region f containing no current sources, the magnetic fields H and H 0 derive, respectively, from the scalar potentials and 0 . The integration of (11) leads to the following volume integral equation [START_REF] Han | Integral equation method using total scalar potential for the simulation of linear or nonlinear 3D magnetostatic field with open boundary[END_REF], [START_REF] Kettunen | Volume integral equations in non-linear 3-D magnetostatics[END_REF]:

(r) + 1 4π f χ ∇ r • r -r |r -r | 3 d = 0 . ( 12
)
The region f is meshed with tetrahedra elements using first order shape functions

= i α i i ( 13
)
where α i and i are, respectively, the shape function and the degree of freedom of the potential associated to the ith node of the mesh. A collocation method at mesh nodes is used to solve [START_REF] Vandevelde | A survey of magnetic force distributions based on different magnetization models and on the virtual work principle[END_REF] and leads to solve the following system of algebraic equations

([I ] + [A (χ)]) = 0 ( 14
)
with

I i j = δ i j (15) A (χ) i j = 1 4π f χ ∇α j r • r i -r |r i -r | 3 d ( 16
)
where r i is the coordinates of the ith node. There are no direct analytical expressions for the source potential 0 . A numerical approach presented in [START_REF] Carpentier | Resolution of nonlinear magnetostatic problems with a volume integral method using the magnetic scalar potential[END_REF] is then used. The volume integral method allows to mesh only the active materials. The main drawback is that the interaction matrix [A] is full. The compression methods can overcome this limitation [START_REF] Carrier | A fast adaptive multipole algorithm for particle simulations[END_REF].

B. Local Magnetic Force Computation

Several methods are available in the literature to compute the local magnetic forces: Maxwell's tensor, virtual works, equivalent source or charge methods, eggshell method [START_REF] Henrotte | The eggsheel approach for the computation of the electromagnetic forces in 2D and 3D[END_REF] etc. These methods are several times compared each other [START_REF] Vandevelde | A survey of magnetic force distributions based on different magnetization models and on the virtual work principle[END_REF]. The virtual work principle is more general [START_REF] Bossavit | Virtual power principle and Maxwell's tensor: Which comes first?[END_REF] and has already been successfully applied to compute magnetic forces inside electro-magnetic MEMS [START_REF] Galopin | Finite element modeling of magnetoelectric sensors[END_REF]. Considering that, an adaptation of the virtual work principle in the framework of the volume integral method is used [START_REF] Carpentier | Application of the virtual work principle to compute magnetic forces with a volume integral method[END_REF].

A virtual displacement u n is applied on the mesh node n. According to the virtual work principle [START_REF] Ren | Local force computation in deformable bodies using edge elements[END_REF], the component along u n of the nodal magnetic force F n u is deduced from

F n u = dω co mag du n I ( 17
)
where I is the electric current which is kept constant and ω co mag the magnetic co-energy. This last could be written [START_REF] Stratton | Electromagnetic Theory[END_REF]

ω co mag = ω co ferro + μ 0 2 ∪ 0 H 0 • H 0 d ( 18
)
with

ω co ferro = μ -μ 0 2 f H • H 0 d . ( 19
)
The expression ( 18) is used in (17) to obtain

F n u = dω co ferro du n I . ( 20
)
The relation ( 20) is well suited for the integral formulations, due to the fact that the domain of the integral ( 19) is limited to the ferromagnetic region.

Using the change of basis between the physical mesh element e and the reference element ˆ described by the Jacobian matrix of the transformation [ J e ], the expression (20) can be developed into

F n u = e μ -μ 0 2 ˆ H • dH 0 du n + dH du n • H 0 |J e | + (H • H 0 ) d |J e | du n d . ( 21
)
The computation of each derivative term has been discussed in [START_REF] Carpentier | Application of the virtual work principle to compute magnetic forces with a volume integral method[END_REF].

Each component of the nodal magnetic force F n is obtained taking the virtual displacement u n , respectively, as the basis vectors (e 1 , e 2 , and e 3 ). Therefore, using the relation (21) the nodal magnetic force F n is given by

F n = 3 k=1 F n e k e k . ( 22
)
This adaptation of the virtual work principle allows computing local magnetic forces without computation in the free space. Its main drawback is the resolution of two systems of algebraic equations for each mesh node and each space dimension.

C. Mechanical Formulation

Let us consider the region is submitted to mechanical displacement on boundary u and to applied mechanical stresses on boundary σ . The corresponding boundary conditions are defined by

[σ ] • n = f σ on σ (23) u = u 0 on u (24)
where σ and u are the parts of the boundary (Fig. 1). The weak formulation associated to the mechanical problem (5) is defined by

w∇ • [σ ] d = -w • f d ∀w ( 25
)
where w is a test function verifying the boundary condition (23). Using the rule product for the divergence, the Stoke's theorem and the boundary conditions ( 23) and ( 24), the weak formulation (25) can be written The region m is meshed with tetrahedra elements using first order shape functions

[s (w)] [σ ] d = w • f d + σ w • f σ d , ∀w. (26)
u = i α i U i ( 27
)
where α i and U i are, respectively, the shape function and the vector degree of freedom of the displacement u associated to the ith node of the mesh. The finite-element method applied to the weak mechanical formulation (26) leads to the resolution of the algebraic equation system

[K ] U = F + Q (28)

with

K i j = s λ j C [s (λ i )] d ( 29
)
F i = λ i • f d (30) Q i = λ i • f σ d (31)
where λ i is the vector shape function associated to the vector degree of freedom i . The vector F is equal to the vector of the nodal magnetic force (22).

IV. APPLICATION

Let us study, a MEMS actuator (Fig. 2) [START_REF] Schneider Electric | Microsystème à commande électromagnétique[END_REF] in which an iron plate is fixed to a beam which is embedded at its extremities [u 0 = 0 for the boundary condition (24)]. No external force is applied to the structure [f σ = 0 for the boundary condition (23)]. The description of the MEMS and the corresponding mesh are presented on the Fig. 3. The magnetic behavior of the MEMS is considered linear with an isotropic homogeneous ferromagnetic material of magnetic permeability 1000 (H/m). Its mechanical behavior is considered elastic with a homogeneous and isotropic material of Young's modulus E = 110 GPa and Poisson's ratio ν = 0.3. The MEMS is actuated by a magnetic field produced by a source. Without magnetostrictive effects, the magnetic forces are the origin of the magneto-mechanical coupling and generate torsional strain of the beam which rotates around its axis.

The previous volume integral method [START_REF] Ren | Local force computation in deformable bodies using edge elements[END_REF] is used to solve the magnetostatic problem. The computed magnetic field is given in Fig. 4. The virtual work method ( 22) is applied to compute the magnetic nodal force distribution which is presented on Fig. 5. These nodal forces become the source term of the mechanical problem which is solved using the finite-element method (28). Fig. 6 represents the computed displacement with the new position of the MEMS and the reference configuration. This is in agreement with the torsion motion of the beam which rotates around its axis.

V. CONCLUSION

A magneto-mechanical coupling problem is solved in the framework of the volume integral and finite-element methods. The resolution of the magnetic behavior with the volume integral method allows to compute efficiently the nodal magnetic force distribution by an adapted virtual work principle formulation. This force distribution is then a source term for the mechanical problem which is achieved by finite-element method. Associated to an accurate magnetic force computation, the advantages of this approach is to consider only the active regions during the magneto-mechanical resolution, and to use the same mesh during the resolution steps. A work in progress is the investigation of the nonlinear case.

Fig. 1 .

 1 Fig. 1. Description of the magneto-mechanical problem.

Fig. 2 .

 2 Fig. 2. Description of the MEMS.

Fig. 3 .

 3 Fig. 3. Model of the MEMS with the boundary conditions.

Fig. 4 .

 4 Fig. 4. Magnetic field distribution.

Fig. 5 .

 5 Fig. 5. Nodal magnetic forces distribution.

Fig. 6 .

 6 Fig. 6. Displacement distribution with the new position of the MEMS.