Unusual magnetoelectric effect in paramagnetic rare-earth langasite

Violation of time reversal and spatial inversion symmetries has profound consequences for elementary particles and cosmology. Spontaneous breaking of these symmetries at phase transitions gives rise to unconventional physical phenomena in condensed matter systems, such as ferroelectricity induced by magnetic spirals, electromagnons, non-reciprocal propagation of light and spin waves, and the linear magnetoelectric (ME) effect—the electric polarization proportional to the applied magnetic field and the magnetization induced by the electric field. Here, we report the experimental study of the holmium-doped langasite, HoxLa3−xGa5SiO14, showing a puzzling combination of linear and highly non-linear ME responses in the disordered paramagnetic state: its electric polarization grows linearly with the magnetic field but oscillates many times upon rotation of the magnetic field vector. We propose a simple phenomenological Hamiltonian describing this unusual behavior and derive it microscopically using the coupling of magnetic multipoles of the rare-earth ions to the electric field.


INTRODUCTION
Novel materials and physical effects are essential for the continuing progress in nanoelectronics 1 . The coupling between electric and magnetic dipoles in magnetoelectric (ME) and multiferroic materials enables electric control of magnetism that can significantly reduce power consumption of magnetic memory and data processing devices [2][3][4][5][6] .
The linear ME effect [7][8][9] occurs in centrosymmetric crystals, such as Cr 2 O 3 10 and LiCoPO 4 11 , in which both time reversal and inversion symmetries are spontaneously broken by an antiferromagnetic spin ordering inducing toroidal, monopole, and other unconventional multipole moments 12 . This effect allows for control of antiferromagnetic domains and gives rise to unidirectional light transmission at THz frequencies associated with electromagnons [13][14][15] .
However, the linear ME effect in centrosymmetric crystals requires a perfect match between magnetic and lattice structures, which narrows the pool of ME materials substantially. Symmetry requirements are less stringent for magnets with noncentrosymmetric crystal lattices, such as boracites 16 , copper metaborates 17 , and the skyrmion material, Cu 2 OSeO 3 18 , in which diverse magnetic orders exhibit ME behavior.
A bi-linear ME effect-electric polarization proportional to the second power of the applied magnetic field-does not even require time-reversal symmetry breaking and can occur in the paramagnetic state of non-centrosymmetric magnets through a variety of miscroscopic mechanisms 19,20 . This higher-order effect is by no means weak: the magnetically induced electric polarization in the recently studied holmium hexaborate is much higher than the polarization of linear MEs at comparable magnetic fields 21 .
The Ho-doped langasite studied in this paper does not seem to fit into this classification of ME effects. Langasite, La 3 Ga 5 SiO 14 , and related compounds have been synthesized as early as the 1980s 22,23 and actively studied since then 24 for their interesting non-linear optical, elastic and piezoelectric 25,26 properties and applications in acousto-and electro-optics 27,28 . The noncentrosymmetric crystallographic space group P321 of langasites allows for optical rotation and piezoelectricity, but not for spontaneous electric polarization due to the presence of orthogonal threefold and twofold symmetry axes. The study of ME properties was focused on Fe-substituted langasites showing non-collinear periodically modulated orders of Fe spins 29 . The interplay between magnetic frustration, structural and spin chiralities gives rise to a rich variety of multiferroic behaviors 30,31 .
The Ho-doped langasite does not order down to lowest temperatures. Yet, its magnetically induced electric polarization oscillates 6 times when the direction of the magnetic field rotates through 360 ∘ in the ab-plane, which is indicative of a highly nonlinear ME response proportional to (at least) the sixth power of the magnetic field. On the other hand, the magnitude of the induced polarization grows linearly with the applied field, which seems to be at odds with the absence of any magnetic ordering in this material. This ME behavior is very different from that of classical MEs, such as Cr 2 O 3 . We show that it can be understood on the basis of langasite symmetry and a hierarchy of energy scales in the spectrum of the magnetic Ho ion.

