Square and rhombic lattices of magnetic skyrmions in a centrosymmetric binary compound

Magnetic skyrmions are topologically stable swirling spin textures with particle-like character, and have been intensively studied as a candidate of high-density information bit. While magnetic skyrmions were originally discovered in noncentrosymmetric systems with Dzyaloshinskii-Moriya interaction, recently a nanometric skyrmion lattice has also been reported for centrosymmetric rare-earth compounds, such as Gd2PdSi3 and GdRu2Si2. For the latter systems, a distinct skyrmion formation mechanism mediated by itinerant electrons has been proposed, and the search of a simpler model system allowing for a better understanding of their intricate magnetic phase diagram is highly demanded. Here, we report the discovery of square and rhombic lattices of nanometric skyrmions in a centrosymmetric binary compound EuAl4, by performing small-angle neutron and resonant elastic X-ray scattering experiments. Unlike previously reported centrosymmetric skyrmion-hosting materials, EuAl4 shows multiple-step reorientation of the fundamental magnetic modulation vector as a function of magnetic field, probably reflecting a delicate balance of associated itinerant-electron-mediated interactions. The present results demonstrate that a variety of distinctive skyrmion orders can be derived even in a simple centrosymmetric binary compound, which highlights rare-earth intermetallic systems as a promising platform to realize/control the competition of multiple topological magnetic phases in a single material.

the main text.

II. Magnetic domains in phases I, II, and IV
In this section, we discuss the magnetic domains in phases I, II, and IV. When the spin texture breaks the symmetry of the original crystal lattice, multiple energetically-degenerated magnetic domains appear. These domains are related to each other by the symmetry elements that are lost during the phase transition process. For example, as discussed in the main text, phase II is the double-Q state satisfying two distinctive fundamental magnetic modulation vectors making an angle of 80° with respect to each other. Its reciprocal-space scattering pattern and corresponding real-space spin texture calculated based on Eq. (2) in the main text are indicated in Supplementary Fig. 3a,d, where the appearance of four fundamental magnetic reflection spots are expected from a single domain (domain 1). Since this spin texture breaks the four-fold symmetry of the original EuAl4 crystal lattice, there will appear another equivalent magnetic domain (i.e. domain 2) that can be obtained by the 90° rotation of domain 1 ( Supplementary Fig. 3b,e). By assuming the equal population of these two domains, the appearance of eight fundamental magnetic reflection spots as shown in Supplementary Fig. 3c is expected. This scenario well explains the observed resonant X-ray scattering and SANS patterns in phase II (Supplementary Fig. 1c and Fig. 2b in the main text). The above argument about the magnetic domain formation can be applied to phases I and IV, in the latter of which broken four-fold and mirror symmetry explains the appearance of eight magnetic reflection spots as seen in Fig. 2d in the main text.

III. Analysis of Q1+Q2 and Q1-Q2 magnetic reflections in phases II and III
In this section, we discuss the Q1+Q2 and Q1-Q2 magnetic reflections in the resonant X-ray scattering profiles. As discussed in the main text, phases II and III are the double-Q states and host the higher-order Q1+Q2 reflections in their SANS profiles (Fig. 2b,c in the main text).
Similar Q1+Q2 reflections are also observed in the corresponding X-ray scattering profiles ( Supplementary Fig. 1c,d). Supplementary Fig. 2d-f indicates the magnetic field dependence of the X-ray scattering intensity, the wavenumber |Q|, and the azimuth angle θQ for the Q1+Q2 reflection. The Q1+Q2 reflections are present in phases II and III but absent in phases I. These features confirm that the observed higher-order reflections originate from the double-Q nature of phases II and III, and not the double-scattering process. Note that the experimentally observed |Q| and θQ values for Q1+Q2 reflections are consistent with the ones theoretically calculated from the fundamental reflections (solid lines in Supplementary Fig. 2e,f).
In Supplementary Fig. 4, schematic illustration of reciprocal-space scattering pattern and the line-scan profiles of the resonant X-ray and neutron scattering measured along the [100] direction are plotted for each phase. In phase II ( Supplementary Fig. 4b,e,h), two magnetic reflection peaks are identified in the line-scan profiles, where the ones with smaller and larger |Q| values represent the Q1-Q2 reflection from domain 1 (Supplementary Fig. 3a) and the Q1+Q2 reflection from domain 2 ( Supplementary Fig. 3b), respectively. They merge into a single Q1+Q2 peak in phase III, where the Q1 and Q2 vectors become exactly orthogonal with respect to each other (Supplementary Fig. 4c,f,i). Note that the peak in phase I represents the fundamental magnetic reflection along the [100] direction, i.e., not the higher-order Q1+Q2 reflection.

