Squeezing the periodicity of Néel-type magnetic modulations by enhanced Dzyaloshinskii-Moriya interaction of 4d electrons

In polar magnets, such as GaV4S8, GaV4Se8 and VOSe2O5, modulated magnetic phases namely the cycloidal and the Néel-type skyrmion lattice states were identified over extended temperature ranges, even down to zero Kelvin. Our combined small-angle neutron scattering and magnetization study shows the robustness of the Néel-type magnetic modulations also against magnetic fields up to 2 T in the polar GaMo4S8. In addition to the large upper critical field, enhanced spin-orbit coupling stabilize cycloidal, Néel skyrmion lattice phases with sub-10 nm periodicity and a peculiar distribution of the magnetic modulation vectors. Moreover, we detected an additional single-q state not observed in any other polar magnets. Thus, our work demonstrates that non-centrosymmetric magnets with 4d and 5d electron systems may give rise to various highly compressed modulated states.

The SkL phase was first observed in the chiral cubic helimagnet MnSi with moderate SOC 24 . In cubic helimagnets the symmetrydictated form of DMI enables Bloch-type magnetic modulations, i.e., spin helices and Bloch-skyrmions. This Bloch-type SkL, comprising q-vectors perpendicular to the direction of the applied field, is stabilized by thermal fluctuations over the energetically favoured longitudinal conical state, only in the close vicinity of the Curie temperature 24,25 . Cubic magnetocrystalline anisotropies that are higher-order terms in the SOC determine the orientation of helical order at zero field and induce small deflections of the SkL planes for fields applied along low-symmetry crystallographic directions 25,26 . The increasing strength of the SOC was studied in Mn(Si 1−x Ge x ) by replacing Si with heavier Ge 27 . As a result, the periodicity of the magnetic modulation decreases and the SkL state is transformed to a hedgehog-lattice state. The enhanced cubic anisotropy can also lead to SkL states at low temperatures, without relying on stabilization by thermal fluctuations 28,29 .
For potential memory and spintronic applications a further reduction in the skyrmion size is desired, which can be achieved through the enhancement of the SOC. In the vast majority of the skyrmion host crystals reported to date, the magnetism is governed by the 3d electrons of transition metals, such as V, Mn, Fe, Co, Cu. Here, we demonstrate the emergence of particularly robust modulated magnetic phases in the 4d polar magnet, GaMo 4 S 8 [30][31][32][33][34] . In the rhombohedral phase, our small-angle neutron scattering (SANS) data (Fig. 1c) show that the strong SOC reduces the periodicity of the magnetic structure to λ ≈ 9.8 nm, being one of the shortest modulation observed in bulk crystals hosting the SkL state due to DMI. Moreover, it modifies the distribution of the q-vectors, as illustrated in Fig. 1d, which we explain by an effective Landau theory containing a cubic anisotropy term in addition to the axial anisotropy, inherent to the rhombohedral state. Our temperature, magnetic field and angular dependent magnetization and electric polarization measurements reveal a complex phase diagram where we attributed phases to the cycloidal and SkL states and observed a third modulated magnetic phase. Moreover, we found that the modulated spin states extend up to fields as high as 2 T, which further support their robustness.
Spin cycloids and Néel SkL have been found recently in the polar phase of the lacunar spinels GaV 4 S 8 and GaV 4 Se 8 [35][36][37][38] , which are narrow-gap multiferroic semiconductors 39,40 . The polar rhombohedral structure (space group R3m) develops via a cooperative Jahn-Teller distortion driven by the unpaired electron of each V 4 cluster with spin S = 1/2 occupying a triply degenerate orbital in the hightemperature cubic state (F43m) 39 . The rhombohedral distortion can occur as elongation of the unit cell along any of the four <111>-type cubic directions, thus, four rhombohedral domain states, that we denote by P 1−4 , are present below the structural transition at~42 K 37,39,41 . (Throughout this paper we use the pseudo-cubic notation.) In contrast to the chiral cubic helimagnets, in these lacunar spinels the cycloidal character of the modulations with q-vectors restricted to the plane perpendicular to the rhombohedral axis (Fig. 1a, b) enhances the stability range of the SkL phase 35,37,42 . In both compounds, SANS experiments revealed 1 magnetic modulations with wavelengths of λ~20 nm [35][36][37]43 , indicating a similar ratio of the exchange interactions and the DMIs.
In the compound GaMo 4 S 8 studied here, the tetrahedral Mo 4 clusters carry an unpaired hole with spin S = 1/2 39 . Correspondingly, the cubic state is rhombohedrally distorted by the compression of the unit cell along one of the <111>-type directions 31 . Although, Rastogi et al., pointed out the importance of the correlation in GaMo 4 S 8 and found a rich magnetic phase diagram below T C = 19 K long ago 30 , the spin ordering patterns of these phases have not been studied.