Experiment
The samples used in this work are diluted rare-earth langasites, Ho x La 3−x Ga 5 SiO 14 , with x = 0.043 ± 0.005. The crystal structure of R 3 Ga 5 SiO 14 langasites 32,33 is shown in Fig. 1a. The rare-earth sites, R, have a rather low C 2 symmetry with the rotation axis along the a-axis of the crystal, allowing for electric dipole moments on R-sites, which in the case of magnetic Ho ions can depend on an applied magnetic field.
The field-dependence of magnetization (Fig. 1b) shows saturation at fields~1 T at low temperatures, which suggests that the magnetism of Ho 3+ ( 5 I 8 ) ions in langasites is dominated by the two lowest-energy levels split by the crystal field from other states of the J = 8 multiplet. The level splitting in this so-called non-Kramers doublet, Δ, is estimated to be 3.1 K. As discussed below, based on symmetry arguments and magnetization fits, the Ising-type magnetic moments of Ho-ions are oriented perpendicular to the local symmetry axis, i.e., in the X i Y i plane, forming an angle γ ≈ 30 ∘ with the crystallographic c-axis (Fig. 1b). γ is the only free magnetic anisotropy parameter in our model that reproduces well the field dependence of the magnetization along the b * and c-axes (see Fig. 1a) but underestimates the magnetization for H∥a, which is indicative of deviations of the magnetic moments from the X i Y i -plane (see discussion below). Figure 2 shows the angular dependence of the electric polarization along the c-axis, P c , for the magnetic field vector, H, rotated in the ac-and ab-planes. The complex angular dependence is indicative of a high-order ME effect resulting from the C 3 symmetry of langasite. Although the expansion of the magnetically induced electric dipole of the Ho-ion in powers of H begins with terms ∝ H 2 , the bilinear ME effect is cancelled in the sum over the dipole moments of the Ho dopants in the three La sublattices.
The lowest-order phenomenological expression for P c allowed by C 3 and C 2 symmetries is where we use the Cartesian coordinates withx ¼â,ŷ ¼b Ã ,ẑ ¼ĉ and a 4 (T) is a material constant. Fourth-order ME effect in crystals with trigonal crystal symmetry was discussed earlier [34][35][36] . Equation (1) implies that P z / sin 2 θ sin 2θ cos 3φ, where θ and φ are, respectively, the polar and azimuthal angles describing the magnetic field direction. This angular dependence is in good agreement with the results of experimental measurements at high temperatures shown in Fig. 2a, whereas at low temperatures the experimental θ-dependence becomes more complex, indicating substantial higher-order contributions to the ME effect.
The fourth-order ME effect described by Eq. (1) gives zero electric polarization for a magnetic field in the basal ab plane. However, the experimentally measured electric polarization shown in Fig. 2b is of the same order as that for out-of-plane fields. The electric polarization shows a periodic dependence on   2 Angular dependence of the static electric polarization along the c-axis. a P c vs the polar angle θ between the applied magnetic field, H, and the c-axis, for H rotated in the ac-plane at several temperatures. b P c vs the azimuthal angle φ for H rotated in the ab-plane. The modulation with the 60 ∘ period is indicative of a high-order ME effect (see Eq. (2)). Experimental data are shown with symbols, and solid lines are calculated within the model described in the text.
the azimuthal angle, φ, with the period of 60 ∘ apparently resulting from a sixth-order ME effect. The lowest-order expression for the electric polarization induced by an in-plane magnetic field is, indeed, of sixth order in H: The sixth-order ME response has not been discussed earlier. The right-hand side of Eq. (2), found using symmetry arguments [34][35][36] , is the lowest-order polynomial of H x and H y that transforms as P z . The saw-tooth shape of the angular dependence of the polarization at low temperatures ( Fig. 2b) indicates contributions of ME effects of yet higher orders.
Equations (1) and (2) imply that the electric polarization is proportional to the 4th or 6th power of the applied magnetic field, depending on the orientation of H. Instead, at low temperatures we obseve a nearly linear dependence of P z for any field direction (see Fig. 3). Importantly, phenomenological Eqs. (1) and (2) only hold for weak magnetic fields, |μ 0 H| ≪ max(Δ, k B T), μ 0 being the saturated magnetic moment of a Ho 3+ ion (these equations are obtained in Methods by expansion in powers of H). At low temperatures, their validity is limited to fields ≲ 1 T (see Fig. 1b). The saturation of the Ho magnetic moment at higher magnetic fields gives rise to higher-order harmonics in the observed angular dependence of the electric polarization.
Counterintuitively, the onset of a more complex θand φdependence correlates with the crossover to a linear dependence of P z on the magnitude of the magnetic field (Fig. 3). In what follows we discuss the theory behind this unusual ME effect.