IV. Phase boundary between phases I and II
Supplementary Fig. 5a,b shows the magnetic-field variations of the SANS integrated intensities, taken for the boxed regions, Area 1 and Area 2 shown in Supplementary Fig. 5c,d, respectively. In phase I, the intensity is finite only in Area 1, reflecting the fundamental magnetic reflection Q = Q1 along the 〈100〉 reciprocal lattice direction. In the intermediate fields between phases I and II indicated by the gray regions in Supplementary Fig. 5a,b (corresponding to those in Fig. 3 in the main text), the intensity of Area 2 starts to increase, reflecting the transition into phase II, while that of Area 1 gradually decreases. Such anticorrelation of intensity between phases I and II implies that phase I is gradually replaced by phase II. During this phase transition process, the intensity of Area 1 should contain the contributions not only from the Q1+Q2 reflection of phase II but also from the Q1 reflection of phase I. Here, the phases I and II are characterized by the screw spin texture ( Supplementary   Fig. 5e) and the rhombic skyrmion lattice ( Supplementary Fig. 5g), respectively, and the spatial coexistence of these two phases is expected in the intermediate state.
In the main text, we assumed that the spin texture ( )= ( )/| ( )| of phase II can be approximately described as which represents the superposition of screw spin helices characterized by two obliquely arranged magnetic modulation vectors Q1 and Q2, respectively ( Supplementary Fig. 5g . For EuAl4, the amplitude of ( ) should be constant due to the localized character of magnetic moments on Eu 2+ sites, and the associated normalization process leads to the appearance of higher-order ̂( 1 + 2 ) and ̂( 1 − 2 ) components. Note that the experimentally observed scattering intensity for the Q1+Q2 reflection is considerably larger than the Q1-Q2 reflection in phase II ( Supplementary Fig. 4e,h), which implies that the introduction of ̂( 1 + 2 ) component may be necessary. In this case, Eq. (S1) is replaced by which represents the superposition of the screw spin helices characterized by magnetic modulation vectors Q1, Q2, and Q1+Q2. In general, skyrmion lattices are created by pinching off the helical magnetic stripes. To follow this picture, we have chosen the relative phase so that the skyrmion core position defined by Q1 and Q2 locates at the middle of negative mz (antiparallel to external magnetic field) region of Q1+Q2 screw spin modulation. The resultant spin texture is shown in Supplementary Fig. 5f, where elliptically deformed skyrmions are further elongated along the direction perpendicular to Q1+Q2 compared with the one based on Eq. (S1) ( Supplementary Fig. 5g). In either case, the phase II represents the rhombic skyrmion lattice state. The detailed spin texture for phase II and the associated magnetic phase transition process are issues for the future study.

V. Analysis of Hall resistivity profile
In this section, we discuss the magnetic field dependence of Hall resistivity in Fig. 1g in the main text. In general, the Hall resistivity = N + A + T consists of normal Hall term N = 0 0 proportional to H, anomalous Hall term A = proportional to M, and topological Hall term T = 0 eff proportional to the emergent magnetic field eff . R0 and Rs are normal and anomalous Hall coefficients, respectively, and P is the spin-polarization ratio of the conduction electron [S2].
When the rigid band structure is assumed, the anomalous Hall term originating from the intrinsic and skew scattering mechanism is generally described as A ∝ 2 and A ∝ , respectively [S3]. In Supplementary Fig. 6c, H-dependence of A estimated from these formula and experimental M and profiles at 4 K are plotted, which shows considerable deviation from the experimental profile ( Supplementary Fig. 6b). This argument is also applied to the Hall conductivity, , calculated as = Fig. 6a).
At this stage, it is not straightforward to conclude the detailed origin of the observed discrepancy. This is because the electronic structure (and associated carrier density and Fermisurface properties) in rare-earth compounds sensitively depends on the magnetic structure [S4], and therefore the normal and anomalous Hall term N and A should be characterized by more complicated behavior than the above formula. Moreover, recent theory has proposed that the reciprocal-space and real-space Berry phases are closely linked, which gives rise to the socalled chiral Hall effect derived from noncollinear magnetism [S5]. In the present compound, this chiral Hall effect may also provide a sizable contribution to the Hall resistivity.
In principle, the appearance of topological Hall term T can be expected in the skyrmion lattice phases, i.e. phases II and III. Here, the emergent magnetic field eff = Φ 0 Φ originates from the quantum-mechanical Berry phase gained by the conduction electrons passing through skyrmion spin textures, and Φ and Φ 0 = ℎ/ represent skyrmion density and flux quantum, respectively [S2]. If we tentatively assume the typical value of spin polarization ratio P ~ 0.01 previously estimated for the other Gd-based centrosymmetric skyrmion-hosting materials [S6,S7], the expected amplitude of topological Hall term T for phases II and III with the modulation period of 3.5 nm would be in order of |∆ T |~0.16 Ωcm . For the full understanding of the experimental profile and the quantitative evaluation of topological Hall term in EuAl4, however, further theoretical analysis and the detailed information on the electronic structure would be necessary.

