Magnetoelectric coupling of domains, domain walls and vortices in a multiferroic with independent magnetic and electric order

Magnetically induced ferroelectrics exhibit rigidly coupled magnetic and electric order. The ordering temperatures and spontaneous polarization of these multiferroics are notoriously low, however. Both properties can be much larger if magnetic and ferroelectric order occur independently, but the cost of this independence is that pronounced magnetoelectric interaction is no longer obvious. Using spatially resolved images of domains and density-functional theory, we show that in multiferroics with separately emerging magnetic and ferroelectric order, the microscopic magnetoelectric coupling can be intrinsically strong even though the macroscopic leading-order magnetoelectric effect is forbidden by symmetry. We show, taking hexagonal ErMnO3 as an example, that a strong bulk coupling between the ferroelectric and antiferromagnetic order is realized because the structural distortions that lead to the ferroelectric polarization also break the balance of the competing superexchange contributions. We observe the manifestation of this coupling in uncommon types of topological defects like magnetoelectric domain walls and vortex-like singularities.

T he different interdependence of magnetic and electric order in multiferroics with jointly ('type-II') and separately ('type-I') emerging magnetic and ferroelectric order 1 manifests prominently on the level of domains. In the former case, magnetic and electric domains are one-to-one correlated, and every magnetic domain wall is a ferroelectric domain wall, too. This leads to potentially useful phenomena like control of magnetic domains by electric fields, and vice versa 2 . In contrast, in the latter case, the magnetic and ferroelectric domains and domain walls are not intrinsically linked to each other because of the independence of the associated order parameters. Magnetoelectric coupling phenomena may still be present, but they are not mandatory. The type-I multiferroic systems of hexagonal ferrites and manganites are model cases for such magnetoelectric coupling phenomena. In ferrites, it was shown that the leading-order, that is, linear magnetoelectric effect dominates the magnetoelectric coupling 3,4 . However, in the isostructural manganites, disregarding any magnetic-field-induced phases and rare-earth ordering 5 , the linear magnetoelectric effect is symmetryforbidden 6 . The demonstration of strong magnetoelectric coupling phenomena in hexagonal manganites, despite this absence of the linear magnetoelectric effect, is at the heart of our work.
We show in experiment and theory that the superexchange interaction drives a pronounced microscopic magnetoelectric interaction. The associated bulk coupling phenomenon is distinctly different from the coupling between 7 or within 8 domain walls proposed earlier. It has significant consequences for the magnetoelectric domain morphology. We furthermore identify three types of magnetic domain walls with different types of magnetic pseudo-vortices at their meeting points. We thus see that the independence of magnetic and electric order in type-I multiferroics can lead to properties that are not open to type-II multiferroics and can thus be beneficial rather than detrimental to their magnetoelectric functionality.