Theory
In the absence of inversion symmetry, the rare-earth ions interact with the electric field, E, through an effective dipole moment operator,d, which in the subspace of the ground-state multiplet can be expressed in terms of the quadrupole and higher-order magnetic multipole moments of the ions 37 . The spectrum of the low-energy states depends then on both electric and magnetic fields, which gives rise to a ME response of the rare-earth ions. The description of the ME effect in Ho-langasite can be simplified by projecting the Hamiltonian describing the ME behavior of the ground-state multiplet on the subspace of the two states, A j i and B j i, forming the non-Kramers doublet (see "Methods").
Alternatively, the Hamiltonian describing the ME coupling of the non-Kramers doublet can be directly deduced from symmetry considerations. Within the unit cell of Ho-langasite, there are three equivalent rare-earth positions with local C 2 symmetry axes along theẐ i (i = 1, 2, 3) directions related by 120 ∘ -rotations around the crystallographic c axis. At site 1,Ẑ 1 k a,Ŷ 1 k c (along the C 3symmetry axis) andX 1 ¼ ½Ŷ 1 Ẑ 1 . In the local coordinate frame, the ME Hamiltonian is, (the index i = 1, 2, 3 labeling the rare-earth ion is omitted for simplicity). Equation (3) is a general symmetry-allowed expression provided one of the states forming the doublet is even and another is odd under . For E X = E Z = 0 and E Y = E c , the total effective Hamiltonian in the subspace spanned by the A j i and B j i states is: where (σ x , σ y , σ z ) are Pauli matrices, Δ is the crystal-field splitting, μ = (μ X , μ Y , 0) and d Y are, respectively, the (transitional) magnetic and electric dipole moments of the Ho 3+ ion. The form ofĤ eff is constrained by 2 Z and time-reversal symmetries of the Ho ions. In addition, 3 c symmetry implies that the (real) coefficients, d Y , μ X , μ Y and g YZ are the same for all Ho-sublattices. The free energy of the doublet, f i (i = 1, 2, 3), can now be easily calculated and the magnetically induced electric polarization is given by where n Ho is the density of Ho ions and M i is related to the average magnetic moment of the i-th Ho ion by 〈μ i 〉 = (μ X , μ Y , 0)M i and equals where +ϵ i and −ϵ i are energies of the two states forming the  section and refs. [38][39][40][41]. The resulting magnetic field-dependence of the electric polarization (dashed lines in Figs. 2 and 3) is in reasonable agreement with experimental data. The largest discrepancy between theory and experiment is found in the angular dependence of the electric polarization shown in Fig. 2. The main reason for this deviation is the intrinsic disorder in the langasite crystal structure (Fig. 1). The four positions in the local surrounding of rare-earth ions are randomly occupied by two Si and two Ga ions 24,32 resulting in six inequivalent Ho ions in each sublattice. The concomitant lattice distortions affect magnetic anisotropy and ME interactions of the ions. Importantly, our simplified model based on average symmetry of rare-earth ions captures the basic physics of the ME response of Ho-langasite. Figure 3 demonstrates that the electric polarization is a linear function of the applied field strength in the broad range of temperatures and external magnetic fields. The region of linear effect agrees well with the regime of saturated magnetic moments shown in Fig. 1b. This can be understood by noting that Eq. (5) can be obtained from the ME coupling, Àg YZ E c H Zi M i , on the i-th Ho sublattice. At low temperatures, βϵ i ≫ 1, M i grows with the applied magnetic field and approaches unity at rather weak fields, |μ 0 H|~Δ, above which the magnetic moment saturates and the ME coupling becomes linear.