VI. SANS patterns for each magnetic phase at 5 K
The SANS data were collected by rotating and tilting the sample such that the magnetic diffraction signal is rotated through the Ewald sphere. The typical SANS patterns at 5 K in Fig.   2a-d in the main text were obtained by combining four distinctive images measured at each rotating (ω) and tilting (χ) angles as shown in Supplementary Fig. 7.

VII. SANS patterns and spin textures for phases IV
In this section, we discuss the detail of the SANS pattern and corresponding spin texture for the temperature (T) -magnetic field (H) region labeled as phase IV. Supplementary Fig. 8 shows ). (S3) It represents the rhombic form of vortex lattice as shown in Supplementary Fig. 8g. At higher-T region, the fundamental modulation vectors become orthogonal to each other and four-fold symmetry is recovered ( Supplementary Fig. 8d,e,f). For the latter state, the polarized neutron scattering results in Fig. 4c in the main text confirms that ̂( ) contains only in-plane spin component. Therefore, its spin texture is also described by Eq. (S3), representing the square vortex lattice state as shown in Supplementary Fig. 8h. As discussed in Supplementary Note IX and Supplementary Fig. 10, our additional theoretical simulation suggests that the above two vortex-lattice spin states with and without four-fold symmetry are energetically almost degenerated. Since we couldn't identify their clear phase boundary from magnetization and electrical transport measurements, we tentatively assigned both of them as phase IV. Some additional discussion on phase IV is also provided in Supplementary Notes VIII and IX.

VIII. Polarized neutron scattering for phases I, IV, V, and VI
To elaborate on the magnetic structures in phases I, IV, V, and VI, additional neutron diffraction measurements with longitudinal polarization analysis were performed at zero field with a triple-axis spectrometer PONTA at JRR-3, Japan. A polarized incident neutron beam with an energy of 13.7 meV was obtained by a Heusler monochromator. The incident neutron The temperature dependences of the q value and the integrated intensity for SF and NSF scatterings for the (q, 0, 0) and (q, q, 0) peaks are summarized in Supplementary Fig. 9c-h. In phases I and V, we observed fundamental magnetic reflections Q1 = (q, 0, 0), and the q value is jumped on the transition between the two phases. For both phases I and V, the intensity of the SF scattering is almost the same as that of the NSF scattering, indicating that ̂( ) possesses both [001] and [010] components normal to Q1. In phases VI and IV, fundamental magnetic peaks appear at Q1 = (q, q, 0), reflecting the reorientation of fundamental magnetic modulation vector. This is in good agreement with the result in Ref. [S8]. The (q, q, 0) magnetic scattering appears only in the NSF channel but not in the SF channel, indicating that ̂( ) possesses [1 ̅ 10] components normal to Q1 but don't have the [001] component. In phase VI, we also identified (q, 0, 0) peak representing the Q1+Q2 reflection. This proves that phase VI is the double-Q state. Since the Q1+Q2 peak appears only in the SF channel but not in the NSF one, it is suggested that ̂( + ) contains the [001] component, but don't have the [010] component in phase VI.
On the basis of the above results, schematic SANS patterns for each phase are illustrated in Supplementary Fig. 9i-l. Since we didn't observe the Q1+Q2 reflections, phases I and V are the single-Q state described by ( ) = 1 cos( • ) + 1 sin( • ) , (S4) corresponding to the screw spin texture with Q1 = (q, 0, 0) ( Supplementary Fig. 9i,j). Phase VI is identified as the double-Q state and can be described by which represents the superposition of four sinusoidal spin orders characterized by magnetic modulation vectors 1 = (q, q, 0), 2 = (q, -q, 0), 1 + 2 , and 1 − 2 (Note that the relative phase of each modulating component is arbitrarily chosen since it cannot be determined from diffraction experiments). The spin texture based on Eq. (S5) is shown in Supplementary   Fig. 9k, which turned out to be identical with the one proposed in a recent theoretical work based on the Heisenberg model on a square lattice with competing interactions (MX-II state in Ref. [S9]). As discussed in Supplementary Note VII, phase IV can be also considered as the double-Q state. Its spin texture at H = 0 can be described by ), (S6) representing the superposition of two sinusoidally modulated spin components with 1 = (q, q, 0) and 2 = (q, -q, 0). The resultant spin texture is shown in Supplementary Fig. 9l, indicating the square vortex-lattice spin state. Here, tiny orthorhombic distortion of crystal structure (phases I & V) from the high-temperature tetragonal structure (phases VI & IV) has been reported in Ref. [S10]. Considering the anisotropy of the spin textures in Supplementary Fig.   9i-l, the reported orthorhombic structural distortion probably reflects the symmetry change of the spin texture.