RESULTS
Wavevector distribution of the magnetic cycloid at zero field We explored the zero-field modulated magnetic states of GaMo 4 S 8 by SANS experiments. (The scattering geometry is shown in Supplementary Fig. 4) The scattered intensity was recorded in zero-field at 2 K upon the 180 ∘ rotation of the sample in 1 ∘ steps, with an acquisition time of 120 s at each angle. The background signal was measured in the paramagnetic phase at 25 K following the same procedure. The scattering images were averaged over a 10 ∘ moving window in the rotation angle to improve the signal-tonoise ratio. Figure 2a-d shows the SANS images obtained on four high-symmetry planes, namely the (111), (110), 112 À Á and (001) planes. A pixel-wise adaptive Wiener filter, assuming Gaussian noise, was applied for better visualization.
The q-dependence of the scattering intensity in the (111) plane, averaged over the polar angles, was fitted by a Gaussian, yielding q j j ¼ 0:64nm −1 for the length of the modulation vectors with a FWHM of 0.2 nm −1 , which corresponds to a real-space periodicity of λ ≈ 9.8 nm. The uncertainty of q j j mainly originates from the broad and anisotropic distribution of the scattering intensity. Whereas the Curie temperature is close to that of GaV 4 Se 8 , suggesting a similar strength of the symmetric exchange (J) in the two compounds, the modulation wavelength in GaMo 4 S 8 is roughly half of that in GaV 4 S 8 35 and GaV 4 Se 8 37 , implying a stronger DMI coupling (D), as λ ∝ J/D.
Over the wide-angle rotation experiment, each scattering image represents a planar cross section of the three-dimensional distribution of the q-vectors, where the azimuthal angle of the vertical slicing plane is varied via stepwise rotation of the sample. The 3D scattering pattern was reconstructed using the whole set of the cross section images. In order to enhance the signal-tonoise ratio and to eliminate the asymmetries of the scattering pattern introduced by imbalances between the populations of the different structural domain states, the 3D scattering pattern was symmetrized for all the symmetry operations of the cubic It is instructive to compare the reciprocal-space q-distributions in GaV 4 S 8 and GaMo 4 S 8 , as shown in Fig. 1a and c, respectively. The neutron scattering data collected in the zero-field cycloidal phase of GaV 4 S 8 at 12 K is reproduced from ref. 36 . In both compounds the cycloidal q-vectors are distributed over four intersecting rings corresponding to the four structural domain states. The ring structure, instead of six well-defined Bragg spots, is due to static orientational disorder of the q-vectors, as discussed in ref. 36 . However, in contrast to GaV 4 S 8 , where the modulation vectors are evenly distributed over rings restricted to the 111 f g-type planes, in GaMo 4 S 8 , the four rings of the q-vectors wave out of the 111 f g planes, crossing them only along the <110>-type directions. This waving pattern of the q-vectors preserves the three-fold rotational symmetry of the rhombohedral structure as highlighted by the schematic image in Fig. 1d.
To explain the distribution of the q-vectors, the following effective Landau potential for the unit vectorq is considered with its x, y, z components defined in the cubic setting, x þq 4 y þq 4 z þ :::: The first term describes the uniaxial anisotropy emerging in the rhombohedral phase with then unit vector parallel to any of the four <111>-type polar axes, and the second term is the lowestorder term compatible with the cubic T d symmetry. The length of the q-vectors is essentially fixed by the DMI, q~D/J, which is consistent with the experimental SANS data. The coefficient α parameterizing the relative strengths of the two terms in Eq. (1) is thus effectively of second order in the SOC. The first term favours the confinement of the q-vectors normal to the polar axes,n, as imposed by the DMI. The waving of the q-vectors out of the {111} planes is captured by the second term. On the microscopic level it represents magnetocrystalline contributions to the Ginzburg-Landau theory for the magnetization that are effectively of fourth order in the SOC.
The minimal-energy solutions to Eq. (1) are sought by parametrizingq in spherical coordinates. The polar angle, Θ is chosen with respect to the polar axisnk<111> of each domain. The azimuthal angle, Φ is enclosed between the in-plane component ofq and one of the corresponding <110> directions. In the limit of weak SOC, α is negligible and the wavevectors are basically in-plane (Θ ≈ π/2), since the spirals are degenerate with respect to the azimuthal angle Φ. This gives rise to an equal distribution on circles, as observed for GaV 4 S 8 in Fig. 1a. Expanding the potential in Θ around π/2 up to the second order and minimizing one obtains for the deviations from the circle that describes the waving of q in GaMo 4 S 8 . The least-square fitting of the data to Eq. (2) yields α = −0.14 ± 0.003 and |q| = 0.64 nm −1 ± 10 −3 (λ = 9.81 nm). As shown in Fig. 2e-h, the fitted model is in excellent agreement with the SANS data, indicating the significant influence of cubic anisotropies in GaMo 4 S 8 as opposed to GaV 4 S 8 in accord with the stronger atomic SOC of Mo. We note that for GaV 4 S 8 we find α = 0 ± 0.04 within the experimental precision. The uncertainty is mainly determined by the image preconditioning rather than the accuracy of the fitting, for details see Supplementary Notes 2 and 3.
The tilting of q out of the {111} planes as well as higher-order SOC represented by additional termsq 6 x þq 6 y þq 6 z in the Landau potential break the degeneracy of q-vectors around the ring and favour either h110i or h112i-type directions. However, within our experimental accuracy, we were not able to resolve an enhanced intensity along any of these directions, see Supplementary Notes 3 and Supplementary Fig. 3.
Magnetic phase diagram at finite magnetic fields We explore the modulated magnetic states stabilized by the interplay of external magnetic fields and the strong SOC of GaMo 4 S 8 . Due to the polar symmetry of this compound, the phase diagram may depend on the angle between the polar axis and the field, thus, we measured longitudinal magnetic and differential magnetoelectric susceptibility at 10 K while the field was rotated in finite steps within the (110) plane between successive field sweeps. The obtained data are respectively shown in Fig. 3a, b as colour maps over the field magnitude-orientation plane, where the angle ϕ was measured from the [111] axis as sketched in Fig.  3c. The critical fields identified as peaks or sharp steps in the susceptibility and magnetocurrent curves are indicated by lines in Fig. 3a, b. (Field scans measured at a few high symmetry directions are shown in Supplementary Figs. 2 and 3.) Compared to the early magnetization study performed on powder samples 30 , in the single crystal sample, we resolved more than two phase transitions strongly depending on the direction of the field. The interpretation of the complicated angulardependent pattern of the anomalies can be simplified by assuming that all four polar domain states are present in the sample. Anomalies and their replica appearing at positions shifted by~109 ∘ are attributed to the P 1 (white lines) and the P 2 (black lines) domains since both polar axes lie within the rotation plane and they span~109 ∘ . The remaining anomaly (grey line) should correspond to both the P 3 and P 4 domains, when the field is rotated in the high-symmetry (110) plane. The field direction [001] is even more special as the polar axes span the same 55 ∘ with the field in all domains, thus, the lines should intersect each other. Such an intersection occurs at~1.5 T where the small deviations likely caused by a slight misorientation of the sample. However, in the same angular range phase boundaries are split around 0.75 T indicated by dotted lines in Fig. 3a. We found that anomalies observed on the different domains collapse into a common phase diagram when plotted on the H ∥ -H ⊥ plane, where H ∥ and H ⊥ represent the field component parallel and perpendicular to the polar axis, respectively (Fig. 3d). This confirms that the sample is in a polar multidomain state and the magnetic state stabilized by the field depends only on the angle between the polar axis and the field. We note that the finite magnetocurrent signal detected in all phases implies that the modulated magnetic states couple to the electric polarization as observed in the sister compounds GaV 4 S 8 and GaV 4 Se 8 38,41,44 . This magnetoelectric coupling is allowed by the polar symmetry of GaMo 4 S 8 41 , and the magnetic field-induced changes of the spin texture modify the polarization texture.
In order to elucidate the nature of the various phases, we explored their properties for certain directions of the applied field in more details. We studied the magnetic field-induced changes in spin textures at 10 K by SANS experiments performed for fields parallel to the [111] axis on the (111) scattering plane. A representative set of SANS images are displayed in Fig. 4a. The smeared wavy ring of intensity present in low magnetic fields, see Fig. 2a, continuously evolves to well-defined Bragg spots in moderate fields pointing to a field-driven enhancement of the magnetic correlation length. To make our analysis more quantitative the field dependence of the scattering intensity and the average length of the q-vectors are plotted in Fig. 4c, d, respectively, in comparison with the susceptibility shown in Fig.  4b. The grey and red curves are obtained by fitting a Gaussian peak on the azimuthally averaged scattering intensity around the h110i and h112i-type directions highlighted by grey and red sectors in Fig. 4a.
The phase transitions are marked by clear anomalies of the scattering parameters. As we have demonstrated above, only the P 1 domain, responsible for the red ring in Figs. 1d and 2, scatters neutrons into the red sections whereas all domains contribute to the intensity detected in the grey regions. Remarkably, our measurements evidence that modulated magnetic structures are extremely robust against a magnetic field, i.e., they extend up to 1.3 and 1.8 T for fields parallel to and inclined at 71 ∘ from the polar axis, respectively. Such exceptional stability of modulated bulk phases has been observed so far only in centrosymmetric materials, hosting nearly atomic-scale skyrmions due to exchange frustration 45,46 . As the modulated phases are expected to be stable up to a critical field of order H FM ∝ D 2 /J, this finding corroborates that the DMI is the strongest in GaMo 4 S 8 among the lacunar spinels known to host SkLs, likely due to enhancement of the SOC from 3d to 4d electrons. Interestingly, after a low-field decrease |q| increases in higher fields in the grey region, which is unusual among skyrmion host spiral magnets. Since |q| does not change above 1.3 T, where only the P 2−4 domains contribute to the SANS intensity, the anomalous field induced shortening of the magnetic periodicity occurs in the P 1 domain, i.e., for field parallel to the polar axis.
Since the DMI pattern dictated by the C 3v symmetry induces cycloidal modulation 9 we assign the zero field ground state of GaMo 4 S 8 to a single-q cycloidal state (orange region in Fig. 3d) in analogy with GaV 4 S 8 and GaV 4 Se 8 . Although it is well established both theoretically and experimentally that the cycloidal phase is the zero field state in polar magnets of C nv symmetry, future polarized neutron diffraction experiments may be performed to directly prove this. The hexagonal SANS pattern with minimal intensity in the red sections, that is observed above~0.75 T, imply the formation of a SkL state or can alternatively manifest the coexistence of cycloidal domains with q-vectors pointing to the different h110i-type directions. Following the correspondence with the phase diagram of the sister compounds we rather attributed the purple phase pocket in Fig. 3d to the SkL state. The robustness of this phase against strong tilting of the magnetic field from the polar axes implies an easy-axis type magnetic anisotropy, thus, it has the same character as in GaV 4 S 8 42,47 , which is in agreement with the predictions of recent first-principles calculations 48,49 .
At finite fields between the cycloidal and the field-polarized FM states for fields perpendicular to the polar axis we detected an additional phase (coloured in blue in Fig. 3d) that is not present in the former material-specific model calculations 42,[48][49][50][51] . We labelled it as 1q state since it is a modulated phase according to SANS,    however, we cannot unambiguously determine its internal spin structure. As it possesses a finite magnetization it is likely to be a fan or a conical phase, which is not present in any other lacunar spinels. The stability range of this intermediate state narrows as the field is rotated toward the polar axis, however, it can separate the cycloidal and SkL phases even in parallel fields. The anomalous increase in |q| above 0.5 T also corresponds to the emergence of this single-q state. Although in the sister compounds the emergence of magnetic states at the polar domain walls (DWs) has been reported 44,52 , we rather assigned the 1q state to a bulk phase. At the onset of the 1q state, we found robust anomalies in all experiments including SANS. Since SANS does not possess the sensitivity to detect DW-confined states with small volume fraction, as found in GaV 4 Se 8 44 , we believe the 1q state corresponds to a new bulk phase.
The temperature-field phase diagram is mapped by susceptibility measurements below T C = 19 K as shown in Fig. 5a Fig. 1.) Schematic phase diagrams deduced for a polar monodomain sample are presented in Fig. 5c and d for fields respectively spanning zero and 55 ∘ with the polar axis. The latter case, H∥[001] likely provides a representative phase diagram for the whole angular range 55 ∘ -90 ∘ as evidenced for 10 K in Fig. 3d. All phases observed at 10 K including the SkL state extend down to the lowest temperatures. In case of strong easy-axis anisotropy, modulated phases are stabilized only by thermal fluctuations as shown for GaV 4 S 8 36 , thus, the presence of cycloidal and SkL phases at low temperatures implies weaker easy-axis anisotropy in GaMo 4 S 8 . The single-q state enters between the cycloidal and the SkL phases only below 13 K for fields applied along the polar axis, and its field stability range grows toward low temperatures. Moreover, below 6 K additional anomalies of the susceptibility appear suggesting the emergence of a new phase not present at 10 K in Fig. 3.