DISCUSSION
Importantly, the ME coupling in Eq. (3) originates from the admixture of higher-energy excited states from the J = 8 multiplet to A j i and B j i in the presence of electric and magnetic fields (see Methods). The linear ME effect is a consequence of a relatively large energy separation, W, of the non-Kramers doublet from the higher-energy states: it is observed for Δ ≪ |μ 0 H| ≪ W, when the dipole magnetic moment of the doublet is saturated while the admixture of higher-energy states remains small. Rear-earth ions with Kramers doublets can also show linear ME response, which makes the whole series of the rare-earth dopands potentially interesting for studies of this unusual phenomenon.
There is an additional ME coupling that results from the matrix elements of the electric and magnetic dipole operators between the A j i and B j i states in Eq. (4) and does not involve high-energy excited states. This coupling, however, is proportional to second and higher powers of the electric field and does not contribute to the magnetically induced electric polarization measured in zero applied electric field, which can be traced back to the fact that the electric dipole is even and magnetic dipole is odd under time reversal.
Next we explain how the linear dependence of P c on the magnetic field strength can be compatible with the complex dependence on the direction of the field (see "Methods" for more details). The six-fold periodicity of the polarization with the azimuthal angle is related to the presence of three rare-earth sublattices. Each sublattice separately gives rise to oscillations with the period of 180 ∘ (see Eq. (17)), since the Ising-like Ho magnetic moment changes sign as the magnetic field rotates in the ab-plane. Similar twofold periodicity is observed in the ME Co 4 Nb 2 O 9 with Heisenberg spins, where the Néel vector describing an antiferromagnetic ordering of Co spins rotates together with the magnetic field vector 42 . In addition, the ME effect in Holangasite is a sum of the responses of the three kinds of Ho positions with magnetic moments rotated through 120 ∘ around the c axis, which gives rise to oscillations with the period of 60 ∘ .
Ho-langasite can show the inverse effect-an electrically induced magnetization ∝ E c -even though the ME coupling is an even function of H. The sensitivity to the magnetic field sign, important for applications of ME materials, is recovered under a bias magnetic field. An additional, electrically induced magnetization changes sign under the electric field reversal. If the bias field is strong enough to saturate Ho magnetic moments, the inverse ME effect is linear (see "Methods" for more details).
We have shown that the magnetic and ME properties of Holangasite may be understood taking into account the two lowest crystal-field levels of Ho 3+ and the interplay of the local and global symmetries. Our theory reconciles a complex oscillating dependence of the polarization on the magnetic field direction with a linear dependence on its magnitude. The electric polarization is enabled by the absence of inversion symmetry at the local Ho positions. In weak magnetic fields, when the Zeeman energy is small compared to the energy splitting in the non-Kramers doublet, the polarization is described by Eqs. (1), (2) and is highly non-linear in H. It becomes linear in the strong-field limit, when the induced magnetic moment of Ho ions saturates, but it remains a complicated function of the magnetic field direction.
Importantly, the doping of non-centrosymmetric materials by magnetic rare-earth ions at low-symmetry positions provides a route to a linear ME response that does not require special magnetic orders and crystal lattices, as it occurs already in the paramagnetic state. Oscillatory dependence of magnetic and ME responses on the direction of the applied magnetic field is governed by crystal symmetry and is a generic property of such systems.

METHODS Experimental
Single crystals of doped holmium langasite Ho x La 3−x Ga 5 SiO 14 were grown by the Czochralski technique. The crystals were proven to be right-handed single domain 24,32 . The handedness of langasites is likely a result of the growth procedure 24 . The orientation and crystallinity of the samples have been controlled using Laue-analysis. The actual concentration of holmium has been obtained using the total reflection X-ray fluorescence spectroscopy. The sample showed homogeneous polarization rotation of about 5 ∘ for the green light, which agrees with literature data 43 .
For the electric polarization experiments a plane-parallel plate with thickness d = 2.13 mm and area A = 106 mm 2 has been cut with the flat surface perpendicular to the c-axis (c-cut). Silver paste electrodes were applied to two surfaces of the sample parallel to the ab-plane. The polarization was measured using an electrometer adapted to a Physical Property Measuring System, with magnetic fields of up to 14 T and temperatures down to 2 K. By changing the orientation of the sample in the cryostat, the direction of the magnetic field relative to the crystal axes was adjusted. The magnetization was measured using a commercial Vibrating Sample Magnetometer with magnetic fields up to 6 T. In these experiments, smaller samples from the same batch were used.