IX. Theoretical simulation of spin configuration
To investigate the microscopic origins of the experimentally observed square and rhombic skyrmion lattice (SkL), we performed simulated annealing for the two-dimensional square lattice system based on the effective spin Hamiltonian derived from the Kondo lattice model [S11], which is given by . We set ̂= 1 as the energy unit of the model, and introduced the anisotropic parameters 1 = 0.85 3 , 2 = 0.95 1 3 , 4 = 0.9 3 , 5 = 0.025 3 , 6 = 3 for 3 = 1 and = 0.95 . The parameters 1 and 2 ( 4 and 5 ) stand for the in-plane bond-dependent anisotropy that fix the spiral plane, while 3 ( 6 ) denotes the easy-axis anisotropy that favors the SkL at 1 and 2 ( 1 and 2 ) [S12]. We set the parameters so as to satisfy 1 = 2 > 1 = 2 by taking = 0.95 , which means that the helix with 1 or 2 has smaller energy than that with 1 and 2 . It is noted that the anisotropic parameters 1 -6 are related to the structure of the Fermi surfaces in addition to the spin-orbit coupling, but we here regard them as phenomenological parameters. The last term of Eq. (S7) represents the Zeeman term under an external magnetic field H along the [001] direction. The energy in Eq. (S7) is minimized by performing simulated annealing from high temperature. The simulations are carried out with the standard Metropolis local updates. The results are obtained for systems with N = 100 2 sites under periodic boundary conditions. In each simulation, we independently performed simulated annealing for different H to find low-energy spin configuration by gradually reducing the temperature with the rate +1 = where is the temperature in the nth step. We take = 0.99995 − 0.99999 and the final temperature is typically taken at T = 0.01, which is reached by performing a total of 10 5 -10 6 Monte Carlo sweeps.
Supplementary Fig. 10 shows the magnetic field dependence of stable spin textures theoretically calculated for H || [001] with ̂= 0 and = 0.95. We have mainly identified four distinct magnetic phases (I, II, III, and IV) before reaching the forced ferromagnetic (FM) phase. The simulated real-space distribution of local magnetization ( ) and the corresponding momentum-space distribution of the spin structure factor |̂( )| 2 for each phase are displayed in Supplementary Fig. 10a-e and f-j, respectively. The zero-field state, phase I, is a single-Q state hosting a screw spin texture ( Supplementary Fig. 10a). By increasing H, the double-Q state without four-fold symmetry is stabilized in phase II, which shows the rhombic SkL spin texture ( Supplementary Fig. 10b). In phase III, the spin texture is characterized by the four-fold-symmetric double-Q structure corresponding to the square SkL state ( Supplementary Fig. 10c) While further fine tuning of the magnetic anisotropy and/or value may be required to fully reproduce all the details of observed phase transitions, our main conclusion here, i.e., frustration between the distinct peaks ( 1 / 1~1 ) stabilizes the two distinct SkL states under the easy-axis magnetic anisotropy, will be valid in general.

X. Charge density wave and crystal structural transition
EuAl4 has been reported to host a charge density wave characterized by the incommensurate ordering vector along the [001] axis below 145 K [S10]. Because the sinusoidal wave is generally centrosymmetric, there must be a point where the inversion center of the original crystal structure coincides with that of the sinusoidal wave in case of the incommensurate order. Therefore, the present system can be considered as centrosymmetric even with the incommensurate charge density wave.
In Ref. [S10], tiny orthorhombic distortion of crystal structure from the high-temperature tetragonal phase has also been reported below 12.2 K at 0 T. Since the system exhibits the incommensurate magnetic modulation in both tetragonal (phases I & V) and orthorhombic (phases VI & IV) crystal structures ( Supplementary Fig. 9), the reported small orthorhombic crystal distortion can be excluded as the origin of incommensurate magnetic modulation. axis, as shown in Supplementary Fig. 1b. d-f, The corresponding data for the higher- and d were measured in the field increasing process after zero field cooling at 5.0 K and 11.0 K, respectively. The data for e and f were obtained in the temperature increasing process at 1.0 T and 0.05 T, respectively. g,h, Rhombic and square vortex lattice spin texture described by Eq. (S3), which correspond to the low-T and high-T region of phase IV.