Results and discussion
Hexagonal ErMnO 3 as representative of the RMnO 3 family with R = Sc, Y, In, Dy-Lu is formed by layers of corner-sharing MnO 5 bipyramids alternating with sheets of Er 3+ ions. A phase transition from the non-polar P6 3 /mmc to the polar P6 3 cm space group occurs at T C ≈ 1430 K 9 . It results from the collective tilting of every three neighbouring MnO 5 trigonal bipyramids whose apical O 2− ions jointly move along the radial direction relative to their central Er 3+ site. The associated lattice mode has K 3 symmetry and is parameterized by an amplitude Q and azimuthal rotation angle Φ, see Fig. 1. The distortion activates a polar displacement of the Er 3+ ions along the c axis, giving rise to improper ferroelectricity with a polarization P / cos 3Φ. This leads to six possible domain states with Φ = n × 60 ∘ (n = 0, 1, ..5), displaying ΔΦ = 60 ∘ and, hence, alternating polarization between neighbouring domains 10,11 . The system is famous for the formation of ferroelectric pseudo-vortices [12][13][14] , which arise due to the breaking of the effectively continuous symmetry of Φ close to the Curie temperature 15-17 . At T N = 77 K, the antiferromagnetic ordering of the Mn 3+ spins occurs 18 according to the magnetic K 2 representation of the non-polar P6 3 /mmc structure 19 . The spin angle between nearest Mn 3+ neighbours in the basal ab planes is 120 ∘ with an arrangement of neighbouring planes along c as depicted in Fig. 2. We describe the local direction of the Mn 3+ spins with the angle Ψ, where a change of Ψ corresponds to an in-phase rotation of all the spins in the ab plane, but with an opposite sense of rotation in the upper and lower half of the unit cell 19 . For Q = 0, we confirm, using density functional theory (DFT), that the energy of the magnetic order in Fig. 2b is independent of Ψ, as required by symmetry. For Q ≠ 0, however, coupling of the magnetic mode to the structural mode described by a free-energy contribution is permitted, where the index Ψ emphasizes that the magnetic order adapts to the already established trimerized-polar order. The reduction in symmetry changes the energy landscape from a continuum of equally low-energy states for Ψ − Φ to only two possible minimum-energy solutions, Ψ − Φ = ±90 ∘ . This constraint may appear surprising on first glance since the centrosymmetry of the Mn 3+ sublattice is not affected by the distortive-ferroelectric transition 10 . Microscopically, however, a correlation between the spin and the polar lattice is established because the magnetic exchange between the Mn 3+ ions, which determines the antiferromagnetic order, is mediated by the O 2− ions 20,21 whose displacement from their original high-symmetry positions is part of the of the structural distortion that leads to the ferroelectric order. Hence, the superexchange projects the polar order into the magnetic system 3 , thus establishing a bulk magnetoelectric coupling. It breaks the continuous symmetry in (Ψ − Φ) and centres the high-symmetry point of the magnetic lattice to that of the trimerized-polar lattice. The presence of this coupling is all the more striking as the linear magnetoelectric effect describing the emergence of a magnetization (electric polarization) proportional to an applied electric (magnetic) field is symmetry-forbidden in ErMnO 3 6 . In this aspect, the hexagonal manganites are strinkingly different from the hexagonal ferrites. In LuFeO 3 , the linear magnetoelectric effect is allowed and involves contributions from a magnetization of the Fe system along the hexagonal axis 3,4 , which is absent in ErMnO 3 .
Symmetry analysis and DFT show that multiferroic structures Fig. 2a have the same energy; these correspond to two antiferromagnetic structures distinguished by a relative reversal of all Mn 3+ spins as sketched in Fig. 2b. To see how the degeneracy of the two structures in Fig. 2b influences the domain formation, we experimentally investigate the spatial distribution of ferroelectric and antiferromagnetic domains in ErMnO 3 . To date, there has been a single such investigation in the system of hexagonal manganites. Zero bulk magnetoelectric coupling was assumed 7,22,23 , and the peculiar vortex-like arrangement of the six associated domain states was not considered and could not be spatially resolved [12][13][14] .
We probe the spatial distribution of the domains by optical second-harmonic generation (SHG), yet with ten times higher optical resolution than in the earlier experiments (see 'Methods'). SHG denotes the doubling of the frequency of a light wave in a material. It is a highly symmetry-sensitive process and can therefore distinguish between the different types of long-range order and domain states 24,25 . ErMnO 3 features two types of SHG contributions 7 , a ferroelectric one proportional to the electric polarization P / cos 3Φ and a multiferroic one proportional to the product PL, where L / sin 3Ψ is chosen such that it phenomenologically associates the two antiferromagnetic domain states depicted in Fig. 2b to opposite signs. As detailed in 'Methods' and the Supplementary Material, interference of the contributions / P and / PL leads to an antiferromagnetic net SHG contribution / L. This approach 24 allows us to distinguish domains with opposite signs of P, L or PL as regions of different brightness, or the domain walls themselves as dark channels 26 .
In Fig. 2c-e, we see spatially resolved SHG images showing the distribution of the domains for P, L and PL in the same area of a c-oriented ErMnO 3 sample (see 'Methods'). First, the SHG light confirms the relation Ψ − Φ = ±90 ∘ derived from Eq. (1) because its polarization reproduces the associated magnetic symmetry 6,7 . The domain configuration in P yields the familiar distribution of six domains of alternating polarization (see Fig. S2) arranged around a point in which all these domains meet. Strikingly, the domain configuration in L reproduces this distribution, showing an arrangement of six domains of alternating orientation of the antiferromagnetic order parameter around a vortex-like meeting point. The distribution of PL is even more striking as it reveals an area of approximately homogeneous brightness without any domain walls, pointing to a PL single-domain state. The system preserves the value of (Ψ − Φ) through simultaneous changes of Φ and Ψ by the same value of ±60 ∘ when crossing a P or The inset with the coordinate system defining Ψ, in the bottom right, takes the Mn 3+ ion in the centre of the three Er 3+ ions in a as origin. c Magnetoelectric bulk coupling energy ΔE as function of the Mn 3+ spin angle Ψ. Calculated using DFT in the absence (orange) and presence (blue) of the distortive-ferroelectric order. d-f Spatially resolved distribution of SHG intensity on the same region of a c-oriented ErMnO 3 sample for the SHG coupling to P, L and PL, respectively. Black lines-highlighted with dashed blue lines-in d indicate ± P domain walls, and dark and bright regions in e distinguish ± L domains. In either case, the brightness difference results from SHG interference processes 24,26 . For the assignment of ± P domains in d, see Fig. S2. SHG images were taken at room temperature (P) or 20 K (L, PL). Scale bar in f is 25 μm.
L domain wall. In this sense, PL is associated to a configuration of hyperdomains at least an order of magnitude larger (see Fig.  S3) than the P and L domains and not showing their topologicalvortex-like distribution.
The preservation of (Ψ − Φ) immediately disproves the model of a piezomagnetic coupling between strained ferroelectric and locally magnetized antiferromagnetic domain walls proposed earlier 7,22 , as the latter would promote the largest possible wall magnetization, which would be accomplished by a change of Ψ by 180 ∘ across the wall.
The preservation of (Ψ − Φ) across the domain walls can be captured in the free energy by adding a gradient term describing the magnetic domain wall energy to Eq. (1) according to where s, A > 0 are coefficients that we calculate from first principles (see 'Methods'). As stated earlier, the second term is minimized by having Ψ − Φ = ±90 ∘ in each domain. The energy cost from the first term increases with the total rotation angle ΔΨ across the domain wall.
Regarding the change of Ψ across the domain wall we have the following possibilities. (i) ΔΨ = ±60 ∘ if Ψ tracks ΔΦ = ±60 ∘ across the wall; in this case both P and L change sign and PL is preserved. (ii) ΔΨ = ∓120 ∘ if only P (and with it PL) changes sign while L is preserved. The PL-preserving domain wall has the lower energy cost, consistent with our observations. (iii) As a third possibility, we may have a purely magnetic domain wall within a single ferroelectric domain in which only L (and with it PL) changes while P is preserved. This requires ΔΨ = 180 ∘ to satisfy Ψ − Φ = ±90 ∘ and is therefore also unfavourable.
We now investigate the occurrence of the three possible types of magnetic domain walls with ΔΨ = ±60 ∘ , ∓120 ∘ or 180 ∘ in the same region of the sample after consecutive annealing cycles above T N . The microscopy images in Fig. 3a-d show the ferroelectric ± P domains and the antiferromagnetic ± L domains along with the analysis in terms of the distribution of Φ and Ψ in Fig. 3e-h. Note that the direct tracking of the orientation of Φ and Ψ across a section through the domain wall is prohibited by the optical resolution limit of 1 μm.
The energetically preferred ΔΨ = ±60 ∘ walls are most common and always appear at the same location because of their clamping to the ΔΦ = ±60 ∘ walls, which are not affected by the heating above T N . In Fig. 3b, f, these are the only domain walls present. In addition, Fig. 3c, g reveals a single ΔΨ = ∓120 ∘ wall. Even though L does not change across this wall, it nevertheless represents a magnetic domain wall because Ψ has to change along with Φ in order to retain L. Finally, Fig. 3d, h shows a single ΔΨ = 180 ∘ wall across which the brightness of each antiferromagnetic domain is reversed. This is the only unclamped wall, across which only L and Ψ change, while P and Φ do not.
The coexistence of the three types of magnetic walls with the topologically protected network of ferroelectric vortex-like domains leads to additional magnetic pseudo-vortices. According to the analysis in Fig. 3f-h we have three types of these. The domain walls with ΔΨ = ±60 ∘ and ΔΦ = ±60 ∘ coincide and therefore exhibit sixfold magnetic pseudo-vortices in both Ψ and L (Fig. 3f). The meeting of a ΔΨ = ±60 ∘ wall and a ΔΨ = 180 ∘ wall can occur as an intersection, which establishes a fourfold magnetic pseudo-vortex in both Ψ and L (Fig. 3h). Alternatively, the meeting of a ΔΨ = ±60 ∘ wall and a ΔΨ = 180 ∘ wall can occur as junction with a ΔΨ = ∓120 ∘ wall, which establishes a threefold pseudo-vortex in Ψ, but not in L (Fig. 3g).
Finally, to confirm that this interpretation is consistent with the configuration of domains and domain walls seen in Fig. 3, we studied the distribution of the ErMnO 3 domains in phase-field simulations, based on our model free energy in Eq. (2) and our parameters calculated by DFT. Figure 4 shows the resulting calculated domain configurations in Φ, Ψ and (Ψ − Φ) and their observable projections P, L and PL. The results are in excellent agreement with the measured data in Figs. 2 and 3. In particular, all the magnetic pseudo-vortices discussed above are obtained and the PL hyperdomains are an order of magnitude larger than the domains in P and L alone. This confirms that the magnetic domain morphology in ErMnO 3 is indeed a direct consequence of the microscopic magnetoelectric bulk coupling.
In conclusion, we see that despite the absence of the linear magnetoelectric effect in the type-I-multiferroic hexagonal manganites, these compounds host a pronounced microscopic bulk magnetoelectric coupling. For our model compound, hexagonal ErMnO 3 , we show that the interaction between the spins and the lattice occurs because of superexchange, specifically because the coupling between the magnetic Mn 3+ ions is mediated by the O 2− ions, which have different positions in different polar domains. The hidden bulk magnetoelectric coupling in this type-I multiferroic leads to phenomena not open to type-II multiferroics, such as a rich variety in the types and topology of domain walls. Structural distortions determine the magnetic order of many multiferroics. Therefore, a microscopic magnetoelectric bulk coupling, demonstrated here for ErMnO 3 , should be present in many of these. In particular, we see that the independence of magnetic and electric order in type-I multiferroics is beneficial rather than detrimental to their magnetoelectric functionality because of the greater freedom in coupling the magnetic to the electric domains.