Microscopic theory of ME effect in Ho-langasite
The symmetry of the rare-earth environment in langasites is described by the point group C 2 . This cyclic group has two one-dimensional representations with respect to the rotation through an angle π around the twofold symmetry axis (local Z i -axis): the symmetric, A, and antisymmetric, B, representations. Ho 3+ with an even number of f-electrons is a non-Kramers ion. The crystal field with C 2 symmetry completely removes the degeneracy of the ground 5 I 8 multiplet that splits into 2J + 1 = 17 singlets with the wave functions transforming according to either A or B irreducible representations of the C 2 group.
In the absence of inversion, the rare-earth ions interact with the electric field, E, through an effective dipole moment operator,d, which in the space of the Ho 3+ ground-state multiplet can be expressed in terms of the quadrupole and higher-order magnetic multipole moments of the ions 37 : whereQ αβ ¼ 1 2 ½Ĵ αĴβ þĴ βĴα À 2 3 JðJ þ 1Þδ αβ is the quadrupole moment operator,Ĵ ¼ ðĴ X ;Ĵ Y ;Ĵ Z Þ being the total angular momentum operator of the Ho 3+ lowest-energy multiplet, a, b and c are phenomenological constants and the higher-order multipoles are omitted. Equation (7) is written in the local coordinate frame, (X i , Y i , Z i ), and the index i = 1, 2, 3 labeling the rare-earth ion is omitted for simplicity. The phases of the wave functions forming the ground-state multiplet can be chosen in such a way that the matrix elements of the electric dipole (even under time-reversal) are real, whereas the matrix elements of the magnetic dipole (odd under time-reversal) are all imaginary.
The ground state of the Ho 3+ ion in the crystal-field potential, V cf , is a quasi-doublet (two close-by singlets with the gap Δ~2 cm −1 ) which are separated by the energy W~30-50 cm −1 from the excited crystal-field states. The wave functions, A j i and B j i, of the quasi-doublet can belong either to the same or to different representations. As discussed in the main text, the c-component of the magnetization measured in the experiment is nonzero only in the latter case.
The response of Ho 3+ ions to the applied electric and magnetic fields is described byV ¼ Àd Á E þ μ B g LĴ Á H, where the first term describes the coupling of the effective dipole operator Eq. (7) to the electric field and the second term is the Zeeman coupling, μ B and g L = 5/4 being, respectively, the Bohr magneton and the Lande factor. To describe the ME response of the Ho-langasite, we project the total Hamiltonian,Ĥ ¼V þV cf , on the quasi-doublet subspace: where i, j = 1, 2 label the states of the non-Kramers doublet, A j i and B j i, and the effect of the higher-energy states from the 5 I 8 multiplet (k = 3, 4, …, 17) is treated in the second order of perturbation theory assuming jhkjVjiij ( W, and W k ¼ E k À ðE1 þ E2Þ 2 . We note that this approach can also be used to analyze the rare-earth contribution to magneto-elastic and optical properties. The matrix elements ofV in the quasi-doublet subspace are constrained by C 2 symmetry (cf. Eq. (4)): Due to the 3 c symmetry of the langasite crystal, the real coefficients d 0 , d X , μ X , etc. are the same for all Ho sites, i.e., are independent of i. The phases of the wavefunctions A j i and B j i are chosen such that the matrix elements of the time-even electric dipole operator involve real matrices, I and σ x , whereas the time-odd magnetic dipole operator is proportional to σ y with imaginary matrix elements.
In the last term of Eq. (8) we only leave the off-diagonal parts proportional to the product of the matrix elements ofd Á E and μ B g LĴ Á H giving rise to the ME coupling, (cf. Eq. (3)). Here, the coefficients g XY , g YX , etc. are real numbers. This interaction is invariant under C 2 and time-reversal symmetries, and since it involves both the electric and magnetic dipole operators, the matrix elements ofĤ me in the quasi-doublet subspace are imaginary. Adding all contributions to the effective two-level Hamiltonian and assuming E Xi ¼ E Zi ¼ 0 and E Yi ¼ E c (i = 1, 2, 3), we obtain Eq. (4). The energies of the states forming the non-Kramers doublet at the i-th Ho 3+ site, obtained by diagonalization of the Hamiltonian Equation (4), equal ±ϵ i , where The free energy is then given by where n Ho is the density of Ho ions and k B is the Boltzmann constant. The average magnetic moment measured at zero applied electric field is where β ¼ 1 kBT , μ 0 is the magnetic moment of the non-Kramers doublet and m i is a unit vector defined by The relation between the unit vectors in the hexagonal, Cartesian and the three local frames (see Fig. 1a Note that since the electric polarization is measured at zero applied electric field, d Y that enters into the expression for the energy levels Eq. (11), drops out from the expression for P c , in which The energy splitting in the non-Kramers doublet by both electric and magnetic fields (see Eqs. (9) and (11)) does not result in the ME coupling linear in electric field, because the electric and magnetic dipoles transform with opposite signs under time reversal. The observed ME effect originates solely from the admixture of higher-energy states in applied electric and magnetic fields (the last term in Eq. ( (8)). In addition, Ho ions carry a "builtin" electric dipole moment along the local Z i axis in the ab plane, which results from the first term in Eq. (9) and is independent of the magnetic field. These electric dipoles on three Ho sites are rotated through 120 ∘ around the c axis with respect to each other and give no net contribution to the electric polarization.
In the limit of strong fields and low temperatures, |μ ⋅ H i | ≫ Δ, k B T, when the magnetic moment saturates, the electric polarization, shows a complex saw-tooth dependence on the direction of the magnetic field in the ab plane (due to the sign-functions) and, at the same time, a simple linear dependence on the strength of the magnetic field, in agreement with our experimental observations. In the opposite limit of weak fields and high temperatures, k B T ≫ |μ 0 H|, Δ, P c % n Ho g YZ μ 3 0 H 4 4 k B T ð Þ 3 sin γcos 2 γ cos θsin 3 θ cos 3φ; (18) where H is the magnitude of the magnetic field and γ ¼ arccosðμ Y =μ 0 Þ, which agrees with Eq. (1) obtained by symmetry and gives a 4 (T) ∝ T −3 . At zero temperature, ðk B TÞ 3 in Eq. (18) is replaced by Δ 3 /12. For magnetic field applied in the ab plane, the electric polarization at high temperatures, P c % n Ho g YZ μ 5 0 H 6 240 k B T ð Þ 5 cos 5 γsin 6 θ sin 6φ ð Þ; is proportional to H 6 and oscillates six times as φ varies by 2π, in agreement with Eq. (2) (a 6 ∝ T −5 ). At zero temperature, ðk B TÞ 5 ! Δ 5 =90. The coefficients a 4 and a 6 in Eqs. (1,2) can be estimated from the lowfield slopes of experimental data shown in Fig. 3. As demonstrated in the insets to Fig. 3, they indeed follow the expected a 4 ∝ T −3 and a 6 ∝ T −5 behavior.
The ME coupling ∝ E c , which leads to Eq. (16) for the magnetically induced polarization, also gives rise to the electrically induced magnetization, δM ¼ À ∂f me ∂H . Figure 4a shows the magnetization at E c = 0 (Eq. (13)) plotted vs the angle φ describing the orientation of the magnetic field in the ab plane. In general, M is not parallel to H. In the saturation regime, it depends on the sign of the magnetic field projection on the easy axes of Ho ions, but not on the magnetic field strength. Figure 4b shows the electrically induced magnetization, δM, divided by E c . M and dM=dE c ð Þ Ec ¼0 have opposite symmetries with respect to φ → − φ: M x is a symmetric function of φ and M y , M z are antisymmetric, whereas for dM/dE c it is the other way around.
Interestingly, the inverse ME response is amplified when the magnetic field is normal to the easy axis of one of the Ho-sublattices. At zero temperature and for μ 0 H ≫ Δ, δM shows sharp peaks at (m i ⋅ H) = 0,

DATA AVAILABILITY
The data that support the findings of this study are available from the authors on reasonable request, see author contributions for specific data sets.