DISCUSSION
The main features of the phase diagram and their assignment are consistent with recent theoretical studies of GaMo 4 S 8 that combines DFT calculations and Monte-Carlo simulations [48][49][50] . These works predict only two modulated structures, the cycloidal and the Néel-type SkL states for field parallel to the polar axis. With a moderate strength of uniaxial anisotropy, these states were found to be stable down to the lowest temperatures. Although the absolute values of the interaction strength largely vary between the different calculations all of them predict a critical temperature consistent with the experiments. The ratio of the relevant DMI and the exchange interactions determining the periodicity of the cycloidal order is similar in refs. 48,50 , and the calculated periodicity (about λ~10 nm) is close to the experimentally observed one. According to Nikolaev and Solovyev 48 , the small periodicity of the cycloidal and SkL phases is not only the consequence of the enhanced DMI due to stronger SOC, but the isotropic exchange term is also reduced from GaV 4 S 8 to GaMo 4 S 8 due to weaker screening. However, further theoretical efforts are required to understand the emergence of the 1q state between the cycloidal and the SkL phases as well as the additional low-temperature phases. An important direction of future research is, in particular, the theoretical investigation of a micromagnetic model that is able to account both for the waving of the wavevector captured by Eq. (1) as well as the magnetic phase diagram observed experimentally.
In conclusion, we studied the magnetically ordered phases of a 4d cluster magnet GaMo 4 S 8 by SANS, magnetization and magneto-current measurements. We found modulated magnetic states with sub-10 nm periodicity that can be attributed to the stronger DMI due to the enhanced SOC of GaMo 4 S 8 , with respect to typical 3d transition metal-based skyrmion hosts. The q-space distribution of the modulation vectors is markedly deformed, which is explained in terms of higher-order anisotropies becoming important in this 4d compound. In finite fields, a series of phase transitions are observed, which is assigned to the transformation of the cycloidal state to the SkL and a new single-q state. Moreover, these modulated spin textures are coupled to the ferroelectric polarization as evidenced by our magneto-current measurements. The exceptional stability of the modulated states against magnetic fields also indicates the importance of SOC in GaMo 4 S 8 . Our findings imply that a remarkable scaling down of the skyrmion size in bulk noncentrosymmetric materials can be achieved by exploring the plethora of 4d and 5d magnets.