Methods
Sample preparation. We used an ErMnO 3 bulk crystal of about 1 × 2-mm 2 lateral dimensions. The sample was grown by the flux technique and oriented using a Laue x-ray diffractometer. The sample was cut perpendicular to the hexagonal c axis using a diamond saw, thinned down to~30 μm by lapping with Al 2 O 3 powder of 9-μm grain size and chemo-mechanically polished with a colloidal silica slurry until a root-mean-square roughness below 50 nm was reached. Ferroelectric domains with tens of micrometres in size are obtained by heating the crystal across T C = 1429 K and re-cooling it at a rate of about 0.01 K/min 9 .
Optical second-harmonic microscopy. A transmission SHG setup described in detail elsewhere 24 is used to acquire the spatially resolved SHG images. We use three experimental configurations (see Table 1) to disentangle the different SHG contributions in ErMnO 3 7,27 with k denoting the propagation direction of the light.
The SHG signal fromχðPÞ is usually recorded at room temperature using an angle of 20 ∘ between the wavevector and the c axis. Spatial maps reveal the ferroelectric domain configuration. At the domain walls, the opposite sign of the order parameter on either side introduces a relative phase shift of 180 ∘ between the respective SHG waves, so that destructive interference distinguishes the domain walls as black lines. The association of opposite P domains is retrieved from linear phase-contrast microscopy measurements as detailed in Fig. S2. The SHG signal fromχðPLÞ is recorded below T N . Spatial maps reveal the multiferroic hyperdomain configuration. The SHG signal fromχðLÞ is obtained from the interference of the Pand PL-related SHG waves. A phase shift between the two contributions and, thus, a change in the total SHG intensity resulting from the interference of the two contributions, can only occur with a change in L. A change in P would affect both SHG contributions and therefore retain their relative phase and the resulting SHG intensity. Images were corrected for the variation of intensity across the laser-beam profile.
For the laser excitation we use a Coherent Ellite Duo laser system (1.55-eV fundamental wave,~120-fs pulse duration, 8-mJ pulse energy, 1-kHz repetition rate). An optical parametric amplifier converts the emission to a photon energy of 1.23 eV and a pulse energy of about 30 μJ. A ×20 long-working-distance microscope objective is used to collect the SHG light emitted from the sample and project it onto a liquid-nitrogen-cooled Jobin Yvon Back Illuminated Deep Depletion CCD camera. Optical filters in the optical path suppress any unwanted frequency components.
First-principles calculations. For our first-principles calculations, we use the PBEsol+U approximation of the exchange correlation potential as implemented in the VASP code [28][29][30][31][32] . For the calculation of the magnetic energy landscape, we use U = 3 eV. In the Er 3+ pseudopotential, the 4f electrons are treated as core electrons, and therefore rare-earth magnetism is neglected. We use a k-point mesh of 6 × 6 × 4 and a plane-wave cutoff energy of 800 eV. All the calculations include spin-orbit coupling. To calculate the energy landscape (Fig. 2), we performed a structural relaxation for each spin configuration with a threshold force of 10 −4 eV/Å 2 for each atom. We neglect the small out-of-plane tilting that is symmetry-allowed for some of the configurations 19 since preliminary calculations showed that such a tilt only has a very small influence on the total energy. To calculate the free energy introduced in the main text in Eq. (2), the anistropy parameter (A) is obtained by performing a least squares fit to calculated spin configurations. The gradient parameter s ¼ 1 in Eq. (2) is calculated using a long-wavelength expansion of the magnetic order parameter, where E KS (k) is the Kohn-Sham energy for the structure with a magnetic order modulated with a wavevector k. Numerically this was done using the spin-spiral implementation in VASP. Normalizing the order parameter Q such that |Q| = 0 in the high-symmetry phase and |Q| = 1 in the ground state, we obtain the following parameters for Eq. (2): s = 742 meVÅ 2 for the gradient parameter and A = 2.13 meV for the anisotropy parameter (per 30 atoms). We note that this is a simplified version of the free energy put forward by Artyukhin et al. 11 , including only magnetic configurations of the magnetic K 2 representation. We also calculate the free-energy parameters for ErMnO 3 using the full expression 11 : þ a sin 2 ðψ 1 À ΦÞ þ sin 2 ðψ 2 À ΦÞ Â Ã À C þ cosðψ 1 þ ψ 2 À 2ΦÞ À C À cosðψ 1 À ψ 2 Þ; ð3Þ for which we obtain s = 742 meVÅ 2 , a = 0.027 meV, C + = −1.05 meV, C − = 0.0023 meV (per 30 atoms).   Order-parameter sensitivity Phase-field modelling. We use the Landau expansion for the free energy FðQ; Φ; P; ψ 1;2 Þ of the hexagonal manganites introduced by Artyukhin et al. 11 to derive the Ginzburg-Landau (GL) equations 33 for the structural and magnetic order parameters. We use the values given by Artyukhin for the couplings of Q, Φ and P and values obtained with our DFT simulations for the relation of ψ 1,2 to Φ and ψ 2,1 (see Eq. (3)). Starting from random inital values for all order parameters, we integrate the GL equations applying a semi-implicit Fourier-spectral solver 34 , using a grid spacing h = 0.1 nm and a time step dt = 0.1. We simulate only one layer in the ab-plane using a square grid consisting of n x × n y = 1024 × 1024 points. This corresponds to a physical length of the simulated system of about 100 nm. Since the magnetic domain walls have an expected width on this order 11 , we have to decrease this value in our simulation to obtain realistic domain patterns. We achieve this by increasing the coupling parameter a in Eq. (3) toã ¼ 10 3 a. This reduces the width of magnetic domain walls such that it becomes comparable to the width of the structural walls without changing the topology of the system. Since the ferroelectric phase transition occurs at higher temperature than the magnetic one, we first iterate the trimerization amplitude Q and the phase Φ as well as the electric polarization P for 10 4 steps. We then iterate the spin angles ψ 1,2 on top of this pattern for 10 4 steps, resulting in the domain pattern shown in Fig. 4.
Reporting summary. Reporting Summary Further Information on research design is available in the Nature Research Reporting Summary linked to this article.

Data availability
Data that support our findings are available upon request.

Code availability
The custom code used for the phase-field simulations is available upon request.