METHODS Synthesis
Single crystals of GaMo 4 S 8 (typical size 0.5-1 mm) were grown by the flux method 53 . A mixture of Mo, MoS 2 and Ga 2 S 3 with the molar ratio of 11:3:4 was pressed into a pellet, loaded in an alumina crucible, sealed in a molybdenum tube by arc melting under argon atmosphere, and heated quickly up to 1873 K and cooled slowly to 1273 K.

Small-angle neutron scattering
SANS experiments were performed at the Oak-Ridge National Laboratory High-Flux Isotope Reactor, using the General-Purpose Small-Angle Neutron Scattering Diffractometer 54,55 . A single crystal with a mass of m = 112 mg was mounted onto a rotatable sample stick with its 110 Â Ã cubic direction parallel to the rotation axis. A neutron wavelength of λ n = 6 Å with Δλ n / λ n = 0.13 broadening was used with the detector set to a distance of 5 m from the sample, employing a collimator of the same length.

Magnetization and polarization measurements
The magnetization measurements were performed using a Quantum Design MPMS. The longitudinal magnetic susceptibility calculated from static magnetization measurements as well as the magnetocurrent were measured at 10 K while the field was rotated in fine steps within the (110) plane between successive field sweeps. The step size of the rotation was 1 and 2.5 degrees in the case of the magnetocurrent and the magnetization measurements, respectively. Magnetocurrent (j = dP/dt) measurements were carried out using a Keysight Electrometer. The sample was contacted on parallel (111) faces, and correspondingly, the magnetic field-induced changes in the electric polarization component parallel to the [111] direction was detected. The magnetoelectric susceptibility was calculated from the magnetocurrent by division through the constant magnetic field sweep rate of 1.2 T/min. The field dependence of the magnetization was measured in decreasing temperature steps.

DATA AVAILABILITY
Data associated with figures are provided with this paper. The data and the data evaluation code used in the tomographic analysis of the zero-field SANS data is publicly available at GitHub: https://github.com/Buadam/SANS_Reciprocal_Space_Tomography. All other data that support the plots within this paper are available from the corresponding authors upon reasonable request.